analyzeFullSkyAzElMapSouthFrequency

This function will produce a pixellated map of the sky in azimuth and elevation. Can be used in conjuction with the full_sky_4degscan_el_0.5dstepSouth.sch schedule file.#

0001 function analyzeFullSkyAzElMapSouthFrequency(start,finish,scanLength,dirname,channels)
0002 %This function will produce a pixellated map of the sky in azimuth and elevation. Can be used in conjuction with the full_sky_4degscan_el_0.5dstepSouth.sch schedule file.#
0003 
0004 % start is the start time
0005 % finish is the end time
0006 % scanLength is the amount of data read in per read
0007 % dirname is the location to save the data file
0008 % channels is the number of channels per freuqency bin
0009 % This function analyzes the full-sky az/el maps and plots the data:
0010 % e.g. if the observer log contains:
0011 % 18-Jun-2011:03:11:42      18-Jun-2011:15:10:16 /home/cbass/gcpCbass/control/sch/full_sky_scan_el_0.5dstep.sch
0012 % you should call:
0013 %
0014 %analyzeFullSkyAzElMapSouthFrequency('29-Jan-2014:01:32:13','29-Jan-2014:08:32:13',0.5,'./',4)
0015 % aplotFrequency('AllSkyMap_29-Jan-2014:01:32:13_29-Jan-2014:08:32:13_data.mat',[0 1000],[-500 500])
0016 
0017 %
0018 % Modified by ogk on 6 April 2012 to use the new alpha routines.
0019 % Modified by CJC 10 Sep 2013 for the Southern antenna
0020 % Modified by CJC 10 Jun 2014 for the Southern antenna using frequency data too
0021 
0022 % scanLength is a value in hours
0023 opengl software
0024 
0025 %%%%%%%%%%%%%%%
0026 dateStart = datenum(start,'dd-mmm-yyyy:HH:MM:SS')
0027 dateEnd = datenum(finish,'dd-mmm-yyyy:HH:MM:SS')
0028 %set up the variables
0029 stepSize=channels %umber of bins per freq
0030 freqRes=(500/stepSize); % resolution of the map
0031 numberMaps=64/stepSize;% number of maps that will be produced per Roach
0032 
0033 
0034 if(dateEnd<dateStart)
0035     disp('analyzeFullSkyAzElMapSouthFrequency:: End date before start date- Correct this before running')
0036     return
0037 end
0038 fprefix = sprintf('AllSkyMap_%s_%s_Res_%d_',start,finish,stepSize);
0039 
0040 currentTimeStart = dateStart;
0041 a=figure('visible','off')
0042 b=figure('visible','off')
0043 set(a,'renderer','zbuffer')
0044 set(b,'renderer','zbuffer')
0045 % Now naievely bin the data:
0046 dEl = 0.75;
0047 dAz = 0.75;
0048 azC = 0:dAz:(360-dAz);
0049 elC = 0:dEl:85;
0050 Cmap = zeros(length(elC),length(azC),8);
0051 Nmap = zeros(length(elC),length(azC),1);
0052 %first we instantiate our frequency storage bins
0053 for i=1:numberMaps
0054     T{i}=Cmap;
0055     Tn{i}=Nmap;
0056 end
0057 
0058 while(currentTimeStart <dateEnd)
0059     
0060     currentTimeEnd = currentTimeStart+scanLength/24; %calculate end time
0061    if (currentTimeEnd >= dateEnd)
0062        currentTimeEnd = dateEnd
0063    end
0064        
0065        
0066     currentTimeStartString = datestr(currentTimeStart,'dd-mmm-yyyy:HH:MM:SS') % convert to C-BASS string
0067     currentTimeEndString = datestr(currentTimeEnd,'dd-mmm-yyyy:HH:MM:SS')% convert to C-BASS string
0068     
0069     % read in the data
0070     try
0071         d = read_arcSouth(currentTimeStartString, currentTimeEndString);
0072         %d = pipe_read([],[],d);
0073         %s = 'FILTERED';
0074         %d = assembleAlphaStreams(d,s);
0075         %d = applyAlpha(d,s);
0076         %d = calculateStokes(d,s);
0077     catch
0078         disp('analyzeFullSkyAzElMapSouthFrequency:: Read_arc failed')
0079        currentTimeStart = currentTimeStart+scanLength/24; 
0080        continue
0081     end
0082    % d=rfiRemovalSouth(d);
0083 %I = d.index.source.fast;
0084 I = (1:length(d.antenna0.servo.apparent));
0085 az = d.antenna0.servo.apparent(I,1);
0086 el = d.antenna0.servo.apparent(I,2);
0087     for i=1:numberMaps
0088         Cmap = T{i}; %read in the relevant Cmap
0089         Nmap = Tn{i}; %read in the relevant NMap
0090         startBin=i*(64/numberMaps)-(64/numberMaps)+1; %start of frequency data
0091         freqBins=[startBin:startBin+stepSize-1]; %frequency data bins -- takes the standard 1 indexing=1 of matlab into account
0092         disp(i)
0093         disp(freqBins)
0094         dataIL1 = mean(d.antenna0.roach1.LL(:,freqBins)-d.antenna0.roach1.load2(:,freqBins),2); %%
0095         dataIL2 = mean(d.antenna0.roach2.LL(:,freqBins)-d.antenna0.roach2.load2(:,freqBins),2);
0096         dataIR1 = mean(d.antenna0.roach1.RR(:,freqBins) -d.antenna0.roach1.load1(:,freqBins),2);
0097         dataIR2 = mean(d.antenna0.roach2.RR(:,freqBins) - d.antenna0.roach2.load1(:,freqBins),2);
0098         dataQ1 = mean(d.antenna0.roach1.Q(:,freqBins),2);
0099         dataU1 = mean(d.antenna0.roach1.U(:,freqBins),2);
0100         dataQ2 = mean(d.antenna0.roach2.Q(:,freqBins),2);
0101         dataU2 = mean(d.antenna0.roach2.U(:,freqBins),2);
0102         features=d.array.frame.features;
0103 
0104         [year,month,day] = mjd2date(tstr2mjd(start));
0105         if ~strcmp(dirname(end),'/')
0106             dirname = [dirname '/'];
0107         end
0108 
0109     %fprefix = sprintf('AllSkyMap-%d-%d-%d-',year,month,day);
0110 
0111 
0112 
0113         % Now do this for each and every sample:
0114         Nsamples = size(az,1);
0115         disp('analyzeFullSkyAzElMapSouthFrequency:: Starting map making...')
0116         for k=1:Nsamples
0117             if mod(k,Nsamples/100)==0
0118             %fprintf('%3.0f%% done with binning\n',k/Nsamples*100)
0119             end
0120             % What bin does this sample belong to?
0121             Iaz = abs(azC-az(k))<dAz/2;
0122             Iel = abs(elC-el(k))<dEl/2;
0123             Cmap(Iel,Iaz,1) = Cmap(Iel,Iaz,1)+dataIL1(k,1); %%LCP Stoke I
0124             Cmap(Iel,Iaz,2) = Cmap(Iel,Iaz,2)+dataQ1(k,1);
0125             Cmap(Iel,Iaz,3) = Cmap(Iel,Iaz,3)+dataU1(k,1);
0126             Cmap(Iel,Iaz,4) = Cmap(Iel,Iaz,4)+dataIR1(k,1); %%RCP Stoke I
0127             Cmap(Iel,Iaz,5) = Cmap(Iel,Iaz,5)+dataIL2(k,1); %%LCP Stoke I
0128             Cmap(Iel,Iaz,6) = Cmap(Iel,Iaz,6)+dataQ2(k,1);
0129             Cmap(Iel,Iaz,7) = Cmap(Iel,Iaz,7)+dataU2(k,1);
0130             Cmap(Iel,Iaz,8) = Cmap(Iel,Iaz,8)+dataIR2(k,1); %%RCP Stoke I
0131             %Cmap(Iel,Iaz,2) = Cmap(Iel,Iaz,2)+data(k,6);
0132             %Cmap(Iel,Iaz,3) = Cmap(Iel,Iaz,3)+data(k,7);
0133             %Cmap(Iel,Iaz,4) = Cmap(Iel,Iaz,4)+data(k,8);
0134 
0135             Nmap(Iel,Iaz,1) = Nmap(Iel,Iaz,1)+1;
0136         end
0137         T{i}=Cmap; %store the Cmap to the structure
0138         Tn{i}=Nmap; %store the Nmap to the structure
0139         disp('analyzeFullSkyAzElMapSouthFrequency:: Done map making!')
0140     end
0141     %save([dirname fprefix 'data.mat'],'fprefix','CI1','','','','','','','azC','elC','climI','climP','climV')
0142     %increment start time by scan length so that the next time will be at a
0143         %later time....
0144 currentTimeStart = currentTimeStart+scanLength/24; 
0145 end
0146 %map = Cmap./repmat(Nmap(:,:,1),[1,1,8,numberMaps]);
0147 
0148 
0149 
0150 save([dirname fprefix 'data.mat'],'fprefix','T','Tn','numberMaps','stepSize','freqRes','dEl','dAz','azC','elC','Nmap')
0151     
0152 end