VisRec.m

% Class to handle records from binary .vis files, as generated by the GPU correlator.
% pep/03Feb16

classdef VisRec < handle 
	% Class to abstract access to visibilities generated by the GPU correlator.
	% Supports the different polarizations, channels, antenna combinations
	% possible, represents a single timeslice of data.
	properties
		nchan       = 0;  % Number of channels in this record
		npol		= 0;  % Number of polarizations
		nelem	    = 0;  % Number of antenna elements
		nbline      = 0;  % Number of baselines
		freq		= 0;  % Frequency of the subband of this visibility set.
        freqflag    = 0;  % Bool to indicate whether channel axis flagging should be carried out.
        timeflag    = 0;  % Bool to indicate whether temporal flagging should be carried out.
		fname		= {}; % Filename of file containing data
		fid			= 0;  % fid of file.
		currec		= 0;  % Current record number in the binary file.
		tfilestart  = 0;  % Timestamp of first record in binary file. -- All ts in MJDsec --
		tfileend	= 0;  % Timestamp of last record in binary file.
		trecstart   = 0;  % Timestamp of start of currently read record.
		trecend	    = 0;  % Timestamp of end of currently read record.
		skip		= 0;  % Number of records to skip while reading data.
		deb			= 0;  % Debug level, for verbosity of messages.
        recbytesize = 0;
        datfloatsize= 0;
        flagsig     = 2.5;% Sigma threshold at which to clip
        xx          = [];
        xy          = [];
        yx          = [];
        yy          = [];
        acm_xx      = [];
        acm_yy      = [];
        dt          = 0;

	end % End of properties.

	methods 

		% Constructor
		% Opens the binary file, initializes internal variables.
		% If metainfo. cannot be determined from the data stream,
		% uses the externally supplied values.
        % Arguments:
        %   fname   : Filename of visibility binary file
        %   info    : Structure containing observational information.
		function  obj = VisRec (fname, info)
			assert (exist (fname) == 2);
			obj.fname = fname;
			try
			    obj.fid  = fopen (obj.fname, 'rb');	
			catch ME
		     	fprintf (2, '### Error in opening file %s.\n', obj.fname);
			end;

			if (isfield (info, 'nchan' ) == 0) obj.nchan =  63; else obj.nchan = info.nchan; end;
			if (isfield (info, 'npol'  ) == 0) obj.npol  =   2; else obj.npol = info.npol; end;
			if (isfield (info, 'nelem' ) == 0) obj.nelem = 288; else obj.nelem = info.nelem; end;
            obj.nbline = obj.nelem * (obj.nelem+1)/2;
			if (isfield (info, 'freqflag') == 0) obj.freqflag= 0; else obj.freqflag = info.freqflag; end;
			if (isfield (info, 'timeflag') == 0) obj.timeflag= 0; else obj.timeflag = info.timeflag; end;
			if (isfield (info, 'freq') == 0) obj.freq = 0; else obj.freq = info.freq; end;
			if (isfield (info, 'skip'  ) == 0) obj.skip  =   0; else obj.skip = info.skip; end;
			if (isfield (info, 'deb'   ) == 0) obj.deb   =   0; else obj.deb = info.deb; end;
						
			% obj.recdat = RecDat(); % Self contained data structure to hold a single records' data.
			% obj.getVisMeta (info); % Initialize internal data based on input data.
		end;


		% Destructor
		function delete (obj)
			if (obj.fid > 0)
				fclose (obj.fid);
			end;
		end;


        % Return the internal state of the VisRec object.
        function [obs] = getObsInfo (obj)
		    if (isfield (obj, 'nchan' ) == 0)    obs.nchan     = obj.nchan; end;    
		    if (isfield (obj, 'npol'  ) == 0)    obs.npol      = obj.npol; end;     
		    if (isfield (obj, 'nelem' ) == 0)    obs.nelem     = obj.nelem; end;    
		    if (isfield (obj, 'nbline') == 0)    obs.nbline    = obj.nbline; end;   
		    if (isfield (obj, 'freqflag')== 0)   obs.freqflag  = obj.freqflag; end; 
		    if (isfield (obj, 'timeflag')== 0)   obs.timeflag  = obj.timeflag; end; 
		    if (isfield (obj, 'freq')   == 0)    obs.freq      = obj.freq; end; 
		    if (isfield (obj, 'skip')   == 0)    obs.skip      = obj.skip; end; 
        end;


		% Examine data stream, generate metadata structure.
		% Called only once to initialize internal data.
		% Arguments:
		%	info	: Structure containing the preinitalized parameters of this
		%			  observation. Note: raw visibility header contains only the
		%			  magic, start and end times, nothing more!
		function info = getVisMeta (obj, info)
			fprintf (1, '<-- Generating observation scan visibility meta data...\n');
			assert (obj.fid > 0);
			hdr = fread (obj.fid, 64, 'double');
			meta = obj.interpretHdr (hdr);
            fprintf (1, '%d', meta.magic);
			if (meta.magic == 999878658) % = 0x3B98F002, correlator magic number.
				fprintf (1, '<-- Found Magic number 0x%x in first record.\n',meta.magic);		
                fprintf (1, '<-- File contains RAW Correlations.\n');
				obj.tfilestart = meta.tstart;
                obj.trecstart = meta.tstart;
				fprintf (1, '<-- First record start time: %s.\n', datestr (mjdsec2datenum(obj.tfilestart)));
                if (isfield (info, 'nchan')) obj.nchan = info.nchan; end;
				if (isfield (info, 'nelem')) obj.nelem = info.nelem; end;
				if (isfield (info, 'npol')) obj.npol  = info.npol; end;
				if (isfield (info, 'freq')) obj.freq  = info.freq; end;
				if (isfield (info, 'sub')) obj.freq  = info.sub*195312.5; end;
                if (isfield (info, 'nbline')) obj.nbline= info.nelem*(info.nelem+1)/2; end;
				% In units of floats, the native output datatype of the
				% correlator.
				obj.datfloatsize = 2*obj.npol*obj.nchan*obj.nbline;
				obj.recbytesize = 512 + obj.datfloatsize*4; % Note: In bytes

            % 0x4141525446414143, calibrated visibilities magic number.
            elseif (meta.magic == 1094799956) % Looking at first 32 bits only
            % elseif (meta.magic == 0x41415254 0x46414143)
                fprintf (1, '<-- Found Magic number 0x%x in first record.\n',meta.magic);
                fprintf (1, '<-- File contains CALIBRATED Correlations.\n');
                obj.tfilestart = meta.tstart;
                obj.trecstart = meta.tstart;
                fprintf (1, '<-- First record start time: %s.\n', datestr(mjdsec2datenum(obj.tfilestart)));
                if (isfield (info, 'nchan')) obj.nchan = info.nchan; end;
                if (isfield (info, 'nelem')) obj.nelem = info.nelem; end;
                if (isfield (info, 'npol')) obj.npol  = info.npol; end;
                if (isfield (info, 'freq')) obj.freq  = info.freq; end;
                if (isfield (info, 'sub')) obj.freq  = info.sub*195312.5; end;
                if (isfield (info, 'nbline'))
                    obj.nbline= info.nelem*info.nelem;
                end;
                % In units of floats, the native output datatype of the
                % correlator.
                obj.datfloatsize = 2*obj.npol*obj.nchan*obj.nbline;
                obj.recbytesize = 512 + obj.datfloatsize*4; % Note: In bytes
			end;
		
            fprintf (1, '<-- Record size : %d Bytes, datasize: %d floats.\n', obj.recbytesize, obj.datfloatsize);

            % Move to next record to obtain integration time
			fseek (obj.fid, 0, 'bof'); % Move to the beginning of the file for further operations.
            fseek (obj.fid, obj.recbytesize, 'cof'); % Move past first record, which can have a 0 in its timestamp.
            
    		hdr = fread (obj.fid, obj.recbytesize/8, 'double'); % Read in the full record as doubles.
    		meta = obj.interpretHdr (hdr);
            obj.tfilestart = meta.tstart; % Rejecting first record.
            obj.trecstart = meta.tstart;
            prevtime = obj.tfilestart;
            for i = 1:3 % Arbitrarily reading 3 records to determine dt.
    			hdr = fread (obj.fid, obj.recbytesize/8, 'double'); % Read in the full record as doubles.
    			meta = obj.interpretHdr (hdr);
                obj.trecstart = meta.tstart;
                dt(i) = meta.tstart - prevtime;
                prevtime = meta.tstart;
    			if (meta.magic == 999878658) % = 0x3B98F002, correlator magic number.
    				fprintf (1, '<-- Found Magic number 0x%x in record %d, time %s, dt %fsec.\n', ...
                             meta.magic, i, datestr (mjdsec2datenum(obj.trecstart)), dt(i));		
                end;
            end;
            obj.dt = mean (dt);
                
            
			% Move to end of the file to obtain last record.
			fseek (obj.fid, -(obj.recbytesize), 'eof');
			hdr = fread (obj.fid, obj.recbytesize/8, 'double'); % Read in hdr of the last record.
			meta = obj.interpretHdr (hdr);
			if (meta.magic == 999878658)
				fprintf (1, '<-- Found Magic number 0x%x in last record. %x\n', meta.magic);		
				obj.tfileend = meta.tstart;
				fprintf (1, '<-- Last record start time: %s.\n', datestr(mjdsec2datenum(obj.tfileend)));		
			end;
			fseek (obj.fid, 0, 'bof'); % Move to the beginning of the file for further operations.
		end;
		

		% Function interprets the binary header and returns meta information.
		% Arguments:
		%	hdr : array of 512 uint8, corresponding to the first 512 bytes of
		%		  the visibility record.
		% Returns :
		%	meta: A structure containing metadata as extracted from the header.
		function meta = interpretHdr (obj, hdr)
			% Set the object magic value
            tmp = typecast (hdr(1), 'uint32'); 
			meta.magic  = tmp(1);
            
            if (meta.magic == 999878658)
                meta.tstart = unixtime2mjdsec (hdr(2));
                meta.tend   = unixtime2mjdsec (hdr(3));
            elseif (meta.magic == 1094799956) % 0x41415254
                meta.tstart = unixtime2mjdsec (hdr(3));
                meta.tend   = unixtime2mjdsec (hdr(4));
            end;
		end;


        % Function to carry out visibility splatting on a common grid, but
        % across frequencies.
        % Arguments:
        %   vis     : Visibilities over different subbands/pols.
        %  freq     : Frequencies of the various subbands.
        %   gparm   : Gridding parameter structure.
        % poslocal  : Locational information on the antennas.
        % Returns:
        %   gvis    : Gridded visibilities.
        function vis = gridFreq (obj, vis, freq, gparm)

        end;


        % Function to carry out visibility based flagging along the freq.
        % axis. Statistics are generated over the channel axis.
        % Arguments:
        %   chans   : The Channels over which to calculate statistics.
        %   dat     : Reshaped float array of visibilities.
        % Returns:
        %   flagdat : visibilities mean after ignoring flagged channels.
        function flagdat = flagFreq (obj, chans, dat)
            flagdat = zeros (2, obj.npol, obj.nbline);

            % Doing it the inefficient way for now.
            for v = 1:obj.nbline
                for p = 1:obj.npol
                    for r = 1:2 % re/im
                        seldat = squeeze (dat (r, p, chans, v));
                        sel = ones (obj.nchan, 1) & ~isnan (seldat);
                        for i = 1:5
			                mvis = mean (seldat(sel == 1));
			                svis = std  (seldat(sel == 1));
			                sel  = (abs (seldat(sel) - mvis) < obj.flagsig*svis);
                        end
                        flagdat (r,p,v) = mean (seldat(sel));
                    end;
                end;
            end;
            
        end;


        % Function to convert data stored in the xx/yy linear arrays into
        % a square symmetric Array Covariance Matrix.
        function dat2acm (obj)
            assert (~isempty (obj.xx));
            assert (~isempty (obj.yy));
            t = triu (ones(obj.nelem));
            obj.acm_xx = zeros (obj.nelem);
            obj.acm_xx (t == 1) = obj.xx; 
            acm_diag = diag (diag (obj.acm_xx)); 
            obj.acm_xx = conj (obj.acm_xx + obj.acm_xx' - acm_diag);

            obj.acm_yy = zeros (obj.nelem);
            obj.acm_yy (t == 1) = obj.yy; 
            acm_diag = diag (diag (obj.acm_yy)); 
            obj.acm_yy = conj (obj.acm_yy + obj.acm_yy' - acm_diag);
        end;

        % Function to skip raw records in the .vis file.
        % Assumes that the observation parameters of data are already
        % initialized.
        % Arguments:
        %   rec : Number of records to skip.
        %   unit: string ('rec', 'dt', 'mjdsec') Whether to skip records or seconds, or to
        %   an absolute time.
        % Returns:
        %   Void.
        function skipRec (obj, rec, unit)
            % If unit is time (sec), determine number of records to skip based
            % on integration time and current time.
            switch lower(unit)
                case 'sec' % Should be in seconds
                    recs2move = int32(rec/obj.dt);
                    fprintf (1, '<-- Moving %d recs for %d seconds.\n', recs2move, rec)
                    stat = fseek (obj.fid, obj.recbytesize*recs2move, 0);
                    if (stat < 0)
                        throw (MException ('VisRec.m:skipRec()', ferror (obj.fid)));
                    end;
                
                case 'mjdsec'
                    assert (obj.tfileend > rec);
                    assert (obj.tfilestart < rec);
                    recs2move = int32((rec-obj.trecstart)/obj.dt);
                    fprintf (1, '<-- Moving %d recs (class %s) for %f seconds.\n', recs2move, class (recs2move), rec - obj.trecstart)
                    stat = fseek (obj.fid, obj.recbytesize*recs2move, 0);
                    if (stat < 0)
                        throw (MException ('VisRec.m:skipRec()', ferror (obj.fid)));
                    end;
                
                case 'rec' 
                    stat = fseek (obj.fid, obj.recbytesize*rec, 0);
                    if (stat < 0)
                        throw (MException ('VisRec.m:skipRec()', ferror (obj.fid)));
                    end;
            
                otherwise
                    fprintf (2, '### Unknown skip option %s.\n', unit);
            end;
        end;

        % Function to write out visibilities as records in the format
        % generated by the C++ pipeline, to be inter-operable.
        %   NOTE : We write out both polarizations, XX followed by YY.
        %          readRec() averages over channels, so we don't allow writing
        %          out individual channels for now, only integrated vis.
        %          This has to be ensured by the caller.
        % Arguments:
        %   vis  : Visibility array to write to disk. Can be raw or calibrated
        %          vis, both are written out as full matrices.
        %   pol  : 0 or 1, indicating vis. is for XX or YY polarity.
        %   fid  : The id of the binary file to which to write. New
        %          records are always appended to the specified file.
        %
        % Returns: Status, true for a successful write.
        %
        % The output header structure matches the one created by the c++
        % pipeline, and replicated below. From
        % aartfaac-calibration-pipeline/src/server/packet.h
        %
        % struct output_header_t
        % {
        %  uint64_t magic;                   ///< magic to determine header
        %  double start_time;                ///< start time (unix)
        %  double end_time;                  ///< end time (unix)
        %  int32_t subband;                  ///< lofar subband
        %  int32_t num_dipoles;              ///< number of dipoles (288 or 576)
        %  int32_t polarization;             ///< XX=0, YY=1
        %  int32_t num_channels;             ///< number of channels (<= 64)
        %  float ateam_flux[5];              ///< Ateam fluxes (CasA, CygA, Tau,
        %  Vir, Sun)
        %  std::bitset<5> ateam;             ///< Ateam active
        %  std::bitset<64> flagged_channels; ///< bitset of flagged channels (8
        %  byte)
        %  std::bitset<576> flagged_dipoles; ///< bitset of flagged dipoles (72
        %  byte)
        %  uint32_t weights[78];             ///< stationweights n*(n+1)/2, n in
        %  {6, 12}
        %  uint8_t pad[48];                  ///< 512 byte block
        % };
        function stat = writeRec (obj,vis, pol, fid)

            % Create metadata hdr in correct format
            % hdr.magic       = uint64 (0x4141525446414143);
            % hdr.magic       = uint64 (0x4141525400000000, 'b');
            hdr.magic       = uint64 (1094799956);
            hdr.start_time  = obj.trecstart;
            hdr.end_time    = obj.trecend;
            hdr.subband     = int32 (obj.freq/195312.5);
            hdr.num_dipoles = int32 (obj.nelem);
            hdr.polarization= int32 (pol); % Write out XX first
            hdr.num_channels= int32 (1);
            hdr.ateam_flux  = single ([0,0,0,0,0]); % Not writing out now.
            hdr.rest        = uint8 (zeros (1, 512-15)); % Set all other fields
                                                         % to 0.

            status = 0;
            try
                fwrite (fid, hdr.magic, '2*uint32');
                fwrite (fid, hdr.start_time, 'double');
                fwrite (fid, hdr.end_time, 'double');
                fwrite (fid, hdr.subband, 'uint32');
                fwrite (fid, hdr.num_dipoles, 'uint32');
                fwrite (fid, hdr.polarization, 'uint32');
                fwrite (fid, hdr.num_channels, 'uint32');
                fwrite (fid, hdr.ateam_flux, '5*single');
                fwrite (fid, hdr.rest, '452*uint8');

                fmt = sprintf ('%d*single', length(vis)*2); % for complex
                fwrite (fid, vis, fmt); % Write out float complex.
            catch ME
                fprintf (2, '### VisRec.m:writeRec(): Error in writing to file id %d.\n', fid);
                status = 1;
            end;

            % return the status
            % return status;
        end;

        % Function to read in a single binary visibility record, and return a
        % subset of channels and pols. Also allows skipping time records.
		% Arguments:
		%	pol  : The pols to  be read in. Bool array, [XX, XY, YX, YY]	
		%	chans: The channel subset to  be read in. Range [1:obj.nchans]	
		% Returns:
		%   flagdat : structure containg the actual data vectors for XX and YY pols, and
		%		  metadata.
		function dat = readRec(obj, pol, chans)
			assert (obj.fid > 0);

            stat = fseek (obj.fid, obj.recbytesize*obj.skip, 0);
            if (stat < 0)
                throw (MException ('VisRec.m:readRec()', ferror (obj.fid)));
            end;

			hdr = fread (obj.fid, 64, 'double'); % Read in the first 512 bytes
            if (feof (obj.fid) > 0)
                throw (MException ('VisRec.m:readRec()', ferror (obj.fid)));
            end;

			% Interpret the header of this data record. Hdr is first 512 bytes.
			meta = obj.interpretHdr (hdr);
            fprintf (1, '<-- Rec time: %s\n', datestr (mjdsec2datenum (meta.tstart)));
			
			assert (length (chans) <= obj.nchan);

			% Update internal state variables of this record
			obj.trecstart = meta.tstart;
			obj.trecend   = meta.tend;

			% dat = fread (obj.fid, obj.datfloatsize, [obj.nbline, obj.nchan, obj.npol, 2], 'single');
			dat = reshape (fread (obj.fid, obj.datfloatsize, 'single'), [2, obj.npol, obj.nchan, obj.nbline]);
            if (feof (obj.fid) > 0)
                throw (MException ('VisRec.m:readRec():', ferror (obj.fid)));
            end;

            if (obj.freqflag == 0)
                flagdat = mean (dat(:,:,chans,:), 3); % Carry out a blind averaging on the channel axis.
            else
                flagdat = obj.flagFreq (chans, dat);
            end;

			if (pol(1)) 
				obj.xx = complex (squeeze(flagdat(1,1,:)), squeeze(flagdat(2,1,:)));
			end;

			if (pol(4)) 
				obj.yy = complex (squeeze(flagdat(1,2,:)), squeeze(flagdat(2,2,:)));
			end;

            if (pol(2) || pol(3))
                fprintf (2, '### Ignoring XY and YX pols not present in incoming data. Ignoring.\n');
            end;

            obj.dat2acm();
		end;


        % Function to read in multiple timeslices of data, sigmaclip in time
        % axis and write the averaged visibility back to the object. For implementing
        % time domain averaging in visibilities before imaging.
        % Calls readRec() underneath.
        % NOTE: The VisRec object no longer corresponds to the raw visibility in
        % the associated file after calling this method! Only the trecstart,
        % freq, nchan, acm_xx and acm_yy are updated. recbytesize, datfloatsize
        % are not updated to allow reading in the next set of raw records.
		% Arguments:
		%	pol  : The pols to  be read in. Bool array, [XX, XY, YX, YY]	
		%	chans: The channel subset to  be read in. Range [1:obj.nchans]	
        %   tavg : The number of seconds (modulo integration time) to average over.
		% Returns:
		%   dat  : Time averaged visibilities, flagged both along the channel and time
		%          axis.
		function dat = readTimeAvgRec (obj, pol, chans, tavg)
            
            % Read in the next record to determine its timestamp.
            obj.readRec (pol, chans);

            % Determine the end timestamp
            tend = obj.trecstart + tavg;

            % Check!
            assert (tend > obj.trecstart);

            % Internally store records
            xxtslices(1, :) = obj.xx;
            yytslices(1, :) = obj.yy;
            trecavg = obj.trecstart - obj.dt + tavg/2;

            if (obj.deb > 2)
                fprintf (2, 'readTimeAvgRec: output time %.2f', dat.trecstart);
            end;

            % Ingest records till avg. is reached.
            ind = 2;
            while (obj.trecstart + obj.dt <= tend)
                obj.readRec (pol, chans);
                xxtslices (ind, :) = obj.xx;
                yytslices (ind, :) = obj.yy;
                ind = ind + 1;
                if (obj.deb > 2)
                    fprintf (2, '%d:%.2f ', ind, obj.trecstart);
                end;
            end;

            if (obj.timeflag)
	            % Selection mask
	            selxx = ones (size (xxtslices));
	            selyy = ones (size (yytslices));
                % ### UNIMPLEMENTED, DONT USE!
	            while 1
	                m = nanmean (xxtslices);
	                s = nanstd  (xxtslices);
	                xxtslices(abs(xxtslices-repmat (m, [obj.nbline, 1])) > repmat (obj.flagsig*abs(s), [obj.nbline, 1])) = nan;
	
	                m = nanmean (yytslices);
	                s = nanstd  (yytslices);
	                yytslices(abs(yytslices-repmat (m, [obj.nbline, 1])) > repmat (obj.flagsig*abs(s), [obj.nbline, 1])) = nan;
	            end;
            else
                obj.xx = mean (xxtslices, 1);
                obj.yy = mean (yytslices, 1);
                obj.trecstart = trecavg;
                % TODO: Check what else needs to be updated in order to keep the
                % the parent VisRec object's metadata consistent.
            end;
        end;

        % Function to generate a calibrated map from the observed
        % visibilities. Returns generated images in the img structure
        % Both XX and YY images are returned.
        %
        % Arguments:
        %   arrayconfig : One of 'lba_outer', 'lba_inner'
        %   flagant     : Index of antennas to be flagged.
        % Returns:
        %   img : structure containing generated images.
        %   sol : Solution structure as generated by the calibration routine.
        function [img, solx, soly] = genCalMap (obj, arrayconfig, flagant)
            load ('poslocal_outer.mat', 'poslocal');
            uloc = meshgrid (poslocal(:,1)) - meshgrid (poslocal (:,1)).';
            vloc = meshgrid (poslocal(:,2)) - meshgrid (poslocal (:,2)).';
            wloc = meshgrid (poslocal(:,3)) - meshgrid (poslocal (:,3)).';
			gparm.type = 'pillbox';
		    gparm.lim  = 0;
		    gparm.duv = 0.5;				% Default, reassigned from freq. of obs. to
											% image just the full Fov (-1<l<1)
		    gparm.Nuv = 512;				% size of gridded visibility matrix
		    gparm.uvpad = 512;				% specifies if any padding needs to be added
		    gparm.fft  = 1;

            img.tobs = obj.trecstart;
            img.freq = obj.freq;


	   	    solx = pelican_sunAteamsub (obj.acm_xx, img.tobs, img.freq, eye(obj.nelem), ... 
				flagant, obj.deb, 1, [], [], 'poslocal_outer.mat', [], []);
	   	    soly = pelican_sunAteamsub (obj.acm_yy, img.tobs, img.freq, eye(obj.nelem), ... 
				flagant, obj.deb, 1, [], [], 'poslocal_outer.mat', [], []);

            [uloc_flag, vloc_flag, ~] = gen_flagged_uvloc (uloc, vloc, wloc, flagant);
   		    [~, img.map_xx, ~, img.l, img.m] = ... 
		        fft_imager_sjw_radec (solx.calvis(:), uloc_flag(:), vloc_flag(:), ... 
					gparm, [], [], img.tobs, img.freq, 0);
   		    [~, img.map_yy, ~, img.l, img.m] = ... 
		        fft_imager_sjw_radec (soly.calvis(:), uloc_flag(:), vloc_flag(:), ... 
					gparm, [], [], img.tobs, img.freq, 0);
        end;


        % Function to generate an uncalibrated map from the observed
        % visibilities. Returns generated images in the img structure
        % Both XX and YY images are returned.
        %
        % Arguments:
        %   arrayconfig : One of 'lba_outer', 'lba_inner'
        % Returns:
        %   img : structure containing generated images.
        function img = genUncalMap (obj, arrayconfig)
            assert (obj.freq ~= 0);
            load ('poslocal_outer.mat', 'poslocal');
            uloc = meshgrid (poslocal(:,1)) - meshgrid (poslocal (:,1)).';
            vloc = meshgrid (poslocal(:,2)) - meshgrid (poslocal (:,2)).';
			gparm.type = 'pillbox';
		    gparm.lim  = 0;
		    gparm.duv = 0.5;				% Default, reassigned from freq. of obs. to
											% image just the full Fov (-1<l<1)
		    gparm.Nuv = 512;				% size of gridded visibility matrix
		    gparm.uvpad = 512;				% specifies if any padding needs to be added
		    gparm.fft  = 1;

            img.tobs = obj.trecstart;
            img.freq = obj.freq;
   		    [~, img.map_xx, ~, img.l, img.m] = ... 
		        fft_imager_sjw_radec (obj.acm_xx(:), uloc(:), vloc(:), ... 
					gparm, [], [], img.tobs, img.freq, 0);
   		    [~, img.map_yy, ~, img.l, img.m] = ... 
		        fft_imager_sjw_radec (obj.acm_yy(:), uloc(:), vloc(:), ... 
					gparm, [], [], img.tobs, img.freq, 0);
        end;
			

	end % End of methods

end % End of classdef