analyzeFullSkyAzElMapSouth

 analyzeFullSkyAzElMapSouth(start,finish,dirname,plotflag)

 This function analyzes the full-sky az/el maps and plots the data:
 e.g. if the observer log contains:
 18-Jun-2011:03:11:42      18-Jun-2011:15:10:16 /home/cbass/gcpCbass/control/sch/full_sky_scan_el_0.5dstep.sch
 you should call:
 
analyzeFullSkyAzElMapSouth_HMH('29-Jan-2014:01:32:13','29-Jan-2014:08:32:13',0.5,'./',1)
 aplot_CJC('AllSkyMap_29-Jan-2014:01:32:13_29-Jan-2014:08:32:13_data.mat',[0 1000],[-500 500])

0001 function analyzeFullSkyAzElMapSouth(start,finish,scanLength,dirname,plotflag)
0002 % analyzeFullSkyAzElMapSouth(start,finish,dirname,plotflag)
0003 %
0004 % This function analyzes the full-sky az/el maps and plots the data:
0005 % e.g. if the observer log contains:
0006 % 18-Jun-2011:03:11:42      18-Jun-2011:15:10:16 /home/cbass/gcpCbass/control/sch/full_sky_scan_el_0.5dstep.sch
0007 % you should call:
0008 %
0009 %analyzeFullSkyAzElMapSouth_HMH('29-Jan-2014:01:32:13','29-Jan-2014:08:32:13',0.5,'./',1)
0010 % aplot_CJC('AllSkyMap_29-Jan-2014:01:32:13_29-Jan-2014:08:32:13_data.mat',[0 1000],[-500 500])
0011 
0012 %
0013 % Modified by ogk on 6 April 2012 to use the new alpha routines.
0014 % Modified by CJC 10 Sep 2013 for the Southern antenna
0015 
0016 % scanLength is a value in hours
0017 opengl software
0018 
0019 %%%%%%%%%%%%%%%
0020 dateStart = datenum(start,'dd-mmm-yyyy:HH:MM:SS')
0021 dateEnd = datenum(finish,'dd-mmm-yyyy:HH:MM:SS')
0022 
0023 if(dateEnd<dateStart)
0024     disp('analyzeFullSkyAzElMapSouth:: End date before start date- Correct this before running')
0025     return
0026 end
0027 fprefix = sprintf('AllSkyMap_%s_%s_',start,finish);
0028 
0029 currentTimeStart = dateStart;
0030 a=figure('visible','off')
0031 b=figure('visible','off')
0032 set(a,'renderer','zbuffer')
0033 set(b,'renderer','zbuffer')
0034 % Now naievely bin the data:
0035 dEl = 0.75;
0036 dAz = 0.75;
0037 azC = 0:dAz:(360-dAz);
0038 elC = 0:dEl:85;
0039 Cmap = zeros(length(elC),length(azC),8);
0040 Nmap = zeros(length(elC),length(azC),1);
0041 
0042 while(currentTimeStart <dateEnd)
0043     
0044     currentTimeEnd = currentTimeStart+scanLength/24; %calculate end time
0045    if (currentTimeEnd >= dateEnd)
0046        currentTimeEnd = dateEnd
0047    end
0048        
0049        
0050     currentTimeStartString = datestr(currentTimeStart,'dd-mmm-yyyy:HH:MM:SS') % convert to C-BASS string
0051     currentTimeEndString = datestr(currentTimeEnd,'dd-mmm-yyyy:HH:MM:SS')% convert to C-BASS string
0052     
0053     % read in the data
0054 try
0055     d = read_arcSouth(currentTimeStartString, currentTimeEndString);
0056     %d = pipe_read([],[],d);
0057     %s = 'FILTERED';
0058     %d = assembleAlphaStreams(d,s);
0059     %d = applyAlpha(d,s);
0060     %d = calculateStokes(d,s);
0061 catch
0062     disp('analyzeFullSkyAzElMapSouth:: Read_arc failed')
0063    currentTimeStart = currentTimeStart+scanLength/24; 
0064    continue
0065 end
0066     d=rfiRemovalSouth(d);
0067 %I = d.index.source.fast;
0068 I = (1:length(d.antenna0.servo.apparent));
0069 az = d.antenna0.servo.apparent(I,1);
0070 el = d.antenna0.servo.apparent(I,2);
0071 % ra = d.antenna0.servo.equa(I,1)*180/pi;
0072 % dc = d.antenna0.servo.equa(I,2)*180/pi;
0073 %dataIL1 = mean(d.antenna0.roach1.LL,2)-mean(d.antenna0.roach1.load2,2);
0074 %dataIL2 = mean(d.antenna0.roach2.LL,2)-mean(d.antenna0.roach2.load2,2);
0075 %dataIR1 = mean(d.antenna0.roach1.RR,2)-mean(d.antenna0.roach1.load1,2);
0076 %dataIR2 = mean(d.antenna0.roach2.RR,2)-mean(d.antenna0.roach2.load1,2);
0077 
0078 dataIL1 = mean(d.antenna0.roach1.LL-d.antenna0.roach1.load2,2);
0079 dataIL2 = mean(d.antenna0.roach2.LL-d.antenna0.roach2.load2,2);
0080 dataIR1 = mean(d.antenna0.roach1.RR -d.antenna0.roach1.load1,2);
0081 dataIR2 = mean(d.antenna0.roach2.RR - d.antenna0.roach2.load1,2);
0082 
0083 %d=convertToKelvin(d);
0084 %dataI = d.antenna0.roach1.LL(:,33); % this looks at data from 1 channel only
0085 %dataI=d.antenna0.roach1.StokesI;
0086 dataQ1 = mean(d.antenna0.roach1.Q,2);
0087 dataU1 = mean(d.antenna0.roach1.U,2);
0088 dataQ2 = mean(d.antenna0.roach2.Q,2);
0089 dataU2 = mean(d.antenna0.roach2.U,2);
0090 features=d.array.frame.features;
0091 
0092 [year,month,day] = mjd2date(tstr2mjd(start));
0093 if ~strcmp(dirname(end),'/')
0094     dirname = [dirname '/'];
0095 end
0096 
0097 %fprefix = sprintf('AllSkyMap-%d-%d-%d-',year,month,day);
0098 
0099 
0100 
0101 % Now do this for each and every sample:
0102 Nsamples = size(az,1);
0103 disp('analyzeFullSkyAzElMapSouth:: Starting map making...')
0104 for k=1:Nsamples
0105     if mod(k,Nsamples/100)==0
0106         %fprintf('%3.0f%% done with binning\n',k/Nsamples*100)
0107     end
0108     % What bin does this sample belong to?
0109     Iaz = abs(azC-az(k))<dAz/2;
0110     Iel = abs(elC-el(k))<dEl/2;
0111     Cmap(Iel,Iaz,1) = Cmap(Iel,Iaz,1)+dataIL1(k,1); %%LCP Stoke I
0112     Cmap(Iel,Iaz,2) = Cmap(Iel,Iaz,2)+dataQ1(k,1);
0113     Cmap(Iel,Iaz,3) = Cmap(Iel,Iaz,3)+dataU1(k,1);
0114     Cmap(Iel,Iaz,4) = Cmap(Iel,Iaz,4)+dataIR1(k,1); %%RCP Stoke I
0115     Cmap(Iel,Iaz,5) = Cmap(Iel,Iaz,5)+dataIL2(k,1); %%LCP Stoke I
0116     Cmap(Iel,Iaz,6) = Cmap(Iel,Iaz,6)+dataQ2(k,1);
0117     Cmap(Iel,Iaz,7) = Cmap(Iel,Iaz,7)+dataU2(k,1);
0118     Cmap(Iel,Iaz,8) = Cmap(Iel,Iaz,8)+dataIR2(k,1); %%RCP Stoke I
0119     %Cmap(Iel,Iaz,2) = Cmap(Iel,Iaz,2)+data(k,6);
0120     %Cmap(Iel,Iaz,3) = Cmap(Iel,Iaz,3)+data(k,7);
0121     %Cmap(Iel,Iaz,4) = Cmap(Iel,Iaz,4)+data(k,8);
0122 
0123     Nmap(Iel,Iaz,1) = Nmap(Iel,Iaz,1)+1;
0124 end
0125 disp('analyzeFullSkyAzElMapSouth:: Done map making!')
0126 
0127 %save([dirname fprefix 'data.mat'],'fprefix','CI1','','','','','','','azC','elC','climI','climP','climV')
0128 %increment start time by scan length so that the next time will be at a
0129     %later time....
0130     currentTimeStart = currentTimeStart+scanLength/24; 
0131 end
0132 Cmap = Cmap./repmat(Nmap(:,:,1),[1,1,8]);
0133 
0134 % the limits on the color bars
0135 climI = [0.5e4 2e4];
0136 climP = [-0.05 0.05];
0137 climV = [-0.05 0.05];
0138 
0139 CILL1 = Cmap(:,:,1); %%LL
0140 CIRR1 = Cmap(:,:,4); %%RR
0141 CV1 = (CILL1-CIRR1)/2; 
0142 CQ1 = Cmap(:,:,2);
0143 CU1 = Cmap(:,:,3);
0144 CP1 = sqrt(CQ1.^2+CU1.^2);
0145 CI1 = (CILL1+CIRR1)/2;
0146 
0147 CILL2 = Cmap(:,:,5);
0148 CIRR2 = Cmap(:,:,8);
0149 CV2 = (CILL2-CIRR2)/2;
0150 CQ2 = Cmap(:,:,6);
0151 CU2 = Cmap(:,:,7);
0152 CP2 = sqrt(CQ2.^2+CU2.^2);
0153 CI2 = (CILL2+CIRR2)/2;
0154 
0155 
0156 save([dirname fprefix 'data.mat'],'fprefix','CILL1','CIRR1','CV1','CQ1','CU1','CP1','CI1','CILL2','CIRR2','CV2','CQ2','CU2','CP2','CI2','azC','elC','climI','climP','climV')
0157 if plotflag
0158     
0159 %     figure('Position',[1 1000 360*3 90*3])
0160 %    % imagesc(azC(:),elC(:),(CI1-repmat(nanmedian(CI1,2),1,size(CI1,2))),climI)
0161 %     imagesc(azC(:),elC(:),(CI1-repmat(nanmedian(CI1,2),1,size(CI1,2))),climI)
0162 %     set(gca,'YDir','normal')
0163 %     colorbar
0164 %     xlabel('Azimuth')
0165 %     ylabel('Elevation')
0166 %     title('I1')
0167 %     set(gcf,'PaperPositionMode','auto')
0168 %     eval(sprintf('print -depsc2 -r200 %s',[dirname fprefix 'I1.eps']))
0169 %
0170 %
0171 %     figure('Position',[1 1000 360*3 90*3])
0172 %    % imagesc(azC(:),elC(:),(CI2-repmat(nanmedian(CI2,2),1,size(CI2,2))),climI)
0173 %     imagesc(azC(:),elC(:),(CI2-repmat(nanmedian(CI2,2),1,size(CI2,2))),climI)
0174 %     set(gca,'YDir','normal')
0175 %     colorbar
0176 %     xlabel('Azimuth')
0177 %     ylabel('Elevation')
0178 %     title('I2')
0179 %     set(gcf,'PaperPositionMode','auto')
0180 %     eval(sprintf('print -depsc2 -r200 %s',[dirname fprefix 'I2.eps']))
0181 %
0182      figure('Position',[1 1000 360*3 90*3])
0183      %imagesc(azC(:),elC(:),(CI-repmat(nanmedian(CI,2),1,size(CI,2))),climI)
0184      imagesc(azC(:),elC(:),(CI1-repmat(nanmedian(CI1,2),1,size(CI1,2))),[0 1000])
0185      set(gca,'YDir','normal')
0186      colorbar
0187      xlabel('Azimuth')
0188      ylabel('Elevation')
0189      title('I=(I1+I2)/2')
0190      set(gcf,'PaperPositionMode','auto')
0191      eval(sprintf('print -depsc2 -r200 %s',[dirname fprefix 'I.eps']))
0192      
0193       figure('Position',[1 1000 360*3 90*3])
0194      %imagesc(azC(:),elC(:),(CI-repmat(nanmedian(CI,2),1,size(CI,2))),climI)
0195      imagesc(azC(:),elC(:),(CI2-repmat(nanmedian(CI2,2),1,size(CI2,2))),[0 1000])
0196      set(gca,'YDir','normal')
0197      colorbar
0198      xlabel('Azimuth')
0199      ylabel('Elevation')
0200      title('I=(I1+I2)/2')
0201      set(gcf,'PaperPositionMode','auto')
0202      eval(sprintf('print -depsc2 -r200 %s',[dirname fprefix 'I2.eps']))
0203 % %
0204 %     figure('Position',[1 1000 360*3 90*3])
0205 %     imagesc(azC(:),elC(:),(CV-repmat(nanmedian(CV,2),1,size(CV,2))),climV)
0206 %     set(gca,'YDir','normal')
0207 %     colorbar
0208 %     xlabel('Azimuth')
0209 %     ylabel('Elevation')
0210 %     title('V=(I1-I2/2)')
0211 %     set(gcf,'PaperPositionMode','auto')
0212 %     eval(sprintf('print -depsc2 -r200 %s',[dirname fprefix 'V.eps']))
0213  
0214 %    figure('Position',[1 1000 360*3 90*3])
0215 %    imagesc(azC(:),elC(:),(CQ-repmat(nanmedian(CQ,2),1,size(CQ,2))),climP)
0216 %     set(gca,'YDir','normal')
0217 %     colorbar
0218 %     xlabel('Azimuth')
0219 %     ylabel('Elevation')
0220 %     title('Q')
0221 %     set(gcf,'PaperPositionMode','auto')
0222 %     eval(sprintf('print -depsc2 -r200 %s',[dirname fprefix 'Q.eps']))
0223  
0224 %     figure('Position',[1 1000 360*3 90*3])
0225 %     imagesc(azC(:),elC(:),(CU-repmat(nanmedian(CU,2),1,size(CU,2))),climP)
0226 %     set(gca,'YDir','normal')
0227 %     colorbar
0228 %     xlabel('Azimuth')
0229 %     ylabel('Elevation')
0230 %     title('U')
0231 %     set(gcf,'PaperPositionMode','auto')
0232 %     eval(sprintf('print -depsc2 -r200 %s',[dirname fprefix 'U.eps']))
0233 %
0234      figure('Position',[1 1000 360*3 90*3])
0235      imagesc(azC(:),elC(:),(CP1-repmat(nanmedian(CP1,2),1,size(CP1,2))))
0236      set(gca,'YDir','normal')
0237      colorbar
0238      xlabel('Azimuth')
0239      ylabel('Elevation')
0240      title('P = sqrt(Q^2+U^2)')
0241      set(gcf,'PaperPositionMode','auto')
0242      eval(sprintf('print -depsc2 -r200 %s',[dirname fprefix 'P.eps']))
0243 end
0244 
0245     
0246 end