0001 function source_raster(start, stop, savefile, plotdir, doLeak, model, dcal)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 warning off
0012
0013
0014 if(nargin<7)
0015
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
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 ind = bitand(d.array.frame.features, 2^15)>0;
0035 if(length(find(ind))~=0)
0036 d.array.frame.features(ind) = d.array.frame.features(ind)-2^15+2^9;
0037 end
0038
0039
0040 dcal = reduceData(d, 'reduc/rasterScript.red');
0041
0042
0043 display('Calculating RA/DEC');
0044 long=-118.2822;
0045 lat=37.2339;
0046 az = dcal.antenna0.servo.apparent(:,1);
0047 el = dcal.antenna0.servo.apparent(:,2);
0048 jd=mjd2jd(dcal.antenna0.receiver.utc);
0049 [equa] = horiz_coo([pi/180*(az) pi/180*(el)],jd,[pi/180*(long) ...
0050 pi/180*(lat)],'e');
0051 dcal.antenna0.servo.equa = equa;
0052 end
0053
0054
0055 ind = bitand(dcal.array.frame.features, 2^9) > 0;
0056 [si ei] = findStartStop(ind);
0057
0058 display(sprintf('There are %d rasters in this data set', length(si)));
0059
0060
0061 xgrid = -5.5:0.1:5.5;
0062 ygrid = xgrid;
0063 [X, Y] = meshgrid(xgrid,ygrid);
0064
0065 mapI = nan(length(si), length(xgrid)-1, length(xgrid)-1);
0066 mapIvar = mapI;
0067 mapV = mapI;
0068 mapVrem = mapI;
0069 mapVvar = mapI;
0070 mapVremvar=mapI;
0071 mapQ = mapI;
0072 mapU = mapI;
0073 mapQrem = mapI;
0074 mapUrem = mapI;
0075 mapQvar = mapI;
0076 mapUvar = mapI;
0077 mapQremvar = mapI;
0078 mapUremvar = mapI;
0079
0080 allpa = nan(1,length(si));
0081 for m=1:length(si)
0082 ind = zeros(size(dcal.array.frame.features));
0083 ind(si(m):ei(m)) = 1;
0084 ind = logical(ind);
0085 dcut = framecut(dcal, ind);
0086
0087
0088 da = abs(deriv(dcut.antenna0.servo.az));
0089 ind = da>0.0094 & da<0.0108;
0090 dcut = framecut(dcut, ind);
0091
0092
0093
0094 perGood = length(find(ind))./length(da)*100;
0095
0096 allGoodPer(m) = perGood;
0097
0098 if(perGood > 15)
0099
0100
0101
0102 I1 = dcut.antenna0.receiver.data(:,1);
0103 I2 = dcut.antenna0.receiver.data(:,8);
0104 V = dcut.antenna0.receiver.data(:,1) - dcut.antenna0.receiver.data(:,8);
0105 Q = dcut.antenna0.receiver.data(:,6);
0106 U = dcut.antenna0.receiver.data(:,7);
0107 t = (dcut.antenna0.receiver.utc - dcut.antenna0.receiver.utc(1))*24;
0108
0109
0110
0111
0112
0113 ra = dcut.antenna0.servo.equa(:,1)*180/pi;
0114 dec = dcut.antenna0.servo.equa(:,2)*180/pi;
0115
0116 rafc = mean(ra);
0117 decfc = mean(dec);
0118 aa = I1+ I2;
0119 mm = find(aa == max(aa));
0120 rafc = dcut.antenna0.servo.equa(mm,1)*180/pi;
0121 decfc = dcut.antenna0.servo.equa(mm,2)*180/pi;
0122 [x y] = radec_to_fieldoff(ra, dec, rafc, decfc);
0123
0124
0125 thisra = rafc*pi/180;
0126 thisdec = decfc*pi/180;
0127 lat = mean(dcut.antenna0.tracker.siteActual(:,2))*pi/180;
0128 lst = dcut.antenna0.tracker.lst(ceil(mm/100))*15*pi/180;
0129 ha = (lst - thisra);
0130 cosel_cospa = sin(lat)*cos(thisdec) - sin(thisdec)*cos(lat)*cos(ha);
0131 cosel_sinpa = sin(ha)*cos(lat);
0132 pa = atan2(cosel_sinpa, cosel_cospa);
0133 pa = pa*180/pi;
0134 allpa(m) = pa;
0135
0136
0137
0138 angdist = pyth(x,y);
0139 indfit = angdist > 1;
0140 [slope, b] = linfit(t(indfit), I1(indfit));
0141 I1rem = I1 - (slope*t+b);
0142 [slope, b] = linfit(t(indfit), I2(indfit));
0143 I2rem = I2 - (slope*t+b);
0144 Vrem = I1rem - I2rem;
0145 [slope, b] = linfit(t(indfit), Q(indfit));
0146 Qrem = Q - (slope*t+b);
0147 [slope, b] = linfit(t(indfit), U(indfit));
0148 Urem = U - (slope*t+b);
0149 Itot = I1rem + I2rem;
0150
0151
0152
0153 sigmaFit = 0.73/(2*sqrt(2*log(2)));
0154 LB = [0 -0.5 -0.5 0.1 0.1];
0155 UB = [inf 0.5 0.5 2 2];
0156 dataCell{1} = x;
0157 dataCell{2} = y;
0158 dataCell{3} = [1 1 1 1 1];
0159
0160 dataCell{4} = [max(I1rem) 0 0 sigmaFit sigmaFit];
0161 [fitVals(m,1,:), resnorm, resid1, exitflag, output, lambda, J1] = ...
0162 lsqcurvefit('gauss2D_fit_nograd', [max(I1rem) 0 0 sigmaFit sigmaFit], ...
0163 dataCell, I1rem, LB, UB);
0164
0165 dataCell{4} = [max(I2rem) 0 0 sigmaFit sigmaFit];
0166 [fitVals(m,2,:), resnorm, resid1, exitflag, output, lambda, J1] = ...
0167 lsqcurvefit('gauss2D_fit_nograd', [max(I2rem) 0 0 sigmaFit sigmaFit], ...
0168 dataCell, I2rem, LB, UB);
0169
0170 dataCell{4} = [max(Itot) 0 0 sigmaFit sigmaFit];
0171 [fitVals(m,3,:), resnorm, resid1, exitflag, output, lambda, J1] = ...
0172 lsqcurvefit('gauss2D_fit_nograd', [max(Itot) 0 0 sigmaFit sigmaFit], ...
0173 dataCell, Itot, LB, UB);
0174
0175 dataCell{4} = [median(Qrem) 0 0 sigmaFit sigmaFit];
0176 LB = [-1 -0.5 -0.5 0.1 0.1];
0177 [QfitVals(m,1,:), resnorm, resid1, exitflag, output, lambda, J1] = ...
0178 lsqcurvefit('gauss2D_fit_nograd', [median(Qrem) 0 0 sigmaFit sigmaFit], ...
0179 dataCell, Qrem, LB, UB);
0180
0181 dataCell{4} = [median(Urem) 0 0 sigmaFit sigmaFit];
0182 [UfitVals(m,1,:), resnorm, resid1, exitflag, output, lambda, J1] = ...
0183 lsqcurvefit('gauss2D_fit_nograd', [median(Urem) 0 0 sigmaFit sigmaFit], ...
0184 dataCell, Urem, LB, UB);
0185
0186
0187 dataCell{4} = [median(Qrem) squeeze(fitVals(m,3,2:5))'];
0188 dataCell{3} = [1 1 1 0 0];
0189 [QfitVals(m,2,:), resnorm, resid1, exitflag, output, lambda, J1] = ...
0190 lsqcurvefit('gauss2D_fit_nograd', dataCell{4}, ...
0191 dataCell, Qrem, LB, UB);
0192
0193 dataCell{4} = [median(Urem) squeeze(fitVals(m,3,2:5))'];
0194 dataCell{3} = [1 1 1 0 0];
0195 [UfitVals(m,2,:), resnorm, resid1, exitflag, output, lambda, J1] = ...
0196 lsqcurvefit('gauss2D_fit_nograd', dataCell{4}, ...
0197 dataCell, Urem, LB, UB);
0198
0199
0200 dataCell{4} = [median(Qrem) squeeze(fitVals(m,3,2:5))'];
0201 dataCell{3} = [1 0 0 1 1];
0202 [QfitVals(m,3,:), resnorm, resid1, exitflag, output, lambda, J1] = ...
0203 lsqcurvefit('gauss2D_fit_nograd', dataCell{4}, ...
0204 dataCell, Qrem, LB, UB);
0205
0206 dataCell{4} = [median(Urem) squeeze(fitVals(m,3,2:5))'];
0207 dataCell{3} = [1 0 0 1 1];
0208 [UfitVals(m,3,:), resnorm, resid1, exitflag, output, lambda, J1] = ...
0209 lsqcurvefit('gauss2D_fit_nograd', dataCell{4}, ...
0210 dataCell, Urem, LB, UB);
0211
0212
0213 dataCell{4} = [median(Qrem) squeeze(fitVals(m,3,2:5))'];
0214 dataCell{3} = [1 0 0 0 0];
0215 [QfitVals(m,4,:), resnorm, resid1, exitflag, output, lambda, J1] = ...
0216 lsqcurvefit('gauss2D_fit_nograd', dataCell{4}, ...
0217 dataCell, Qrem, LB, UB);
0218
0219 dataCell{4} = [median(Urem) squeeze(fitVals(m,3,2:5))'];
0220 dataCell{3} = [1 0 0 0 0];
0221 [UfitVals(m,4,:), resnorm, resid1, exitflag, output, lambda, J1] = ...
0222 lsqcurvefit('gauss2D_fit_nograd', dataCell{4}, ...
0223 dataCell, Urem, LB, UB);
0224
0225
0226 if(abs(mean(fitVals(m,:,2))) < 0.3)
0227 xoff(m) = mean(fitVals(m,:,2));
0228 else
0229 xoff(m) = xoff(m-1);
0230 end
0231 if(abs(mean(fitVals(m,:,3))) < 0.3)
0232 yoff(m) = mean(fitVals(m,:,3));
0233 else
0234 yoff(m) = yoff(m-1);
0235 end
0236
0237
0238
0239 x = x-xoff(m);
0240 y = y-yoff(m);
0241
0242
0243 mapI1 = nan(length(xgrid)-1, length(xgrid)-1);
0244 mapI2 = nan(length(xgrid)-1, length(xgrid)-1);
0245 mapV1 = nan(length(xgrid)-1, length(xgrid)-1);
0246 mapV1rem= nan(length(xgrid)-1, length(xgrid)-1);
0247 mapQ1 = nan(length(xgrid)-1, length(xgrid)-1);
0248 mapU1 = nan(length(xgrid)-1, length(xgrid)-1);
0249
0250
0251 mapI1var = nan(length(xgrid)-1, length(xgrid)-1);
0252 mapI2var = nan(length(xgrid)-1, length(xgrid)-1);
0253 mapV1var = nan(length(xgrid)-1, length(xgrid)-1);
0254 mapV1remvar = nan(length(xgrid)-1, length(xgrid)-1);
0255 mapQ1var = nan(length(xgrid)-1, length(xgrid)-1);
0256 mapU1var = nan(length(xgrid)-1, length(xgrid)-1);
0257
0258 x4label = xgrid+0.1/2;
0259 x4label(length(x4label)) = [];
0260 for mm=2:length(xgrid)
0261 if(mm==2)
0262 fx = x>=xgrid(1) & x<xgrid(2);
0263 else
0264 fx = x>=xgrid(mm-1) & x<xgrid(mm);
0265 end
0266 for mmm=2:length(ygrid)
0267 if(mmm==2)
0268 fy = y>=ygrid(1) & y<ygrid(2);
0269 else
0270 fy = y>=ygrid(mmm-1) & y<ygrid(mmm);
0271 end
0272 f = find(fx&fy);
0273 mapI1(mm-1,mmm-1) = nanmean(I1rem(f));
0274 mapI2(mm-1,mmm-1) = nanmean(I2rem(f));
0275 mapV1(mm-1,mmm-1) = nanmean(V(f));
0276 mapV1rem(mm-1,mmm-1) = nanmean(Vrem(f));
0277 mapQ1(mm-1,mmm-1) = nanmean(Q(f));
0278 mapU1(mm-1,mmm-1) = nanmean(U(f));
0279 mapQ1rem(mm-1,mmm-1) = nanmean(Qrem(f));
0280 mapU1rem(mm-1,mmm-1) = nanmean(Urem(f));
0281 mapI1var(mm-1,mmm-1) = nanvar(I1rem(f));
0282 mapI2var(mm-1,mmm-1) = nanvar(I2rem(f));
0283 mapV1var(mm-1,mmm-1) = nanvar(V(f));
0284 mapV1remvar(mm-1,mmm-1) = nanvar(Vrem(f));
0285 mapQ1var(mm-1,mmm-1) = nanvar(Q(f));
0286 mapU1var(mm-1,mmm-1) = nanvar(U(f));
0287 mapQ1remvar(mm-1,mmm-1) = nanvar(Qrem(f));
0288 mapU1remvar(mm-1,mmm-1) = nanvar(Urem(f));
0289 end
0290 end
0291
0292
0293 mapI(m,:,:) = mapI1 + mapI2;
0294 mapIvar(m,:,:) = mapI1var + mapI2var;
0295 mapV(m,:,:) = mapV1;
0296 mapVvar(m,:,:) = mapV1var;
0297 mapVrem(m,:,:) = mapV1rem;
0298 mapVremvar(m,:,:) = mapV1remvar;
0299 mapQ(m,:,:) = mapQ1;
0300 mapU(m,:,:) = mapU1;
0301 mapQvar(m,:,:) = mapQ1var;
0302 mapUvar(m,:,:) = mapU1var;
0303 mapQrem(m,:,:) = mapQ1rem;
0304 mapUrem(m,:,:) = mapU1rem;
0305 mapQremvar(m,:,:) = mapQ1remvar;
0306 mapUremvar(m,:,:) = mapU1remvar;
0307
0308 figure('visible', 'off')
0309 imagesc(x4label, x4label, squeeze(mapI(m,:,:)));
0310 axis image; colorbar; xlabel('x-offset (degrees)');
0311 ylabel('y-offset (degrees)');
0312 txt = sprintf('I, pa=%3.2f, Scan #%d of %d', pa, m, length(si) );
0313 title(txt);
0314 filename = sprintf('%s/I_%02d.png', plotdir, m);
0315 eval(sprintf('print -dpng %s', filename));
0316
0317
0318
0319 imagesc(x4label, x4label, squeeze(mapQ(m,:,:)));
0320 axis image; colorbar; xlabel('x-offset (degrees)');
0321 ylabel('y-offset (degrees)');
0322 txt = sprintf('Q, pa=%3.2f, Scan #%d of %d', pa, m, length(si) );
0323 title(txt);
0324 filename = sprintf('%s/Q_%02d.png', plotdir, m);
0325 eval(sprintf('print -dpng %s', filename));
0326
0327 imagesc(x4label, x4label, squeeze(mapU(m,:,:)));
0328 axis image; colorbar; xlabel('x-offset (degrees)');
0329 ylabel('y-offset (degrees)');
0330 txt = sprintf('U, pa=%3.2f, Scan #%d of %d', pa, m, length(si) );
0331 title(txt);
0332 filename = sprintf('%s/U_%02d.png', plotdir, m);
0333 eval(sprintf('print -dpng %s', filename));
0334
0335 imagesc(x4label, x4label, squeeze(mapQrem(m,:,:)));
0336 axis image; colorbar; xlabel('x-offset (degrees)');
0337 ylabel('y-offset (degrees)');
0338 txt = sprintf('Q baseline removed, pa=-%3.2f, Scan #%d of %d', pa, m, length(si) );
0339 title(txt);
0340 filename = sprintf('%s/Qrem_%02d.png', plotdir, m);
0341 eval(sprintf('print -dpng %s', filename));
0342
0343 imagesc(x4label, x4label, squeeze(mapUrem(m,:,:)));
0344 axis image; colorbar; xlabel('x-offset (degrees)');
0345 ylabel('y-offset (degrees)');
0346 txt = sprintf('U baseline removed, pa=-%3.2f, Scan #%d of %d', pa, m, length(si) );
0347 title(txt);
0348 filename = sprintf('%s/Urem_%02d.png', plotdir, m);
0349 eval(sprintf('print -dpng %s', filename));
0350
0351 imagesc(x4label, x4label, squeeze(mapV(m,:,:)));
0352 axis image; colorbar; xlabel('x-offset (degrees)');
0353 ylabel('y-offset (degrees)');
0354 txt = sprintf('V, pa=%3.2f, Scan #%d of %d', pa, m, length(si) );
0355 title(txt);
0356 filename = sprintf('%s/V_%02d.png', plotdir, m);
0357 eval(sprintf('print -dpng %s', filename));
0358
0359 imagesc(x4label, x4label, squeeze(mapVrem(m,:,:)));
0360 axis image; colorbar; xlabel('x-offset (degrees)');
0361 ylabel('y-offset (degrees)');
0362 txt = sprintf('V, baseline removed, pa=%3.2f, Scan #%d of %d', pa, m, length(si) );
0363 title(txt);
0364 filename = sprintf('%s/Vrem_%02d.png', plotdir, m);
0365 eval(sprintf('print -dpng %s', filename));
0366
0367 end
0368 end
0369
0370 eval(sprintf('save %s/%s', plotdir, savefile));
0371
0372
0373
0374
0375
0376
0377
0378
0379 for m=1:length(si)
0380 aa = mapIvar(m,:,:);
0381 medval = nanmedian(aa(:));
0382 maxval = nanmax(aa(:));
0383 f = find(mapIvar(m,:,:)==0);
0384 mapIvar(m,f) = mean([maxval, medval]);
0385
0386 aa = mapQvar(m,:,:);
0387 maxval = nanmax(aa(:));
0388 f = find(aa == 0);
0389 mapQvar(m,f) = maxval;
0390
0391 aa = mapQremvar(m,:,:);
0392 maxval = nanmax(aa(:));
0393 f = find(aa == 0);
0394 mapQremvar(m,f) = maxval;
0395
0396 aa = mapUvar(m,:,:);
0397 maxval = nanmax(aa(:));
0398 f = find(aa == 0);
0399 mapUvar(m,f) = maxval;
0400
0401 aa = mapUremvar(m,:,:);
0402 maxval = nanmax(aa(:));
0403 f = find(aa == 0);
0404 mapUremvar(m,f) = maxval;
0405
0406 aa = mapVvar(m,:,:);
0407 maxval = nanmax(aa(:));
0408 f = find(aa == 0);
0409 mapVvar(m,f) = maxval;
0410
0411 aa = mapVremvar(m,:,:);
0412 maxval = nanmax(aa(:));
0413 f = find(aa == 0);
0414 mapVremvar(m,f) = maxval;
0415 end
0416
0417
0418 [I_weight I_weight_var] = vwstat(mapI, mapIvar,1);
0419 [Q_weight Q_weight_var] = vwstat(mapQ, mapQvar, 1);
0420 [Qrem_weight Qrem_weight_var] = vwstat(mapQrem, mapQremvar, 1);
0421 [U_weight U_weight_var] = vwstat(mapU, mapUvar, 1);
0422 [Urem_weight Urem_weight_var] = vwstat(mapUrem, mapUremvar, 1);
0423 [V_weight V_weight_var] = vwstat(mapV, mapVvar, 1);
0424 [Vrem_weight Vrem_weight_var] = vwstat(mapVrem, mapVremvar, 1);
0425
0426
0427 I_mean = nanmean(mapI,1);
0428 Q_mean = nanmean(mapQ,1);
0429 Qrem_mean = nanmean(mapQrem, 1);
0430 U_mean = nanmean(mapU,1);
0431 Urem_mean = nanmean(mapUrem, 1);
0432 V_mean = nanmean(mapV,1);
0433 Vrem_mean = nanmean(mapVrem,1);
0434
0435
0436
0437 imagesc(x4label, x4label, squeeze(I_weight));
0438 axis image; colorbar; xlabel('x-offset (degrees)');
0439 ylabel('y-offset (degrees)');
0440 title('Full I');
0441 filename = sprintf('%s/I_weighted.png', plotdir);
0442 eval(sprintf('print -dpng %s', filename));
0443
0444 imagesc(x4label, x4label, squeeze(Q_weight));
0445 axis image; colorbar; xlabel('x-offset (degrees)');
0446 ylabel('y-offset (degrees)');
0447 title('Full Q');
0448 filename = sprintf('%s/Q_weighted.png', plotdir);
0449 eval(sprintf('print -dpng %s', filename));
0450
0451 imagesc(x4label, x4label, squeeze(Qrem_weight));
0452 axis image; colorbar; xlabel('x-offset (degrees)');
0453 ylabel('y-offset (degrees)');
0454 title('Full Q - baseline removed');
0455 filename = sprintf('%s/Qrem_weighted.png', plotdir);
0456 eval(sprintf('print -dpng %s', filename));
0457
0458 imagesc(x4label, x4label, squeeze(U_weight));
0459 axis image; colorbar; xlabel('x-offset (degrees)');
0460 ylabel('y-offset (degrees)');
0461 title('Full U');
0462 filename = sprintf('%s/U_weighted.png', plotdir);
0463 eval(sprintf('print -dpng %s', filename));
0464
0465 imagesc(x4label, x4label, squeeze(Urem_weight));
0466 axis image; colorbar; xlabel('x-offset (degrees)');
0467 ylabel('y-offset (degrees)');
0468 title('Full U - baseline removed');
0469 filename = sprintf('%s/Urem_weighted.png', plotdir);
0470 eval(sprintf('print -dpng %s', filename));
0471
0472
0473 imagesc(x4label, x4label, squeeze(V_weight));
0474 axis image; colorbar; xlabel('x-offset (degrees)');
0475 ylabel('y-offset (degrees)');
0476 title('Full V');
0477 filename = sprintf('%s/V_weighted.png', plotdir);
0478 eval(sprintf('print -dpng %s', filename));
0479
0480 imagesc(x4label, x4label, squeeze(Vrem_weight));
0481 axis image; colorbar; xlabel('x-offset (degrees)');
0482 ylabel('y-offset (degrees)');
0483 title('Full V - baseline removed');
0484 filename = sprintf('%s/Vrem_weighted.png', plotdir);
0485 eval(sprintf('print -dpng %s', filename));
0486
0487
0488 imagesc(x4label, x4label, squeeze(I_mean));
0489 axis image; colorbar; xlabel('x-offset (degrees)');
0490 ylabel('y-offset (degrees)');
0491 title('Full I');
0492 filename = sprintf('%s/I_meaned.png', plotdir);
0493 eval(sprintf('print -dpng %s', filename));
0494
0495 imagesc(x4label, x4label, squeeze(Q_mean));
0496 axis image; colorbar; xlabel('x-offset (degrees)');
0497 ylabel('y-offset (degrees)');
0498 title('Full Q');
0499 filename = sprintf('%s/Q_meaned.png', plotdir);
0500 eval(sprintf('print -dpng %s', filename));
0501
0502 imagesc(x4label, x4label, squeeze(Qrem_mean));
0503 axis image; colorbar; xlabel('x-offset (degrees)');
0504 ylabel('y-offset (degrees)');
0505 title('Full Q - baseline removed');
0506 filename = sprintf('%s/Qrem_meaned.png', plotdir);
0507 eval(sprintf('print -dpng %s', filename));
0508
0509 imagesc(x4label, x4label, squeeze(U_mean));
0510 axis image; colorbar; xlabel('x-offset (degrees)');
0511 ylabel('y-offset (degrees)');
0512 title('Full U');
0513 filename = sprintf('%s/U_meaned.png', plotdir);
0514 eval(sprintf('print -dpng %s', filename));
0515
0516 imagesc(x4label, x4label, squeeze(Urem_mean));
0517 axis image; colorbar; xlabel('x-offset (degrees)');
0518 ylabel('y-offset (degrees)');
0519 title('Full U - baseline removed');
0520 filename = sprintf('%s/Urem_meaned.png', plotdir);
0521 eval(sprintf('print -dpng %s', filename));
0522
0523
0524 imagesc(x4label, x4label, squeeze(V_mean));
0525 axis image; colorbar; xlabel('x-offset (degrees)');
0526 ylabel('y-offset (degrees)');
0527 title('Full V');
0528 filename = sprintf('%s/V_meaned.png', plotdir);
0529 eval(sprintf('print -dpng %s', filename));
0530
0531 imagesc(x4label, x4label, squeeze(Vrem_mean));
0532 axis image; colorbar; xlabel('x-offset (degrees)');
0533 ylabel('y-offset (degrees)');
0534 title('Full V - baseline removed');
0535 filename = sprintf('%s/Vrem_meaned.png', plotdir);
0536 eval(sprintf('print -dpng %s', filename));
0537
0538
0539
0540
0541
0542
0543
0544
0545 for m=1:size(mapI,1)
0546 Imid(m) = nanmean(nanmean(mapI(m,55:56,55:56)));
0547 Qmid(m) = (nanmean(nanmean(mapQrem(m,55:56,55:56))));
0548 Umid(m) = (nanmean(nanmean(mapUrem(m,55:56,55:56))));
0549 end
0550
0551 poldata = [Imid;Qmid;Umid;allpa];
0552 poldata = poldata';
0553 poldata(isnan(poldata(:,4)),:) = [];
0554 poldata(poldata(:,1)<0.5,:) = [];
0555 model = [615 0.07 150];
0556 calmat_map = polcal(poldata, model);
0557
0558
0559
0560
0561 Ifit = fitVals(:,3,1);
0562 Qfit = QfitVals(:,1,1);
0563 Ufit = UfitVals(:,1,1);
0564 poldata = [Ifit Qfit Ufit allpa'];
0565 poldata(isnan(poldata(:,4)),:) = [];
0566 poldata(poldata(:,1)<0.5,:) = [];
0567 calmat_fit1 = polcal(poldata, model);
0568
0569
0570 Ifit = fitVals(:,3,1);
0571 Qfit = QfitVals(:,2,1);
0572 Ufit = UfitVals(:,2,1);
0573 poldata = [Ifit Qfit Ufit allpa'];
0574 poldata(isnan(poldata(:,4)),:) = [];
0575 poldata(poldata(:,1)<0.5,:) = [];
0576 calmat_fit2 = polcal(poldata, model);
0577
0578
0579 Ifit = fitVals(:,3,1);
0580 Qfit = QfitVals(:,3,1);
0581 Ufit = UfitVals(:,3,1);
0582 poldata = [Ifit Qfit Ufit allpa'];
0583 poldata(isnan(poldata(:,4)),:) = [];
0584 poldata(poldata(:,1)<0.5,:) = [];
0585 calmat_fit3 = polcal(poldata, model);
0586
0587
0588
0589 Ifit = fitVals(:,3,1);
0590 Qfit = QfitVals(:,4,1);
0591 Ufit = UfitVals(:,4,1);
0592 poldata = [Ifit Qfit Ufit allpa'];
0593 poldata(isnan(poldata(:,4)),:) = [];
0594 poldata(poldata(:,1)<0.5,:) = [];
0595 calmat_fit4 = polcal(poldata, model);
0596
0597 eval(sprintf('save %s/%s', plotdir, savefile));
0598
0599
0600
0601
0602 newsavefile = strcat(savefile, '_leak_map');
0603 newplotdir = strcat(plotdir, '/leak/map');
0604 source_raster_remLeak_v2([],[], calmat_map, newsavefile, newplotdir, dcal);
0605
0606 newsavefile = strcat(savefile, '_leak_fit1');
0607 newplotdir = strcat(plotdir, '/leak/fit1');
0608 source_raster_remLeak_v2([],[], calmat_fit1, newsavefile, newplotdir, dcal);
0609
0610 newsavefile = strcat(savefile, '_leak_fit2');
0611 newplotdir = strcat(plotdir, '/leak/fit2');
0612 source_raster_remLeak_v2([],[], calmat_fit2, newsavefile, newplotdir, dcal);
0613
0614 newsavefile = strcat(savefile, '_leak_fit3');
0615 newplotdir = strcat(plotdir, '/leak/fit3');
0616 source_raster_remLeak_v2([],[], calmat_fit3, newsavefile, newplotdir, dcal);
0617
0618 newsavefile = strcat(savefile, '_leak_fit4');
0619 newplotdir = strcat(plotdir, '/leak/fit4');
0620 source_raster_remLeak_v2([],[], calmat_fit4, newsavefile, newplotdir, dcal);
0621
0622 calmat_act = [554.8911 20.3628 -22.6129; -10.3716 599.9700 547.3359; ...
0623 -157.4338 451.8572 -616.9229];
0624
0625 newsavefile = strcat(savefile, '_leak_act');
0626 newplotdir = strcat(plotdir, '/leak/act');
0627 source_raster_remLeak_v2([],[], calmat_act, newsavefile, newplotdir, dcal);
0628
0629 return;