TauARaster_quick2

 A close look at the VirgoA raster data
 Makes Az El plots

0001 function TauARaster_quick2(nchan,plot_labels,save_path,source_name,start_date,IntensityType,calmat)
0002 % A close look at the VirgoA raster data
0003 % Makes Az El plots
0004 [home,installeddir] = where_am_i();
0005 maindir= [home,'/',installeddir,'/',save_path,'/'];
0006 filename = ([maindir,start_date,'_',source_name,'_',IntensityType]);
0007 load([filename,'_dgood.mat'])
0008 %load([filename,'_gfit.mat']);
0009 load([filename,'_PA_IQUV.mat']);
0010    
0011 [antenna_no antenna_name] = identifyTelescope(dgood);
0012 telescope_constants
0013     
0014 %load('./test_fits/07-Nov-2012:07:22:02_TauA_dgood.mat')
0015 %load('./test_fits/07-Nov-2012:07:22:02_TauA_gfit.mat')
0016 %load('./test_fits/06-Nov-2012:09:34:32_moon_dgood.mat')
0017 %load('./test_fits/06-Nov-2012:09:34:32_moon_gfit.mat')
0018 % if(exist('calmat'))
0019 %     for i=1:length(dgood.antenna0.receiver.data(:,1))
0020 %         dgood.antenna0.receiver.data(i,1:3)=dgood.antenna0.receiver.data(i,1:3)*calmat;  % Mulitply I,Q,U by the calibration matrix
0021 %     end
0022 %
0023 % end
0024 % Fix the offsets in the data so that it can be binned coherently
0025 n=0
0026 cell_size=0.1
0027 clims{1} = [-1 1]*620;
0028 clims{2} = [-0.1 0.1]*620;
0029 clims{3} = [-0.1 0.1]*620;
0030 clims{4} = [-1 1]*620;
0031 
0032 [s e] = findStartStop(dgood.index.radio_point_scan.fast);
0033 %pols = ['I','Q','U','V','P','P'];
0034 % Set up some pointing arrays to save corrected Az and El Offsets
0035 azOffSaveP = zeros(size(dgood.antenna0.receiver.data(:,1)));
0036 elOffSaveP = zeros(size(dgood.antenna0.receiver.data(:,1)));
0037 dataP = zeros(size(dgood.antenna0.receiver.data(:,1)));
0038 AZ= zeros(size(dgood.antenna0.receiver.data(:,1)));
0039 EL= zeros(size(dgood.antenna0.receiver.data(:,1)));
0040 DEC= zeros(size(dgood.antenna0.receiver.data(:,1)));
0041 RA = zeros(size(dgood.antenna0.receiver.data(:,1)));
0042 points = zeros(size(dgood.antenna0.receiver.data(:,1)));
0043 %Identify the maximum of I,Q,U,V to use as a normalization for each raster
0044 %points = (dgood.index.radio_point_scan.fast==1);
0045 
0046 % contour_max(1) = max(dgood.antenna0.receiver.data(points,1));
0047 % contour_max(2) = max(dgood.antenna0.receiver.data(points,2));
0048 % contour_max(3) = max(dgood.antenna0.receiver.data(points,3));
0049 % contour_max(4) = max(dgood.antenna0.receiver.data(points,4));
0050 
0051 %Make sure that we subtract the offsets (as fit by a gaussian) from
0052 %all the az and el data, so it can be binned coherently
0053 dgood.azApp = interp1(dgood.antenna0.tracker.utc,dgood.antenna0.tracker.horiz_topo(:,1), dgood.antenna0.receiver.utc);
0054 dgood.azOffSave = dgood.azApp - dgood.antenna0.servo.apparent(:,1);
0055 dgood.azOffSave = -dgood.azOffSave;
0056 dgood.azOffSave = dgood.azOffSave.*cos(dgood.antenna0.servo.el*pi/180);% Scale azimuth by cos(elevation)???CHECK
0057 
0058 
0059 
0060 dgood.elApp = interp1(dgood.antenna0.tracker.utc,dgood.antenna0.tracker.horiz_topo(:,2), dgood.antenna0.receiver.utc);
0061 dgood.elOffSave = dgood.elApp - dgood.antenna0.servo.apparent(:,2);
0062 dgood.elOffSave = -dgood.elOffSave;
0063 angle = sqrt(dgood.elOffSave.^2 + dgood.azOffSave.^2);
0064 
0065 %for m= 1:length(s)
0066 for m=good_rasters;
0067 ind = zeros(size(dgood.index.radio_point_scan.fast));
0068 ind(s(m)+1000:e(m)) = 1; %% This bit is important - removes rogue repeat points stuck on source!!
0069 ind = logical(ind);
0070 dcut  = framecut(dgood, ind);
0071 % Make sure we have all the data for just the detrended raster bits.
0072 % Currently uses gfit for channel I1 % Need to check if big differences
0073 % between I channels - shouldnt be!!
0074 azOffSaveP(ind)=dcut.azOffSave - gfit{m,1}.x0*cell_size;
0075 elOffSaveP(ind)=dcut.elOffSave - gfit{m,1}.y0*cell_size;
0076 dataP(ind,1:nchan)=dcut.antenna0.receiver.data(:,1:nchan); 
0077 timeP(ind) = dcut.antenna0.receiver.utc;
0078 AZ(ind) = deg2rad(dcut.antenna0.servo.apparent(:,1)); % radians
0079 EL(ind) = deg2rad(dcut.antenna0.servo.apparent(:,2)); % radians
0080 
0081 % Going to be makeing Az-El plots at current epoch so make sure we have
0082 % that info in equa
0083 dcut=calcRADECcurrent(dcut); 
0084 DEC(ind) = (dcut.antenna0.servo.equa(:,2));% radians
0085 RA(ind)  = (dcut.antenna0.servo.equa(:,1));% radians
0086 %DEC(ind) = deg2rad(mean((dcut.antenna0.tracker.equat_geoc(:,2))));% radians
0087 points(ind)=1;
0088 %Plot each raster
0089 % figure(m)
0090 % for(chan=1:4)
0091 % subplot(2,2,chan)
0092 % xp=azOffSaveP(ind);
0093 % yp=elOffSaveP(ind);
0094 % data = dataP(ind,chan);
0095 % a = (isnan(xp) | isnan(yp));
0096 % [x y z h] = gridcbass2(xp,yp,data,0.05,2,1,1);
0097 % xx{chan}=x;
0098 % yy{chan}=y;
0099 % zz{chan}=z;
0100 % clear x y z;
0101 % zz{chan}(zz{chan}==0)=NaN;
0102 % bb{chan}=inpaint_nans(zz{chan},1);
0103 %
0104 % imagesc(bb{chan},clims{chan})
0105 %
0106 % colorbar
0107 % axis equal
0108 % end
0109 end
0110 
0111 points = logical(points);
0112 clear dcut
0113 clear dgood ind
0114 %%
0115 % Now make binned beam maps of I Q U P
0116 
0117 %%%
0118 %Bit to plot the timestream data at this point
0119 
0120 % figure(99)
0121 % for kk=1:3;
0122 % subplot(3,1,kk);
0123 % plot(dataP(:,kk));
0124 % mode_abs = mode(abs(dataP(:,kk)))
0125 % gtitle(['Time stream data read in before calibration -', source_name]);
0126 %  print('-f99','-dpng',[maindir,'Time_Stream_preFlag_',source_name,'_',start_date,'.png']);
0127 % end
0128 close all
0129 
0130 % Plot after flagging data .gt. abs(0.5)  unit in Q and U channels
0131 %dataP((abs(dataP(:,2))>0.1),:)=nan;
0132 %dataP((abs(dataP(:,3))>0.1),:)=nan;
0133 
0134 % figure(99)
0135 % for kk=1:3;
0136 % subplot(3,1,kk);
0137 % plot(dataP(:,kk));
0138 % gtitle(['Time stream data read in before calibration (flagged) -', source_name]);
0139 %  print('-f99','-dpng',[maindir,'Time_Stream_postFlag_',source_name,'_',start_date,'.png']);
0140 % end
0141 close all
0142 
0143 
0144 %%%NEED TO CHANGE THIS FOR DIFFERENT SIZE CAL MATRIX and use inverse from
0145 %%%what is used in descart
0146 if(exist('calmat'))
0147     calmat
0148     
0149     calmat1=calmat(:,[1,3,5]); % calibration matrix for I1,Q1,U1
0150     calmat2=calmat(:,[2,4,6]); % calibration matrix for I2,Q2,U2
0151     
0152     invcalmat1 = inv(calmat1); % Should calibrate raw observed data --> corrected data
0153     invcalmat2 = inv(calmat2);
0154     
0155     % Apply the inverse matrices to dataP (I1,Q1,U1) and dataP(I2,Q2,U2)
0156     % dataP is [I1,Q1,U1,Q2,U2,Q3,U3,I2,V]
0157     
0158     for i=1:length(dataP(:,1));
0159         dataP(i,[1,2,3])=dataP(i,[1,2,3])*invcalmat1;  % Mulitply I,Q,U by the calibration matrix
0160         dataP(i,[8,4,5])=dataP(i,[8,4,5])*invcalmat2;
0161     end
0162 end
0163 %%%
0164 
0165 % Also make the polarization, P1, P2, P from Q1U1, Q2U2, mena QU
0166 if(nchan==9)
0167    % re-creat V using calibrated data
0168    dataP(:,9) = (dataP(:,1)-dataP(:,8))/2 ; % calibrated V
0169    meanQ =  (dataP(:,2)+dataP(:,4))/2;
0170    meanU =  (dataP(:,3)+dataP(:,5))/2;
0171    meanI =  (dataP(:,1)+dataP(:,8))/2;
0172    dataP(:,10) = meanI;
0173    dataP(:,11) = meanQ;
0174    dataP(:,12) = meanU;
0175    dataP(:,13) = sqrt(dataP(:,2).^2 + dataP(:,3).^2); % P1 is pol from QU1
0176    dataP(:,14) = sqrt(dataP(:,4).^2 + dataP(:,5).^2); % P2 is pol from QU2
0177    dataP(:,15) = sqrt(meanU.^2 + meanQ.^2);% How do we want to form P (from mean Q1Q2,U1,U2??)
0178    labels=plot_labels;
0179    labels(10:18) = {'I','Q','U','P1 (timestream)','P2 (timestream)','P (timestream)','P1 (from Q1,U1 maps)','P2 (from Q2,U2 maps)','P (from Q,U maps)'},
0180    chanstoplot=[1,2,3,13,8,4,5,14,10,11,12,15]; % plot in order I1,Q1,U1,P1...I2,Q2,U2,etc
0181    plotnums =  [1,2,3,4 ,6,7,8,9 ,11,12,13,14];
0182     chanstofit=[1,8]; % fit I1 and I2 for gaussian beam
0183 else
0184     disp(['Need to fix this for CBASS-S'])
0185 end
0186 clear meanI meanQ meanU
0187 
0188 % % % Now select the data we want to use
0189 %
0190 % points = (dgood.index.radio_point_scan.fast==1);
0191 %
0192 % % % Also remove some rfi from this 'points' selected data
0193 %
0194 % points(1.69e6:2.545e6)=0;
0195 % points(2.32e5:2.352e5)=0;
0196 
0197 % Grab the data and bin it for each channel
0198 close all
0199 figure
0200 %labels={'I','Q','U','V','P (timestream)','P (from Q,U maps)'};
0201 num=0
0202 for (chan=chanstoplot);
0203     num=num+1;
0204     plotnum=plotnums(num);
0205  subplot(3,5,plotnum);
0206  xp=azOffSaveP(points);
0207  yp=elOffSaveP(points);
0208  data=dataP(points,chan);
0209  a = (isnan(xp) | isnan(yp));
0210  [x y z h] = gridcbass2(xp,yp,data,0.05,2,1,1);
0211 
0212  clear xp yp
0213 
0214 xx{chan}=x;
0215 yy{chan}=y;
0216 zz{chan}=z;
0217 clear x y z;
0218 
0219 zz{chan}(zz{chan}==0)=NaN;
0220 bb{chan}=inpaint_nans(zz{chan},1);
0221 %imagesc(bb{chan},clims{chan})
0222 imagesc(bb{chan});
0223 colorbar;
0224 xlabel(labels{chan});
0225 axis equal
0226 end
0227 plotnum=5;
0228 matnum=15;
0229 for chan=[2,4,11];
0230 matnum=matnum+1;
0231 bb{matnum}=(sqrt(bb{chan}.^2 +bb{chan+1}.^2)); % Make P map by sqrt(sum of Q and U maps squared)
0232 subplot(3,5,plotnum);
0233 imagesc(bb{matnum});
0234 colorbar;
0235 xlabel(labels{matnum});
0236 axis equal;
0237 plotnum=plotnum+5;
0238 if(exist('calmat'));
0239     gtitle(['Calibrated, NO PA rotation, binned maps -', source_name,'_',IntensityType])
0240     print('-f1','-dpng',[maindir,'CALMAT_noPArot_binned_map_',source_name,'_',IntensityType,'_',start_date,'.png'])
0241 else
0242     gtitle(['Un-calibrated, NO PA rotation, binned maps -', source_name,'_',IntensityType])
0243     print('-f1','-dpng',[maindir,'Uncal_noPArot_binned_map_',source_name,'_',IntensityType,'_',start_date,'.png'])
0244 end
0245     %close all
0246 end
0247 
0248 %%
0249 % Make an azimuthally binned beam map.
0250 
0251 % % Having got XI,YI,ZI then shift the data to be centred by useing datatick
0252 % to export centre:
0253 %figure
0254 pnum=0;
0255 
0256 for chan=chanstofit;
0257     pnum=pnum+1;
0258 %cursor_info
0259 XII=xx{1};
0260 %YII=yy{chan}-0.25;
0261 YII=yy{1};
0262 ZII=bb{chan};
0263 %surface(XII,YII,bb{chan})
0264 % Or fit a gaussian and grab the fitted centre
0265 [sf gof t] = gauss2D_fit_act(XII,YII,ZII);
0266 big_gfit{chan} = sf;
0267 x_off = big_gfit{chan}.x0;
0268 y_off = big_gfit{chan}.y0;
0269 XII = XII-x_off;
0270 YII = YII-y_off;
0271 
0272 angle = sqrt(XII.^2 + YII.^2);
0273 ang=reshape(angle,1,81*81); % Used to be 81*81
0274 zinew = reshape(bb{chan},1,81*81);% Used to be 81*81
0275 %extrap_max = interp1(ang,zinew,0,'pchip','extrap') %extrapolated beam max at angle=0
0276 extrap_max =1;
0277 figure;
0278 %plot(ang)
0279 %plot(ang,zinew)
0280 
0281 % Now we want to bin the data
0282 x = ang;
0283 y = zinew./extrap_max;
0284 
0285 topEdge = 5; % define limits in degrees
0286 botEdge = 0; % define limits
0287 %numBins = topEdge/0.05; % define number of bins
0288 binWidth = 0.05;
0289 numBins = topEdge/binWidth; % define number of bins
0290 
0291 binEdges = linspace(botEdge, topEdge, numBins+1);
0292 [h,whichBin] = histc(x, binEdges);
0293 
0294 for i = 1:numBins
0295     flagBinMembers = (whichBin == i);
0296     binMembers     = y(flagBinMembers);
0297     binMean(i)     = nanmean(binMembers);
0298     binStd(i)      = nanstd(binMembers);
0299     xbinMembers    = x(flagBinMembers);
0300     xbinMean(i)    = nanmean(xbinMembers);
0301 end
0302 
0303 % Save max of I beam for normalisation of other beams
0304 if(chan==1)
0305     Ibeam_max = max(binMean);
0306 end
0307 binCentres = binEdges +(binWidth)/2;
0308 %plot(binCentres(1:end-1),binMean./max(binMean))
0309 %plot(binEdges(1:end-1),binMean./max(binMean))
0310 %extrap_max = interp1(xbinMean(1:70),binMean(1:70),0,'pchip','extrap') %extrapolated beam max at angle=0
0311 %extrap_max = interp1(ang,zinew,0,'pchip','extrap') %extrapolated beam max at angle=0
0312 
0313 %plot(xbinMean,binMean)
0314 % Grab the simulated beam
0315 figure(19)
0316 subplot(1,2,pnum);
0317 n_cuts = 2;
0318 object='TauA';
0319 fit_lim = 0.8;
0320 [mean_power_ZP mean_power_P sim_beam sim_data sim_angles phi vals_sim FWHM_sim sim_fit_angles simfile] = read_sim_beam_freq(n_cuts, fit_lim,object);
0321 %errorbar(binEdges(1:end-1),binMean./max(binMean),binStd(1:end))
0322 
0323 plot(xbinMean,binMean./Ibeam_max);
0324 
0325 if(chan==1)
0326     errorbar(xbinMean,binMean./Ibeam_max,binStd(1:end)./Ibeam_max);
0327 hold all;
0328 plot(sim_beam(:,1),sim_beam(:,2)/max(sim_beam(:,2)));
0329 end
0330 grid on;
0331 xlim([0 4]);
0332 
0333 
0334 % %subplot(2,2,chan)
0335 % binCentres = binEdges +(binWidth)/2;
0336 %
0337 % plot(binEdges(1:end-1),binMean./max(binMean),'.')
0338 legend(labels{chan})
0339 end
0340 
0341 ylim([0 max(binMean./Ibeam_max)*1.2]);
0342 
0343    
0344 if(exist('calmat'))
0345 gtitle(['Calibrated,azimuthally binned beam profile (I) - ', source_name,'_',IntensityType]);
0346 
0347 print('-f19','-dpng',[maindir,'CALMAT_AzimuthalBeam_',source_name,'_',IntensityType,'_',start_date,'.png'])
0348 else
0349     gtitle(['Un-calibrated,azimuthally binned beam profile (I) - ', source_name,'_',IntensityType]);
0350     print('-f19','-dpng',[maindir,'AzimuthalBeam_',source_name,'_',IntensityType,'_',start_date,'.png'])
0351 end
0352 close all
0353 %% Plot things in terms of parallactic angle
0354 
0355 %dall = framecut(dgood,dgood.index.radio_point_scan.fast);
0356 %dall.flags.fast(1.345e5:1.989e5)=1; % in case you need more flags
0357 %dall.flags.fast(2.29e6:2.321e6)=1;
0358 %dall= framecut(dall,~dall.flags.fast(:,1));
0359 %plot(dall.antenna0.receiver.data(:,1),'r')
0360 
0361 I= dataP(points,10);
0362 Q= dataP(points,11);
0363 U= dataP(points,12);
0364 V= dataP(points,9);
0365 P= dataP(points,15); % P from timestream data
0366 time_cut = timeP(points);
0367 
0368 %AZ = deg2rad(dall.antenna0.servo.apparent(:,1)); % radians
0369 %EL = deg2rad(dall.antenna0.servo.apparent(:,2)); % radians
0370 
0371 %%%%%%%%%%%%
0372 % Calculate the parallactic angle
0373 %AZ = AZ(points);%radians
0374 %EL = EL(points);%radians
0375 DEC= DEC(points);%radians Current epoch
0376 RA = RA(points); % radians Current Epoch
0377 
0378 % Grab latitude and longitude of antenna
0379 
0380     
0381 LAT_D = antenna_latitude(antenna_no)
0382 LON_D = antenna_longitude(antenna_no)
0383 
0384 LAT = deg2rad(LAT_D);
0385 LON = deg2rad(LON_D);
0386 JD = mjd2jd(time_cut);
0387 LST=lst(JD,deg2rad(LON_D),'m'); % Uses East longitude in radians
0388                                     % Gives LST in fraction of days
0389     
0390 clear dataP timeP xx yy bb zz;
0391 clear mean_power_ZP mean_power_P sim_beam sim_data sim_angles phi vals_sim FWHM_sim sim_fit_angles simfile;
0392 clear JD;
0393 % Clear the working space and manage memeory
0394 cwd = pwd;
0395 cd(tempdir);
0396 pack
0397 cd(cwd)
0398 %
0399 whos
0400 [PA]=parallactic_angle(((RA)),((DEC)),(LST),(LAT)); % Everything must be in radians, output is in radians
0401 PA_D = rad2deg(PA); % gets definition in same way as real data. %% CHeck sign here - 13/8/2014 - ACT
0402 %%%%%%
0403 
0404 
0405 %%%%%%%%%%%%%%%
0406 % Rotate Q and U and calculate P
0407 Qc = Q.*cos(2*PA) - U.*sin(2*PA) ;
0408 Uc = U.*cos(2*PA) + Q.*sin(2*PA) ;
0409 Pc = sqrt(Qc.^2 + Uc.^2);
0410 
0411 % Now plot the data before and after parallactic sky rotation
0412 scrsz = get(0,'ScreenSize');
0413 
0414 
0415 close all
0416 figure(22)
0417 
0418 ax(1)=subplot(3,2,1);
0419 h(1)=plot(time_cut,I,'r');
0420 xlabel('I')
0421 if(~exist('calmat'))
0422 ylim([0 2])
0423 end
0424 
0425 ax(4)= subplot(3,2,2);
0426 h(4)=plot(time_cut,V,'k');
0427 xlabel('V ')
0428 
0429 
0430 ax(2)=subplot(3,2,3);
0431 h(2)=plot(time_cut,Q,'g');
0432 xlabel('Q - no PA correction')
0433 if(~exist('calmat'))
0434 ylim([-0.1 0.1])
0435 end
0436 
0437 ax(3)=subplot(3,2,4);
0438 h(3)=plot(time_cut,U,'b');
0439 xlabel('U - no PA correction')
0440 if(~exist('calmat'))
0441 ylim([-0.1 0.1])
0442 end
0443 
0444 ax(5)= subplot(3,2,5);
0445 h(5)=plot(time_cut,P,'k');
0446 xlabel('P  - no PA correction');
0447 if(~exist('calmat'))
0448 ylim([0 0.1])
0449 end
0450 
0451 ax(6)=subplot(3,2,6);
0452 h(6)=plot(time_cut,PA_D);
0453 xlabel('Parallactic angle.deg')
0454 
0455 linkaxes(ax,'x');
0456 
0457 if(exist('calmat'))
0458     gtitle('Calibrated - no PA rotation')
0459     print('-f22','-dpng',[maindir,'CALMAT_IQUV_1_',source_name,'_',IntensityType,'_',start_date,'.png'])
0460 else
0461     gtitle('No calibration - no PA rotation')
0462     print('-f22','-dpng',[maindir,'IQUV_1_',source_name,'_',IntensityType,'_',start_date,'.png'])
0463 end
0464 %close all
0465 
0466 figure(23)
0467 ax(1)=subplot(3,2,1);
0468 h(1)=plot(time_cut,I,'r');
0469 xlabel('I')
0470 if(~exist('calmat'))
0471 ylim([0 2])
0472 end
0473 
0474 ax(4)= subplot(3,2,2);
0475 h(4)=plot(time_cut,V,'k');
0476 xlabel('V ')
0477 
0478 ax(2)=subplot(3,2,2);
0479 h(2)=plot(time_cut,Qc,'g');
0480 xlabel('Q')
0481 if(~exist('calmat'))
0482 ylim([-0.1 0.1])
0483 end
0484 
0485 ax(3)=subplot(3,2,3);
0486 h(3)=plot(time_cut,Uc,'b');
0487 xlabel('U')
0488 if(~exist('calmat'))
0489 ylim([-0.1 0.1])
0490 end
0491 
0492 ax(4)= subplot(3,2,4);
0493 h(4)=plot(time_cut,Pc,'k');
0494 xlabel('P')
0495 if(~exist('calmat'))
0496 ylim([0 0.1])
0497 end
0498 
0499 ax(5)=subplot(3,2,5);
0500 h(5)=plot(time_cut,PA_D);
0501 xlabel('Parallactic angle, deg')
0502 linkaxes(ax,'x');
0503 
0504 if(exist('calmat'))
0505     gtitle('Calibrated - PA rotated')
0506 print('-f23','-dpng',[maindir,'CALMAT_IQUV_2_',source_name,'_',IntensityType,'_',start_date,'.png'])
0507 else
0508     gtitle('No calibration - PA rotated')
0509     print('-f23','-dpng',[maindir,'IQUV_2_',source_name,'_',IntensityType,'_',start_date,'.png'])
0510 end
0511 
0512 %%
0513 %Make gridded maps of Qc,Uc,Pc
0514 close all
0515 
0516 datac = [I Qc Uc V Pc]; % Timestream data
0517 labels_small={'I','Q','U','V','P timestream','P maps'};
0518 for (chan=1:5)
0519 subplot(2,3,chan);
0520  xp=azOffSaveP(points);
0521  yp=elOffSaveP(points);
0522  data = datac(:,chan);
0523  clear datac;
0524  a = (isnan(xp) | isnan(yp));
0525  [x y z h] = gridcbass2(xp,yp,data,0.1,2,1,1);
0526 %xp=azOffSave(points);
0527 %yp=elOffSave(points);
0528 %data = Qc;
0529 %a = (isnan(xp) | isnan(yp));
0530 %[x y z h] = gridcbass2(xp,yp,Qc,0.05,2,1,1);
0531 
0532 
0533 
0534 xx{chan}=x;
0535 yy{chan}=y;
0536 zz{chan}=z;
0537 clear x y z;
0538 
0539 zz{chan}(zz{chan}==0)=NaN;
0540 bb{chan}=inpaint_nans(zz{chan},1);
0541 %imagesc(bb{chan},clims{chan})
0542 imagesc(bb{chan});
0543 xlabel(labels_small{chan});
0544 colorbar;
0545 axis equal
0546 
0547 end
0548 
0549 % Also map the P map made by combining the Q and U maps
0550 bb{6}=(sqrt(bb{2}.^2 +bb{3}.^2)) ;% Make P map by sqrt(sum of Q and U maps squared)
0551 chan=6;
0552 subplot(2,3,chan);
0553 imagesc(bb{chan});
0554 colorbar;
0555 xlabel(labels_small{chan});
0556 axis equal;
0557 
0558 if(exist('calmat'))
0559     gtitle(['Calibrated and PA rotated binned maps - ', source_name,'_',IntensityType]);
0560     print('-f1','-dpng',[maindir,'CALMAT_Binned_PArotated_map_',source_name,'_',IntensityType,'_',start_date,'.png'])
0561 else
0562     gtitle(['Un-calibrated but PA rotated binned maps - ', source_name,'_',IntensityType]);
0563    print('-f1','-dpng',[maindir,'Binned_PArotated_map_',source_name,'_',IntensityType,'_',start_date,'.png'])   
0564 end
0565 close all
0566 
0567 %Polarization angle plot
0568 ZI_Q = bb{2};
0569 ZI_U = bb{3};
0570 ZI_P = bb{5};
0571 
0572 figure(100)
0573  Chi = 0.5*atan(ZI_U./ZI_Q);
0574  [u,v] = pol2cart(Chi,ZI_P);
0575 % quiver(XI,YI,u,v);
0576 %
0577 % figure
0578  contour(ZI_P);
0579  hold on
0580  quiver(u,v);
0581  %title('P(from Q and U maps) and Polarization angle')
0582  
0583  if(exist('calmat'))
0584     gtitle(['Calibrated and PA rotated Polarization map - ', source_name,'_',IntensityType]);
0585     print('-f100','-dpng',[maindir,'CALMAT_Binned_PArotated_Polmap_',source_name,'_',IntensityType,'_',start_date,'.png'])
0586 else
0587     gtitle(['Un-calibrated but PA rotated Polarization map - ', source_name,'_',IntensityType]);
0588    print('-f100','-dpng',[maindir,'Binned_PArotated_Polmap_',source_name,'_',IntensityType,'_',start_date,'.png'])   
0589 end
0590     
0591  
0592 
0593 
0594 %Make final plot P
0595 % P = sqrt(dgood.antenna0.receiver.data(:,3).^2 + dgood.antenna0.receiver.data(:,2).^2);
0596 % dataP = P(points);
0597 % subplot(2,2,4)
0598 %  xp=azOffSaveP(points);
0599 %  yp=elOffSaveP(points);
0600 %  data = dataP;
0601 %  a = (isnan(xp) | isnan(yp));
0602 %  [x y z h] = gridcbass2(xp,yp,data,0.05,2,1,1);
0603 %
0604 % xx{chan}=x;
0605 % yy{chan}=y;
0606 % zz{chan}=z;
0607 % clear x y z;
0608 %
0609 % zz{chan}(zz{chan}==0)=NaN;
0610 % bb{chan}=inpaint_nans(zz{chan},1);
0611 % %imagesc(bb{chan},clims{chan})
0612 % imagesc(bb{chan})
0613 % colorbar
0614 % axis equal
0615 
0616 
0617 
0618 
0619 % %%
0620 % %Should now produce a calibrated parallactic angle rotated data set
0621 %   ti=-3:cell_size:3;
0622 %  [XI,YI]=meshgrid(ti,ti);
0623 %  a = (isnan(azOffSave) | isnan(elOffSave));
0624 %  ZI_P = griddata(azOffSave(~a),elOffSave(~a),Pc(~a),XI,YI);
0625 %  ZI_Q = griddata(azOffSave(~a),elOffSave(~a),Qc(~a),XI,YI);
0626 %  ZI_U = griddata(azOffSave(~a),elOffSave(~a),Uc(~a),XI,YI);
0627 %  ZI_I = griddata(azOffSave(~a),elOffSave(~a),dall.antenna0.receiver.data(~a,1),XI,YI);
0628 % figure
0629 % subplot(2,2,1)
0630 % contour(XI,YI,ZI_I)
0631 % title('Intensity')
0632 % subplot(2,2,2)
0633 % contour(XI,YI,ZI_P)
0634 % title('Polarization')
0635 % subplot(2,2,3)
0636 % contour(XI,YI,ZI_Q)
0637 % title('Q')
0638 % subplot(2,2,4)
0639 % contour(XI,YI,ZI_U)
0640 % title('U')