Description of source

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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
