source_raster_remLeak_v2

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  function source_raster(start, stop)

    offset should be in sample numbers

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function dcal = source_raster(start, stop, leak_matrix, savefile, plotdir, dcal)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 %  function source_raster(start, stop)
0006 %
0007 %    offset should be in sample numbers
0008 %
0009 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0010 
0011 if(nargin==2)
0012   leak_matrix = [554.8911 20.3628 -22.6129; -10.3716  599.9700 547.3359; -157.4338 451.8572 -616.9229];    
0013   savefile = 'inprogess_testN.mat';
0014   plotdir  = '/home/sjcm/cbass_analysis/';
0015   offset = 0;
0016 end
0017 
0018 
0019 if(nargin<6)
0020   % first we read in the data
0021   if(iscell(start))
0022     for m=1:length(start)
0023       if(m==1)
0024     d = read_arc(start{m}, stop{m});
0025       else
0026     d1 = read_arc(start{m}, stop{m});
0027     d  = catstruct(1, [d d1]);
0028       end
0029     end
0030   else
0031     d = read_arc(start, stop);
0032   end
0033   
0034   % next we have to change the feature in the schedule
0035   ind = bitand(d.array.frame.features, 2^7)>0;
0036   if(length(find(ind))~=0)
0037     d.array.frame.features(ind) = d.array.frame.features(ind)-2^7+2^9;
0038   end
0039   ind = bitand(d.array.frame.features, 2^15)>0;
0040   if(length(find(ind))~=0)
0041     d.array.frame.features(ind) = d.array.frame.features(ind)-2^15+2^9;
0042   end
0043 
0044   % next we do some pipeline
0045   dcal = reduceData(d, 'reduc/rasterScript.red');
0046   
0047   % we will want the RA/DEC for all data
0048   display('Calculating RA/DEC');
0049   long=-118.2822;
0050   lat=37.2339;
0051   az = dcal.antenna0.servo.apparent(:,1);
0052   el = dcal.antenna0.servo.apparent(:,2);
0053   jd=mjd2jd(dcal.antenna0.receiver.utc);
0054   [equa] = horiz_coo([pi/180*(az) pi/180*(el)],jd,[pi/180*(long) ...
0055     pi/180*(lat)],'e');
0056   dcal.antenna0.servo.equa = equa;
0057 end
0058 
0059 % next, we want to split up each individual "raster scan"
0060 ind = bitand(dcal.array.frame.features, 2^9) > 0;
0061 [si ei] = findStartStop(ind);
0062 
0063 display(sprintf('There are %d rasters in this data set', length(si)));
0064 
0065 % define grids and things that take up the most memory
0066 xgrid = -5.5:0.1:5.5;
0067 ygrid = xgrid;
0068 [X, Y] = meshgrid(xgrid,ygrid);
0069 
0070 mapI     = nan(length(si), length(xgrid)-1, length(xgrid)-1);
0071 mapIvar  = mapI; 
0072 mapV     = mapI;  
0073 mapVrem  = mapI;  
0074 mapVvar  = mapI; 
0075 mapVremvar=mapI;
0076 mapQ     = mapI;  
0077 mapU     = mapI;  
0078 mapQrem  = mapI;  
0079 mapUrem  = mapI;  
0080 mapQvar  = mapI; 
0081 mapUvar  = mapI; 
0082 mapQremvar = mapI; 
0083 mapUremvar = mapI; 
0084 mapQp1   = mapI; 
0085 mapUp1   = mapI; 
0086 mapQp1var= mapI; 
0087 mapUp1var= mapI; 
0088 mapQp2   = mapI; 
0089 mapUp2   = mapI; 
0090 mapQp2var= mapI; 
0091 mapUp2var= mapI; 
0092 
0093 
0094 for m=1:length(si)
0095   ind = zeros(size(dcal.array.frame.features));
0096   ind(si(m):ei(m)) = 1;
0097   ind = logical(ind);
0098   dcut = framecut(dcal, ind);
0099 
0100   % next let's cut out when we're between scans
0101   da = abs(deriv(dcut.antenna0.servo.az));
0102   ind = da>0.0094 & da<0.0108;
0103   dcut = framecut(dcut, ind);
0104 
0105   %if the telescope is near a wrap, then it'll be just slewing around, and ...
0106       %  not doing us any good.
0107   perGood = length(find(ind))./length(da)*100;
0108 
0109   allGoodPer(m) = perGood;
0110   
0111   if(perGood > 15)
0112   % do the rest.
0113 
0114   % remove the baseline in intensity
0115   I1 = dcut.antenna0.receiver.data(:,1);
0116   I2 = dcut.antenna0.receiver.data(:,8);
0117   V  = dcut.antenna0.receiver.data(:,1) - dcut.antenna0.receiver.data(:,8);
0118   Q  = dcut.antenna0.receiver.data(:,6);
0119   U  = dcut.antenna0.receiver.data(:,7);
0120   t = (dcut.antenna0.receiver.utc - dcut.antenna0.receiver.utc(1))*24;
0121   
0122   % given that there might be pointing offsets, we want to:
0123   % 1.  make a map to find the center
0124   % 2.  calculate the sky offsets based on that center point
0125   % 3.  save each map so we can then display them.
0126   ra  = dcut.antenna0.servo.equa(:,1)*180/pi;
0127   dec = dcut.antenna0.servo.equa(:,2)*180/pi;
0128   % convert to sky/xy
0129   rafc = mean(ra);
0130   decfc = mean(dec);
0131   aa   = I1+ I2;
0132   mm   = find(aa == max(aa));
0133   rafc = dcut.antenna0.servo.equa(mm,1)*180/pi;
0134   decfc = dcut.antenna0.servo.equa(mm,2)*180/pi;
0135   [x y] = radec_to_fieldoff(ra, dec, rafc, decfc);
0136   
0137   % also need to calculate the parallactic angle.
0138   thisra = rafc*pi/180;
0139   thisdec = decfc*pi/180;
0140   lat     = mean(dcut.antenna0.tracker.siteActual(:,2))*pi/180;
0141   lst     = dcut.antenna0.tracker.lst(ceil(mm/100))*15*pi/180;
0142   ha  = (lst - thisra);
0143   cosel_cospa = sin(lat)*cos(thisdec) - sin(thisdec)*cos(lat)*cos(ha);
0144   cosel_sinpa = sin(ha)*cos(lat);
0145   pa  = atan2(cosel_sinpa, cosel_cospa);
0146   pa = pa*180/pi;
0147   allpa(m) = pa;
0148   
0149   % remove the baseline from points which are not within 1 degree of
0150   % source
0151   angdist = pyth(x,y);
0152   indfit  = angdist > 1;
0153   [slope, b] = linfit(t(indfit), I1(indfit));
0154   I1rem = I1 - (slope*t+b);
0155   [slope, b] = linfit(t(indfit), I2(indfit));
0156   I2rem = I2 - (slope*t+b);
0157   Vrem  = I1rem - I2rem;
0158   [slope, b] = linfit(t(indfit), Q(indfit));
0159   Qrem  = Q - (slope*t+b);
0160   [slope, b] = linfit(t(indfit), U(indfit));
0161   Urem  = U - (slope*t+b);  
0162   Itot  = I1rem + I2rem;
0163   
0164   % correct for leakage
0165   datamatrix = [Itot Qrem Urem];
0166   dataLeak   = datamatrix*leak_matrix;
0167   
0168   Irem = dataLeak(:,1);
0169   Qrem = dataLeak(:,2);
0170   Urem = dataLeak(:,3);
0171     
0172   % next we have to rotate Q and U by the parallactic angle (2X, actually)
0173   % % MOD TO -2X FOR TEST.
0174   Qprime1 = cosd(-2*pa)*Qrem - sind(-2*pa)*Urem;
0175   Uprime1 = sind(-2*pa)*Qrem + cosd(-2*pa)*Urem;
0176 
0177   % now that things are in sky offsets, we want to fit a 2D gaussian to find
0178   % the peak.
0179   % now that things are in sky offsets, we want to fit a 2D gaussian to find
0180   % the peak.
0181   sigmaFit = 0.73/(2*sqrt(2*log(2)));
0182   LB = [0 -0.5  -0.5 0.1 0.1];
0183   UB = [inf 0.5  0.5 2 2];
0184   dataCell{1} = x; 
0185   dataCell{2} = y;
0186   dataCell{3} = [1 1 1 1 1];
0187   
0188   dataCell{4} = [max(I1rem) 0 0 sigmaFit sigmaFit]; % initGuess
0189   [fitVals(m,1,:), resnorm, resid1, exitflag, output, lambda, J1] = ...
0190       lsqcurvefit('gauss2D_fit_nograd', [max(I1rem) 0 0 sigmaFit sigmaFit], ...
0191       dataCell, I1rem, LB, UB);
0192   
0193   dataCell{4} = [max(I2rem) 0 0 sigmaFit sigmaFit]; % initGuess
0194   [fitVals(m,2,:), resnorm, resid1, exitflag, output, lambda, J1] = ...
0195       lsqcurvefit('gauss2D_fit_nograd', [max(I2rem) 0 0 sigmaFit sigmaFit], ...
0196       dataCell, I2rem, LB, UB);
0197   
0198   dataCell{4} = [max(Itot) 0 0 sigmaFit sigmaFit]; % initGuess
0199   [fitVals(m,3,:), resnorm, resid1, exitflag, output, lambda, J1] = ...
0200       lsqcurvefit('gauss2D_fit_nograd', [max(Itot) 0 0 sigmaFit sigmaFit], ...
0201       dataCell, Itot, LB, UB);
0202   
0203   if(abs(mean(fitVals(m,:,2))) < 0.3)
0204     xoff(m) = mean(fitVals(m,:,2));
0205   else
0206     xoff(m) = xoff(m-1);
0207   end
0208   if(abs(mean(fitVals(m,:,3))) < 0.3)
0209     yoff(m) = mean(fitVals(m,:,3));
0210   else
0211     yoff(m) = yoff(m-1);
0212   end
0213     
0214   
0215   % Next we shift everything so that we're centered on the source.
0216   x = x-xoff(m);
0217   y = y-yoff(m);
0218 
0219 
0220   % now we put it on a grid and grid it.
0221   mapI1  = nan(length(xgrid)-1, length(xgrid)-1);
0222   mapV1  = nan(length(xgrid)-1, length(xgrid)-1);
0223   mapV1rem= nan(length(xgrid)-1, length(xgrid)-1);
0224   mapQ1  = nan(length(xgrid)-1, length(xgrid)-1);
0225   mapU1  = nan(length(xgrid)-1, length(xgrid)-1);
0226   mapQ1p1 = nan(length(xgrid)-1, length(xgrid)-1);
0227   mapV1p1 = nan(length(xgrid)-1, length(xgrid)-1);
0228   mapQ1rem= nan(length(xgrid)-1, length(xgrid)-1);
0229   mapU1rem= nan(length(xgrid)-1, length(xgrid)-1);
0230   
0231   
0232   mapI1var  = nan(length(xgrid)-1, length(xgrid)-1);
0233   mapV1var  = nan(length(xgrid)-1, length(xgrid)-1);
0234   mapV1remvar  = nan(length(xgrid)-1, length(xgrid)-1);
0235   mapQ1var  = nan(length(xgrid)-1, length(xgrid)-1);
0236   mapU1var  = nan(length(xgrid)-1, length(xgrid)-1);
0237   mapQ1p1var = nan(length(xgrid)-1, length(xgrid)-1);
0238   mapV1p1var = nan(length(xgrid)-1, length(xgrid)-1);
0239   mapQ1remvar= nan(length(xgrid)-1, length(xgrid)-1);
0240   mapU1remvar= nan(length(xgrid)-1, length(xgrid)-1);
0241   
0242   x4label = xgrid+0.1/2;
0243   x4label(length(x4label)) = [];
0244   for mm=2:length(xgrid)
0245     if(mm==2)
0246       fx = x>=xgrid(1) & x<xgrid(2);
0247     else
0248       fx = x>=xgrid(mm-1) & x<xgrid(mm);
0249     end
0250     for mmm=2:length(ygrid)
0251       if(mmm==2)
0252     fy = y>=ygrid(1) & y<ygrid(2);
0253       else
0254     fy = y>=ygrid(mmm-1) & y<ygrid(mmm);
0255       end
0256       f = find(fx&fy);
0257       mapI1(mm-1,mmm-1)   = nanmean(Irem(f));
0258       mapV1(mm-1,mmm-1)   = nanmean(V(f));    
0259       mapV1rem(mm-1,mmm-1)   = nanmean(Vrem(f));    
0260       mapQ1(mm-1,mmm-1)   = nanmean(Q(f));
0261       mapU1(mm-1,mmm-1)   = nanmean(U(f));
0262       mapQ1p1(mm-1,mmm-1)  = nanmean(Qprime1(f));
0263       mapU1p1(mm-1,mmm-1)  = nanmean(Uprime1(f));
0264 %      mapQ1p2(mm-1,mmm-1)  = nanmean(Qprime2(f));
0265 %      mapU1p2(mm-1,mmm-1)  = nanmean(Uprime2(f));
0266       mapQ1rem(mm-1,mmm-1)  = nanmean(Qrem(f));
0267       mapU1rem(mm-1,mmm-1)  = nanmean(Urem(f));
0268       mapI1var(mm-1,mmm-1)   = nanvar(Irem(f));
0269       mapV1var(mm-1,mmm-1)   = nanvar(V(f));
0270       mapV1remvar(mm-1,mmm-1)   = nanvar(Vrem(f));
0271       mapQ1var(mm-1,mmm-1)   = nanvar(Q(f));
0272       mapU1var(mm-1,mmm-1)   = nanvar(U(f));      
0273       mapQ1p1var(mm-1,mmm-1)  = nanvar(Qprime1(f));
0274       mapU1p1var(mm-1,mmm-1)  = nanvar(Uprime1(f));
0275       mapQ1remvar(mm-1,mmm-1)  = nanvar(Qrem(f));
0276       mapU1remvar(mm-1,mmm-1)  = nanvar(Urem(f));
0277       %     mapQ1p2var(mm-1,mmm-1)  = nanvar(Qprime2(f));
0278  %     mapU1p2var(mm-1,mmm-1)  = nanvar(Uprime2(f));
0279     end
0280   end
0281   
0282   mapI(m,:,:) = mapI1; 
0283   mapIvar(m,:,:) = mapI1var;
0284   mapV(m,:,:) = mapV1;
0285   mapVvar(m,:,:) = mapV1var;
0286   mapVrem(m,:,:) = mapV1rem;
0287   mapVremvar(m,:,:) = mapV1remvar;
0288   mapQ(m,:,:) = mapQ1;
0289   mapU(m,:,:) = mapU1;
0290   mapQvar(m,:,:) = mapQ1var;
0291   mapUvar(m,:,:) = mapU1var;
0292   mapQrem(m,:,:)   = mapQ1rem;
0293   mapUrem(m,:,:)   = mapU1rem;
0294   mapQp1(m,:,:) = mapQ1p1;
0295   mapUp1(m,:,:) = mapU1p1;
0296   mapQp1var(m,:,:) = mapQ1p1var;
0297   mapUp1var(m,:,:) = mapU1p1var;
0298 %  mapQp2(m,:,:) = mapQ1p2;
0299 %  mapUp2(m,:,:) = mapU1p2;
0300 %  mapQp2var(m,:,:) = mapQ1p2var;
0301 %  mapUp2var(m,:,:) = mapU1p2var;
0302   mapQremvar(m,:,:)   = mapQ1remvar;  
0303   mapUremvar(m,:,:)   = mapU1remvar;  
0304   
0305   
0306   figure('visible', 'off')
0307   imagesc(x4label, x4label, squeeze(mapI(m,:,:)));
0308   axis image; colorbar; xlabel('x-offset (degrees)');
0309   ylabel('y-offset (degrees)');
0310   txt = sprintf('I, pa=%3.2f, Scan #%d of %d', pa, m, length(si) );
0311   title(txt);
0312   filename = sprintf('%s/I_%02d.png', plotdir, m);
0313   eval(sprintf('print -dpng %s', filename));
0314 
0315   
0316   imagesc(x4label, x4label, squeeze(mapQ(m,:,:)));
0317   axis image; colorbar; xlabel('x-offset (degrees)');
0318   ylabel('y-offset (degrees)');
0319   txt = sprintf('Q, pa=%3.2f, Scan #%d of %d', pa, m, length(si) );
0320   title(txt);
0321   filename = sprintf('%s/Q_%02d.png', plotdir, m);
0322   eval(sprintf('print -dpng %s', filename));
0323   
0324   imagesc(x4label, x4label, squeeze(mapU(m,:,:)));
0325   axis image; colorbar; xlabel('x-offset (degrees)');
0326   ylabel('y-offset (degrees)');
0327   txt = sprintf('U, pa=%3.2f, Scan #%d of %d', pa, m, length(si) );
0328   title(txt);
0329   filename = sprintf('%s/U_%02d.png', plotdir, m);
0330   eval(sprintf('print -dpng %s', filename));
0331   
0332   imagesc(x4label, x4label, squeeze(mapQp1(m,:,:)));
0333   axis image; colorbar; xlabel('x-offset (degrees)');
0334   ylabel('y-offset (degrees)');
0335   txt = sprintf('Q rotated by pa=%3.2f, Scan #%d of %d', pa, m, length(si) );
0336   title(txt);
0337   filename = sprintf('%s/Qp1_%02d.png', plotdir, m);
0338   eval(sprintf('print -dpng %s', filename));
0339   
0340   imagesc(x4label, x4label, squeeze(mapUp1(m,:,:)));
0341   axis image; colorbar; xlabel('x-offset (degrees)');
0342   ylabel('y-offset (degrees)');
0343   txt = sprintf('U rotated by pa=%3.2f, Scan #%d of %d', pa, m, length(si) );
0344   title(txt);
0345   filename = sprintf('%s/Up1_%02d.png', plotdir, m);
0346   eval(sprintf('print -dpng %s', filename));
0347 
0348 %  figure(2)
0349 %  imagesc(x4label, x4label, squeeze(mapQp2(m,:,:)));
0350 %  axis image; colorbar; xlabel('x-offset (degrees)');
0351 %  ylabel('y-offset (degrees)');
0352 %  txt = sprintf('Q rotated by pa=-%3.2f, Scan #%d of %d', pa, m, length(si) );
0353 %  title(txt);
0354 %  filename = sprintf('%s/Qp2_%02d.png', plotdir, m);
0355 %  eval(sprintf('print -dpng %s', filename));
0356 %
0357 %  figure(3)
0358 %  imagesc(x4label, x4label, squeeze(mapUp2(m,:,:)));
0359 %  axis image; colorbar; xlabel('x-offset (degrees)');
0360 %  ylabel('y-offset (degrees)');
0361 %  txt = sprintf('U rotated by pa=-%3.2f, Scan #%d of %d', pa, m, length(si) );
0362 %  title(txt);
0363 %  filename = sprintf('%s/Up2_%02d.png', plotdir, m);
0364 %  eval(sprintf('print -dpng %s', filename));
0365   
0366   imagesc(x4label, x4label, squeeze(mapQrem(m,:,:)));
0367   axis image; colorbar; xlabel('x-offset (degrees)');
0368   ylabel('y-offset (degrees)');
0369   txt = sprintf('Q baseline removed, pa=-%3.2f, Scan #%d of %d', pa, m, length(si) );
0370   title(txt);
0371   filename = sprintf('%s/Qrem_%02d.png', plotdir, m);
0372   eval(sprintf('print -dpng %s', filename));
0373   
0374   imagesc(x4label, x4label, squeeze(mapUrem(m,:,:)));
0375   axis image; colorbar; xlabel('x-offset (degrees)');
0376   ylabel('y-offset (degrees)');
0377   txt = sprintf('U baseline removed, pa=-%3.2f, Scan #%d of %d', pa, m, length(si) );
0378   title(txt);
0379   filename = sprintf('%s/Urem_%02d.png', plotdir, m);
0380   eval(sprintf('print -dpng %s', filename));
0381 
0382   imagesc(x4label, x4label, squeeze(mapV(m,:,:)));
0383   axis image; colorbar; xlabel('x-offset (degrees)');
0384   ylabel('y-offset (degrees)');
0385   txt = sprintf('V, pa=%3.2f, Scan #%d of %d', pa, m, length(si) );
0386   title(txt);
0387   filename = sprintf('%s/V_%02d.png', plotdir, m);
0388   eval(sprintf('print -dpng %s', filename));
0389   
0390   imagesc(x4label, x4label, squeeze(mapVrem(m,:,:)));
0391   axis image; colorbar; xlabel('x-offset (degrees)');
0392   ylabel('y-offset (degrees)');
0393   txt = sprintf('V, baseline removed, pa=%3.2f, Scan #%d of %d', pa, m, length(si) );
0394   title(txt);
0395   filename = sprintf('%s/Vrem_%02d.png', plotdir, m);
0396   eval(sprintf('print -dpng %s', filename));
0397 
0398 end  % per good clause
0399 end 
0400 
0401 eval(sprintf('save %s/%s', plotdir, savefile));
0402 
0403 % now that we have all of those images generated, let's make the composite
0404 % and investigate the leakage, shall we?
0405 
0406 % make the composite map.  Need to do it as a weighted mean.
0407 % for the ones which have a variance of 0 (only one entry), we should set it
0408 % to a typical value of the variance in that data.
0409 
0410 for m=1:length(si)
0411   aa = mapIvar(m,:,:);
0412   medval = nanmedian(aa(:));
0413   maxval = nanmax(aa(:));
0414   f = find(mapIvar(m,:,:)==0);
0415   mapIvar(m,f) = mean([maxval, medval]);
0416   
0417   aa = mapQvar(m,:,:);
0418   maxval = nanmax(aa(:));
0419   f = find(aa == 0);
0420   mapQvar(m,f) = maxval;
0421 
0422   aa = mapQp1var(m,:,:);
0423   maxval = nanmax(aa(:));
0424   f = find(aa == 0);
0425   mapQp1var(m,f) = maxval;
0426   
0427   aa = mapQremvar(m,:,:);
0428   maxval = nanmax(aa(:));
0429   f = find(aa == 0);
0430   mapQremvar(m,f) = maxval;
0431 
0432 %  aa = mapQp2var(m,:,:);
0433 %  maxval = nanmax(aa(:));
0434 %  f = find(aa == 0);
0435 %  mapQp2var(m,f) = maxval;
0436   
0437   aa = mapUvar(m,:,:);
0438   maxval = nanmax(aa(:));
0439   f = find(aa == 0);
0440   mapUvar(m,f) = maxval;
0441   
0442   aa = mapUremvar(m,:,:);
0443   maxval = nanmax(aa(:));
0444   f = find(aa == 0);
0445   mapUremvar(m,f) = maxval;
0446 
0447   aa = mapUp1var(m,:,:);
0448   maxval = nanmax(aa(:));
0449   f = find(aa == 0);
0450   mapUp1var(m,f) = maxval;
0451 
0452 %  aa = mapUp2var(m,:,:);
0453 %  maxval = nanmax(aa(:));
0454 %  f = find(aa == 0);
0455 %  mapUp2var(m,f) = maxval;
0456   
0457   aa = mapVvar(m,:,:);
0458   maxval = nanmax(aa(:));
0459   f = find(aa == 0);
0460   mapVvar(m,f) = maxval;
0461 
0462   aa = mapVremvar(m,:,:);
0463   maxval = nanmax(aa(:));
0464   f = find(aa == 0);
0465   mapVremvar(m,f) = maxval;
0466 end
0467 
0468 % Next we make our maps
0469 [I_weight I_weight_var] = vwstat(mapI, mapIvar,1);
0470 [Q_weight Q_weight_var]     = vwstat(mapQ, mapQvar, 1);
0471 [Qrem_weight Qrem_weight_var]     = vwstat(mapQrem, mapQremvar, 1);
0472 [Qp1_weight Qp1_weight_var] = vwstat(mapQp1, mapQp1var, 1);
0473 %[Qp2_weight Qp2_weight_var] = vwstat(mapQp2, mapQp2var, 1);
0474 [U_weight U_weight_var]     = vwstat(mapU, mapUvar, 1);
0475 [Urem_weight Urem_weight_var]     = vwstat(mapUrem, mapUremvar, 1);
0476 [Up1_weight Up1_weight_var] = vwstat(mapUp1, mapUp1var, 1);
0477 [Up2_weight Up2_weight_var] = vwstat(mapUp2, mapUp2var, 1);
0478 [V_weight V_weight_var]     = vwstat(mapV, mapVvar, 1);
0479 [Vrem_weight Vrem_weight_var]     = vwstat(mapVrem, mapVremvar, 1);
0480 
0481 I_mean = nanmean(mapI,1);
0482 I_var  = nanvar(mapI,1);
0483 Q_mean = nanmean(mapQ,1);
0484 Qrem_mean = nanmean(mapQrem, 1);
0485 Qp1_mean  = nanmean(mapQp1,1);
0486 Qp1_var  = nanvar(mapQp1,1);
0487 U_mean = nanmean(mapU,1);
0488 Urem_mean = nanmean(mapUrem, 1);
0489 Up1_mean  = nanmean(mapUp1,1);
0490 Up1_var  = nanvar(mapUp1,1);
0491 V_mean = nanmean(mapV,1);
0492 Vrem_mean = nanmean(mapVrem,1);
0493 
0494 
0495 Ptot = sqrt(Up1_mean.^2 + Qp1_mean.^2);
0496 Ptot_weight = sqrt(Up1_weight.^2 + Qp1_weight.^2);
0497 Ptot_rms    = sqrt(Qp1_var + Up1_var);
0498 Ptot_weight_rms    = sqrt(Qp1_weight_var + Up1_weight_var);
0499 
0500 
0501 Psig = squeeze(Ptot./Ptot_rms);
0502 Psig_weight = squeeze(Ptot_weight./Ptot_weight_rms);
0503 
0504 chi1 = sum(sum(Psig(40:70,40:70).^2));
0505 chi1_weight = sum(sum(Psig_weight(40:70,40:70).^2));
0506 
0507 eval(sprintf('save %s/%s', plotdir, savefile));
0508 
0509 
0510 
0511 
0512 % FIGURES!!!
0513 
0514 imagesc(x4label, x4label, squeeze(I_weight));
0515 axis image; colorbar; xlabel('x-offset (degrees)');
0516 ylabel('y-offset (degrees)');
0517 title('Full I');
0518 filename = sprintf('%s/I_weighted.png', plotdir);
0519 eval(sprintf('print -dpng %s', filename));
0520 
0521 imagesc(x4label, x4label, squeeze(Q_weight));
0522 axis image; colorbar; xlabel('x-offset (degrees)');
0523 ylabel('y-offset (degrees)');
0524 title('Full Q');
0525 filename = sprintf('%s/Q_weighted.png', plotdir);
0526 eval(sprintf('print -dpng %s', filename));
0527 
0528 imagesc(x4label, x4label, squeeze(Qrem_weight));
0529 axis image; colorbar; xlabel('x-offset (degrees)');
0530 ylabel('y-offset (degrees)');
0531 title('Full Q - baseline removed');
0532 filename = sprintf('%s/Qrem_weighted.png', plotdir);
0533 eval(sprintf('print -dpng %s', filename));
0534 
0535 imagesc(x4label, x4label, squeeze(Qp1_weight));
0536 axis image; colorbar; xlabel('x-offset (degrees)');
0537 ylabel('y-offset (degrees)');
0538 title('Full Q (rotated by PA)');
0539 filename = sprintf('%s/Qp1_weighted.png', plotdir);
0540 eval(sprintf('print -dpng %s', filename));
0541 
0542 imagesc(x4label, x4label, squeeze(U_weight));
0543 axis image; colorbar; xlabel('x-offset (degrees)');
0544 ylabel('y-offset (degrees)');
0545 title('Full U');
0546 filename = sprintf('%s/U_weighted.png', plotdir);
0547 eval(sprintf('print -dpng %s', filename));
0548 
0549 imagesc(x4label, x4label, squeeze(Urem_weight));
0550 axis image; colorbar; xlabel('x-offset (degrees)');
0551 ylabel('y-offset (degrees)');
0552 title('Full U - baseline removed');
0553 filename = sprintf('%s/Urem_weighted.png', plotdir);
0554 eval(sprintf('print -dpng %s', filename));
0555 
0556 imagesc(x4label, x4label, squeeze(Up1_weight));
0557 axis image; colorbar; xlabel('x-offset (degrees)');
0558 ylabel('y-offset (degrees)');
0559 title('Full U (rotated by PA)');
0560 filename = sprintf('%s/Up1_weighted.png', plotdir);
0561 eval(sprintf('print -dpng %s', filename));
0562 
0563 
0564 imagesc(x4label, x4label, squeeze(V_weight));
0565 axis image; colorbar; xlabel('x-offset (degrees)');
0566 ylabel('y-offset (degrees)');
0567 title('Full V');
0568 filename = sprintf('%s/V_weighted.png', plotdir);
0569 eval(sprintf('print -dpng %s', filename));
0570 
0571 imagesc(x4label, x4label, squeeze(Vrem_weight));
0572 axis image; colorbar; xlabel('x-offset (degrees)');
0573 ylabel('y-offset (degrees)');
0574 title('Full V - baseline removed');
0575 filename = sprintf('%s/Vrem_weighted.png', plotdir);
0576 eval(sprintf('print -dpng %s', filename));
0577 
0578 imagesc(x4label, x4label, squeeze(I_mean));
0579 axis image; colorbar; xlabel('x-offset (degrees)');
0580 ylabel('y-offset (degrees)');
0581 title('Full I');
0582 filename = sprintf('%s/I_meaned.png', plotdir);
0583 eval(sprintf('print -dpng %s', filename));
0584 
0585 imagesc(x4label, x4label, squeeze(Q_mean));
0586 axis image; colorbar; xlabel('x-offset (degrees)');
0587 ylabel('y-offset (degrees)');
0588 title('Full Q');
0589 filename = sprintf('%s/Q_meaned.png', plotdir);
0590 eval(sprintf('print -dpng %s', filename));
0591 
0592 imagesc(x4label, x4label, squeeze(Qrem_mean));
0593 axis image; colorbar; xlabel('x-offset (degrees)');
0594 ylabel('y-offset (degrees)');
0595 title('Full Q - baseline removed');
0596 filename = sprintf('%s/Qrem_meaned.png', plotdir);
0597 eval(sprintf('print -dpng %s', filename));
0598 
0599 imagesc(x4label, x4label, squeeze(Qp1_mean));
0600 axis image; colorbar; xlabel('x-offset (degrees)');
0601 ylabel('y-offset (degrees)');
0602 title('Full Q (rotated by PA)');
0603 filename = sprintf('%s/Qp1_meaned.png', plotdir);
0604 eval(sprintf('print -dpng %s', filename));
0605 
0606 imagesc(x4label, x4label, squeeze(U_mean));
0607 axis image; colorbar; xlabel('x-offset (degrees)');
0608 ylabel('y-offset (degrees)');
0609 title('Full U');
0610 filename = sprintf('%s/U_meaned.png', plotdir);
0611 eval(sprintf('print -dpng %s', filename));
0612 
0613 imagesc(x4label, x4label, squeeze(Urem_mean));
0614 axis image; colorbar; xlabel('x-offset (degrees)');
0615 ylabel('y-offset (degrees)');
0616 title('Full U - baseline removed');
0617 filename = sprintf('%s/Urem_meaned.png', plotdir);
0618 eval(sprintf('print -dpng %s', filename));
0619 
0620 imagesc(x4label, x4label, squeeze(Up1_mean));
0621 axis image; colorbar; xlabel('x-offset (degrees)');
0622 ylabel('y-offset (degrees)');
0623 title('Full U (rotated by PA)');
0624 filename = sprintf('%s/Up1_meaned.png', plotdir);
0625 eval(sprintf('print -dpng %s', filename));
0626 
0627 
0628 imagesc(x4label, x4label, squeeze(V_mean));
0629 axis image; colorbar; xlabel('x-offset (degrees)');
0630 ylabel('y-offset (degrees)');
0631 title('Full V');
0632 filename = sprintf('%s/V_meaned.png', plotdir);
0633 eval(sprintf('print -dpng %s', filename));
0634 
0635 imagesc(x4label, x4label, squeeze(Vrem_mean));
0636 axis image; colorbar; xlabel('x-offset (degrees)');
0637 ylabel('y-offset (degrees)');
0638 title('Full V - baseline removed');
0639 filename = sprintf('%s/Vrem_meaned.png', plotdir);
0640 eval(sprintf('print -dpng %s', filename));
0641 
0642 
0643 imagesc(x4label, x4label, squeeze(sqrt(Qp1_mean.^2 + Up1_mean.^2)));
0644 axis image; colorbar; xlabel('x-offset (degrees)');
0645 ylabel('y-offset (degrees)');
0646 title('Full Polarization - rotated');
0647 filename = sprintf('%s/FullPol_meaned.png', plotdir);
0648 eval(sprintf('print -dpng %s', filename));
0649 
0650 
0651 figure('visible', 'off')
0652 subplot(2,2,1);
0653 imagesc(x4label,x4label, squeeze(I_mean))
0654 axis image; colorbar; xlabel('x-offset (degrees)');
0655 ylabel('y-offset (degrees)');
0656 title('I');
0657 subplot(2,2,2);
0658 imagesc(x4label,x4label, squeeze(Qp1_mean))
0659 axis image; colorbar; xlabel('x-offset (degrees)');
0660 ylabel('y-offset (degrees)');
0661 title('Q');
0662 subplot(2,2,3);
0663 imagesc(x4label,x4label, squeeze(Up1_mean))
0664 axis image; colorbar; xlabel('x-offset (degrees)');
0665 ylabel('y-offset (degrees)');
0666 title('U');
0667 subplot(2,2,4);
0668 imagesc(x4label,x4label, squeeze(sqrt(Up1_mean.^2 + Qp1_mean.^2)));
0669 axis image; colorbar; xlabel('x-offset (degrees)');
0670 ylabel('y-offset (degrees)');
0671 title('sqrt(Q^2 + U^2)');
0672 filename = sprintf('%s/all_stokes.png', plotdir);
0673 eval(sprintf('print -dpng %s', filename));
0674 
0675 
0676 return;