calcTnoise_v2

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

  function tnoiseVals = calcTnoise(d)

    This is a test method to calculate the noise diode temperature and the
    opacity from our calibrator observations.  

  sjcm

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

0001 function tnoise = calcTnoise(d)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 %  function tnoiseVals = calcTnoise(d)
0006 %
0007 %    This is a test method to calculate the noise diode temperature and the
0008 %    opacity from our calibrator observations.
0009 %
0010 %  sjcm
0011 %
0012 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0013 
0014 % constants!
0015 % bandwidths.  have to change when we bump up the gain on LNA3 again.
0016 if(d.array.frame.utc > date2mjd(2011, 06, 01))
0017   bw = [763 696]*1e6;
0018 else
0019   bw = [1000 1000]*1e6;
0020 end
0021 Tcmb = 2.728;
0022 
0023 % we need to know if we are observing any sources.
0024 sources = unique(d.antenna0.tracker.source);
0025 for m=1:length(sources)
0026   [sourceNum(m) sourceFlux(m)] = sourceCorrespondance(sources(m));
0027 end
0028 if(all(sourceNum==99))
0029   display('calcTnoise_v2:: No Calibrator Source Observations in this schedule');
0030   display('calcTnoise_v2:: Can not calculate system temperatures');
0031   tnoise = [];
0032   return;
0033 end
0034 
0035 % calculate the pointing corrections for our calibrator observations
0036 if(issubfield(d, 'correction', 'allpoint'))
0037   [pointoff1 pointoff2] = calculatePointingCorrection(d);
0038 else
0039   display('calcTnoise_v2:: No pointing corrections found for your sources');
0040   display('calcTnoise_v2:: Can not calculate the noise diode temperature');
0041   tnoise = [];
0042   return;
0043 end
0044 
0045 % calculate the slope of the skydip
0046 if(~isempty(find(d.index.skydip.slow)))
0047   opacityInfo = calculateOpacityInfo(d);
0048 else
0049   display('calcTnoise_v2:: No skydips found in your data');
0050   display('calcTnoise_v2:: Can not continue calculation');
0051   tnoise = [];
0052   return;
0053 end
0054 
0055 % now we get all the power measurements on our calibrators
0056 cal = obtainPowerMeasurements(d);
0057 % check if we actually have any source
0058 % 18-Mar-2011 (MAS) -- Modified, this crashed if structct() rendered cal no
0059 % longer a structure (which is possible!).
0060 if(~isstruct(cal) || length(cal.whatSource)<1)
0061   display('calcTnoise_v2:: No Calibrator Source Observations in this schedule');
0062   display('calcTnoise_v2:: Can not calculate system temperatures');
0063   tnoise = [];
0064   return;
0065 end
0066 
0067 % now we want the source information:  y-factors, sourcename, and elevation
0068 cal2 = calculateYFactors(cal, d.correction.point_out);
0069 
0070 % pointing offset
0071 cal2 = calcSourceTemp(cal2, pointoff2);
0072 
0073 % next, we match up the opacity info with each source as well.
0074 cal2 = matchOpacity(cal2, opacityInfo);
0075 
0076 if(~all(isnan(cal2.slope)))
0077   % next we can calculate tau from the skydips and srcs
0078   tau  =  calcTauFromSrcDips(cal2, d);
0079 end
0080 
0081   
0082 if(all(isnan(cal2.Tamb)))
0083   display('calcTnoise_v2:: No unflagged Opacity values in your data');
0084   display('calcTnoise_v2:: Can not calculate noise diode');
0085   tnoise = [];
0086   return;
0087 end
0088 
0089 % and now, we are equipped with slope, Tamb, Tcmb, ycrazy, elevation, Tsrc
0090 % (all in cal2), and we can solve for opacity and TND
0091 % solver sucks!  we can do analytic solution.
0092 %  1.  tauo = log(y*Tsrc/Tnd)*sin(el_src)
0093 %  2.  tauo = slope*Tnd/(Tamb-Tcmb)
0094 
0095 % this is method 2.1 in abscal_v2.pdf
0096 
0097 % FIRST WAY:  USES ON/OFF INTEGRATIONS TO MEASURE PSRC-POFF
0098 for m=1:length(cal2.source)
0099   for mm=1:2
0100     thiscal.ycrazy = cal2.ycrazy(m,mm);
0101     thiscal.slope  = cal2.slope(m,mm);
0102     thiscal.Tamb   = cal2.Tamb(m);
0103     thiscal.Tcmb   = Tcmb;
0104     thiscal.elev   = cal2.elev(m);
0105     thiscal.Tsrc   = cal2.Tsrc(m);
0106     x0 = [0.01 3];
0107 
0108     alpha = 2*thiscal.slope./(thiscal.Tamb - thiscal.Tcmb); 
0109     beta  = thiscal.ycrazy.*thiscal.Tsrc.*sin(thiscal.elev);
0110     gamma = sin(thiscal.elev);
0111     
0112     tnd  =  (-gamma + sqrt(gamma.^2+4*alpha*beta))./(2*alpha);
0113     tauo1  = tnd.*alpha;
0114     tauo2  = beta./tnd - gamma;
0115     tau01(m,mm) = tauo1;
0116     Tnd(m,mm)  = tnd;
0117   end
0118 end
0119 cal2.calcTau = tau01;
0120 cal2.Tnd     = Tnd;
0121 cal2.point   = pointoff2;
0122 cal2.opacity = opacityInfo;
0123 cal2.tauSrcDip = tau;
0124 
0125 
0126 
0127 % second way:
0128 % USES THE POINTING CROSSES TO GET A MEASURE OF PSRC-POFF
0129 for m=1:length(cal2.source)
0130   for mm=1:2
0131     thiscal.ycrazy = cal2.ycrazy2(m,mm);
0132     thiscal.slope  = cal2.slope(m,mm);
0133     thiscal.Tamb   = cal2.Tamb(m);
0134     thiscal.Tcmb   = Tcmb;
0135     thiscal.elev   = cal2.elev(m);
0136     thiscal.Tsrc   = cal2.Tsrc_o(m);
0137     x0 = [0.01 3];
0138 
0139     alpha = 2*thiscal.slope./(thiscal.Tamb - thiscal.Tcmb); 
0140     beta  = thiscal.ycrazy.*thiscal.Tsrc.*sin(thiscal.elev);
0141     gamma = sin(thiscal.elev);
0142     
0143     tnd  =  (-gamma + sqrt(gamma.^2+4*alpha*beta))./(2*alpha);
0144     tauo1  = tnd.*alpha;
0145     tauo2  = beta./tnd - gamma;
0146     tau01(m,mm) = tauo1;
0147     Tnd(m,mm)  = tnd;
0148   end
0149 end
0150 cal2.calcTau2 = tau01;
0151 cal2.Tnd2     = Tnd;
0152 cal2.point2   = pointoff2;
0153 cal2.opacity2 = opacityInfo;
0154 cal2.tauSrcDip2 = tau;
0155 
0156 
0157 
0158 % need to calculate the errors
0159 cal2.TndErr  = zeros(size(cal2.Tnd));
0160 cal2.calcTauErr = zeros(size(cal2.calcTau));
0161 
0162 tnoise = cal2;
0163 
0164 return;
0165 
0166 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0167 %
0168 %  SUB-FUNCTIONS
0169 %
0170 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0171 
0172   
0173 
0174 function [pointoff1 pointoff2] = calculatePointingCorrection(d)
0175 
0176 % just takes our data and puts it in the format we want for the calculations
0177 
0178 % FWHM = 2sqrt(2ln(2)) call it beta
0179 beta = 0.73/(2*sqrt(2*log(2)));
0180 
0181 % old way:
0182 point = d.correction.pointing;
0183 pointoff1 = [];
0184 while(~isempty(point))
0185   dt = (point(:,1) - point(1,1))*24*60;
0186   ind = find(dt<6 & point(:,2) == point(1,2));
0187   angoff = spaceangle(point(ind,3), point(ind,4), ...
0188       point(ind,3)+point(ind,5), point(ind,4)+point(ind,6), 'deg');
0189   ptmean = nanmean(point(ind,:),1);
0190   pointoff1 = [pointoff1; [ptmean([1 2]) nanmean(angoff) nanstd(angoff)]];
0191   point(ind,:) = [];
0192 end
0193 correction = gauss([1 0 beta], pointoff1(:,3));
0194 pointoff1(:,5) = correction;
0195 
0196 % new way
0197 xoff      = d.correction.allpoint(:, [7 8]);
0198 yoff      = d.correction.allpoint(:, [9 10]);
0199 badpoint  = logical(d.correction.allpoint(:, [5 6]));
0200 time      = d.correction.allpoint(:, 1);
0201 source    = d.correction.allpoint(:, 2);
0202 errx      = d.correction.allpoint(:, [11 12]);
0203 erry      = d.correction.allpoint(:, [13 14]);
0204 %
0205 xoff(badpoint) = nan;
0206 yoff(badpoint) = nan;
0207 errx(badpoint) = inf;
0208 erry(badpoint) = inf;
0209 [xoff vxoff]   = vwstat(xoff, errx.^2, 2);
0210 [yoff vyoff]   = vwstat(yoff, erry.^2, 2);
0211 totaloff       = pyth(xoff, yoff);
0212 totalerroff    = pyth(sqrt(vxoff), sqrt(vyoff));
0213 % next, take a mean of values for the same source, within a short amount
0214 % of time.
0215 pointoff2 = [];
0216 while(~isempty(time))
0217   dt = (time(:,1) - time(1,1))*24*60;
0218   ind = find(dt<6 & source(:) == source(1));
0219   [ptmean pterr] =  vwstat(totaloff(ind), totalerroff(ind));
0220   pointoff2 = [pointoff2; [mean(time(ind)) mean(source(ind)) ptmean ...
0221       pterr]];
0222   time(ind) = [];
0223   source(ind) = [];
0224   totaloff(ind) = [];
0225   totalerroff(ind) = [];
0226 end
0227 correction = gauss([1 0 beta], pointoff2(:,3));
0228 pointoff2(:,5) = correction;
0229 
0230 return;
0231 
0232   
0233 
0234 function opacityInfo = calculateOpacityInfo(d)
0235 
0236 % calculates what we need from the opacity, namely the slope and Tamb
0237 
0238 [si ei] = findStartStop(d.index.skydip.slow);
0239 
0240 for m=1:length(si)
0241   ind = zeros(size(d.array.frame.features));
0242   ind(si(m):ei(m)) = 1;
0243   ind = logical(ind);
0244   
0245   dc = framecut(d, ind);
0246   time(m,1)    = mean(dc.array.frame.utc);
0247   Tamb(m,1)    = mean(dc.array.weather.airTemperature + 273.15);
0248   TambErr(m,1) = std(dc.array.weather.airTemperature + 273.15);
0249 
0250   % next to calculate the slope, let us just get the data where we are
0251   % actually scanning in elevation % scan down at 0.5deg/s
0252   del = deriv(dc.antenna0.servo.el);
0253   ind = del < -0.4/100;
0254   elvals = dc.antenna0.servo.el(ind);
0255   Ivals  = dc.antenna0.receiver.data(ind, [1 8]);
0256 
0257 
0258   thisflag = [ 0 0 ];
0259   indflag = dc.flags.fast(ind,[1 3]);
0260   for mm=1:2
0261     thisInd = ~indflag(:,mm);
0262     elfits = 1./sind(elvals(thisInd));
0263     Ifit   = Ivals(thisInd,mm);
0264     if(length(Ifit) < 10)
0265       mx = nan;
0266       b = nan;
0267       thisflag(mm) = 1;
0268     else
0269       % fit for slope of the line
0270       [mx b] = linfit(elfits, Ifit);
0271     end
0272     slopes(m,mm) = mx;
0273     if( (max(elfits) - min(elfits)) < 0.12)
0274       thisflag(mm) = 1;
0275     end
0276     if(min(elfits) > 60)
0277       thisflag(mm) = 1;
0278     end
0279     % flag if the fit residuals are bad too
0280     fitline = mx.*elfits + b;
0281     resid   = Ifit - fitline;
0282 
0283     % if the receiver is balanced, we expect the I values to be near zero,
0284     % then we just test for a straight-up value in the residual.
0285     % otherwise, we test for a ratio.
0286     if(abs(mean(Ivals(:,mm)))<0.5)
0287       % check the absolute residuals
0288       if(mean(abs(resid))*100>3)
0289     thisflag(mm) = 1;
0290       end
0291       if(rms(resid)*100>3)
0292     thisflag(mm) = 1;
0293       end
0294     else
0295       % then we test for a ratio of rms. (old way)
0296       if(rms(resid)*100>1)
0297     residrat= abs(rms(resid)./mean(Ivals)*100);
0298     if(residrat>0.5)
0299       thisflag(mm) = 1;
0300     end
0301       end
0302     end
0303   end
0304   flag(m,:) = thisflag;
0305 end
0306 
0307 opacityInfo = [time slopes Tamb TambErr flag];
0308 
0309 return;
0310 
0311 
0312   
0313 function cal = obtainPowerMeasurements(d)
0314 
0315 % splits our data up, and does the proper averaging
0316 numNoiseObs = size(d.correction.alpha.indices,1);
0317 cal.whatSource = zeros(numNoiseObs,1);
0318 cal.onSource   = zeros(numNoiseObs,1);
0319 cal.time       = zeros(numNoiseObs,1);
0320 cal.tau        = zeros(numNoiseObs,1);
0321 cal.sourceFlux = zeros(numNoiseObs,1);
0322 indices = d.correction.alpha.indices;
0323 cols = [1 8 6 7];
0324 
0325 for m=1:numNoiseObs
0326   ind = zeros(size(d.antenna0.receiver.utc));
0327   ind(indices(m,1):indices(m,3)) = 1;
0328   ind = logical(ind);
0329   dc = framecut(d, ind);
0330   % we should have 4 intensity sky values, states on/off for chans 1/2
0331   
0332   % apply the flags
0333   dataNoiseOn = d.antenna0.receiver.data((indices(m,4)+5):indices(m,2), cols);
0334   dataNoiseOff= d.antenna0.receiver.data([indices(m,1):(indices(m,4)-5) ...
0335     (indices(m,2)+5):indices(m,3)], cols);
0336 
0337   flagsOn  = d.flags.fast((indices(m,4)+5):indices(m,2), [1 3 2 2]);
0338   flagsOff = d.flags.fast([indices(m,1):(indices(m,4)-5) ...
0339     (indices(m,2)+5):indices(m,3)], [1 3 2 2]);
0340 
0341   dataNoiseOn(flagsOn>0) = nan;
0342   dataNoiseOff(flagsOff>0) = nan;
0343 
0344   % calculate the opacities from theory
0345   tau = calcOpacity(5, dc.array.weather.airTemperature+273.15, ...
0346       dc.array.weather.relativeHumidity, ...
0347       length(dc.array.weather.relativeHumidity), 1);
0348   cal.tau(m,1) = mean(tau);
0349   
0350   % make sure this obs is on a calibrator
0351   thisSource = unique(dc.antenna0.tracker.source);
0352   if(length(thisSource)>1)
0353     cal.whatSource(m,1) = 99;
0354     cal.sourceFlux(m,1) = nan;
0355     cal.fluxErr(m,1)    = nan;
0356     cal.isoff(m,1)      = 0;
0357   else
0358     [cal.whatSource(m,1) cal.sourceFlux(m,1) cal.fluxErr(m,1)] = ...
0359     sourceCorrespondance(thisSource, dc.array.frame.utc(1));
0360     if(isempty(strfind(thisSource{1}, 'ff')))
0361       cal.isoff(m,1) = 0;
0362     else
0363       cal.isoff(m,1) = 1;
0364     end
0365   end
0366   
0367   if(cal.whatSource(m) ~= 99)
0368     % we should now determine whether it is on or off the source.
0369     offsets = unique(dc.antenna0.tracker.horiz_off(:,1));
0370     if(offsets==0)
0371       cal.onSource(m,1) = 1;
0372     else
0373       cal.onSource(m,1) = 0;
0374     end
0375     cal.time(m,1) = dc.antenna0.receiver.utc(1);
0376     cal.PnoiseOn(m,:)  = nanmean(dataNoiseOn);
0377     cal.PnoiseOff(m,:) = nanmean(dataNoiseOff);
0378     cal.varNoiseOn(m,:) = nanvar(dataNoiseOn);
0379     cal.varNoiseOff(m,:) = nanvar(dataNoiseOff);
0380     cal.datanoiseon{m,1} = dataNoiseOn;
0381     cal.datanoiseoff{m,1} = dataNoiseOff;
0382   else
0383     cal.onSource(m,1) = 0;
0384     cal.PnoiseOn(m,1:4) = nan;
0385     cal.PnoiseOff(m,1:4) = nan;
0386     cal.varNoiseOn(m,1:4) = nan;
0387     cal.varNoiseOff(m,1:4) = nan;
0388     cal.datanoiseon{m,1} = nan;
0389     cal.datanoiseoff{m,1} = nan;
0390   end
0391 
0392   % get info for parallactic angle
0393   az  = mean(dc.antenna0.servo.az)*pi/180;
0394   el  = mean(dc.antenna0.servo.el)*pi/180;
0395   ra  = mean(dc.antenna0.tracker.equat_geoc(:,1))*15*pi/180;
0396   dec = mean(dc.antenna0.tracker.equat_geoc(:,2))*pi/180;
0397   lat = mean(dc.antenna0.tracker.siteActual(:,2))*pi/180;
0398   lst = mean(dc.antenna0.tracker.lst)*15*pi/180;
0399   ha  = (lst - ra);
0400   cosel_cospa = sin(lat)*cos(dec) - sin(dec)*cos(lat)*cos(ha);
0401   cosel_sinpa = sin(ha)*cos(lat);
0402   pa  = atan2(cosel_sinpa, cosel_cospa);
0403   cal.pa(m,1) = pa;
0404   cal.el(m,1) = el;
0405   cal.az(m,1) = az;
0406   
0407   % get info on the temperatures
0408   cal.Tamb(m,1) = mean(dc.array.weather.airTemperature);
0409 end
0410 % make sure that our off-source values are correction
0411 cal.onSource = ~cal.isoff & cal.onSource;
0412 
0413 % last we need to check that if PnoiseOn is flagged, so should PnoiseOff
0414 flag = isnan(cal.PnoiseOn) | isnan(cal.PnoiseOff);
0415 cal.PnoiseOn(flag) = nan;
0416 cal.PnoiseOff(flag) = nan;
0417 cal.varNoiseOn(flag) = nan;
0418 cal.varNoiseOff(flag) = nan;
0419 
0420 % next we get rid of any observations that are not on a recognized
0421 % calibrator.
0422 ind = cal.whatSource == 99;
0423 cal = structcut(cal, ~ind);
0424 
0425 return;
0426 
0427 function cal2 = calculateYFactors(cal, ptout)
0428 
0429 % calcualte the y-factors
0430 
0431 % ok.  now we start with the first one, find the corresponding match, and
0432 % calculate our y-factors.
0433 % y-factor for our measurements:
0434 % ( PblankND - PblankOFF) ./ (Psrc - PblanckOff)
0435 
0436 % let us try this using the pointing offsets too!!
0437 index = 1;
0438 for m=1:(length(cal.time)/2)
0439 
0440   % first we find the correspondance
0441   ind =  cal.whatSource==cal.whatSource(2*m-1) & ...   % match the source
0442       cal.onSource == ~cal.onSource(2*m-1);          % if it is on, match an off
0443   f = find(ind);
0444   % search for the closest one in time.
0445   ff = find( abs(cal.time(f) - cal.time(2*m-1)) == min(abs(cal.time(f) - ...
0446       cal.time(2*m-1))));
0447   matchVal = f(ff);
0448 
0449   cal2.source(index,1)     = cal.whatSource(2*m-1);
0450   cal2.time(index,1)       = cal.time(2*m-1);
0451   cal2.sourceFlux(index,1) = cal.sourceFlux(2*m-1);
0452   cal2.fluxErr(index,1)    = cal.fluxErr(2*m-1);
0453   cal2.theorTau(index,1)   = cal.tau(2*m-1);
0454   
0455   % ok...so we have our match, [2*m-1, matchVal]
0456   if(cal.onSource(2*m-1))
0457     Psrc_noiseon    = cal.PnoiseOn(2*m-1,:);
0458     Psrc_noiseoff   = cal.PnoiseOff(2*m-1,:);
0459     Pblank_noiseon  = cal.PnoiseOn(matchVal,:);
0460     Pblank_noiseoff = cal.PnoiseOff(matchVal,:);
0461     Vsrc_noiseon    = cal.varNoiseOn(2*m-1,:);
0462     Vsrc_noiseoff   = cal.varNoiseOff(2*m-1,:);
0463     Vblank_noiseon  = cal.varNoiseOn(matchVal,:);
0464     Vblank_noiseoff = cal.varNoiseOff(matchVal,:);
0465   else
0466     Psrc_noiseon    = cal.PnoiseOn(matchVal,:);
0467     Psrc_noiseoff   = cal.PnoiseOff(matchVal,:);
0468     Pblank_noiseon  = cal.PnoiseOn(2*m-1,:);
0469     Pblank_noiseoff = cal.PnoiseOff(2*m-1,:);
0470     Vsrc_noiseon    = cal.varNoiseOn(matchVal,:);
0471     Vsrc_noiseoff   = cal.varNoiseOff(matchVal,:);
0472     Vblank_noiseon  = cal.varNoiseOn(2*m-1,:);
0473     Vblank_noiseoff = cal.varNoiseOff(2*m-1,:);
0474   end    
0475   Pon_Poff   = Pblank_noiseon - Pblank_noiseoff;
0476   Psrc_Pon   = Psrc_noiseoff - Pblank_noiseon;
0477   Psrc_Poff  = Psrc_noiseoff - Pblank_noiseoff;
0478   Psrc_onoff = Psrc_noiseon  - Psrc_noiseoff;
0479 
0480   cal2.blank_noiseon(index,:)  = Pblank_noiseon;
0481   cal2.blank_noiseoff(index,:) = Pblank_noiseoff;
0482   cal2.src_noiseon(index,:)    = Psrc_noiseon;
0483   cal2.src_noiseoff(index,:)   = Psrc_noiseoff;
0484   cal2.var_blank_noiseon(index,:)  = Vblank_noiseon;
0485   cal2.var_blank_noiseoff(index,:) = Vblank_noiseoff;
0486   cal2.var_src_noiseon(index,:)    = Vsrc_noiseon;
0487   cal2.var_src_noiseoff(index,:)   = Vsrc_noiseoff;
0488 
0489   % let us find the second measure of Psrc - Psky using the pointing crosses.
0490   ft = find(abs(ptout.time - cal2.time(index,1)) == min(abs(ptout.time - ...
0491       cal2.time(index,1))));
0492   Pdiff2 = [ptout.peakheight(ft,:) 0 0];
0493   
0494   nonlin_ratio = Psrc_onoff./Pon_Poff;
0495   
0496   ynoisesky   = Pblank_noiseoff./Pblank_noiseon;
0497   
0498   ycrazysky = (Pblank_noiseon - Pblank_noiseoff)./(Psrc_noiseoff - ...
0499       Pblank_noiseoff);
0500 
0501   ycrazysky2  = (Pblank_noiseon - Pblank_noiseoff)./(Pdiff2);
0502   
0503   cal2.ynoise(index,:) = ynoisesky;
0504   cal2.ycrazy(index,:) = ycrazysky;
0505   cal2.ycrazy2(index,:)= ycrazysky2;
0506   cal2.elev(index,:)   = mean(cal.el([2*m-1 matchVal]));
0507   cal2.pa(index,:)     = mean(cal.pa([2*m-1 matchVal]));
0508   cal2.az(index,:)     = mean(cal.az([2*m-1 matchVal]));
0509   cal2.nonlin(index,:) = nonlin_ratio;
0510   
0511   index = index+1;
0512 end
0513 
0514 % let us flag any data that has naned values
0515 indbad = isnan(cal2.src_noiseon) | isnan(cal2.src_noiseoff) | ...
0516     isnan(cal2.blank_noiseon) | isnan(cal2.blank_noiseoff);
0517 indbad = sum(indbad,2) > 0;
0518 
0519 if(length(find(indbad))>0)
0520   display(sprintf('Had to cut out %d measurements', length(find(indbad))));
0521   cal2 = structcut(cal2, ~indbad);
0522 end
0523 
0524 return;
0525 
0526 
0527 
0528 function tauVals = calcTauFromSrcDips(cal, d)
0529 
0530 % Uses Method 2 (notebook Pages 94-95) to calculate tau
0531 % section 2.1 of abscal_v2.pdf
0532 
0533 
0534 % ok.  we select the sky dips, and then look at
0535 % ( Pskdip - PblankOFF) ./ (Psrc - PblanckOff)
0536 Tcmb = 2.728;
0537 
0538 [si ei] = findStartStop(d.index.skydip.slow);
0539 
0540 for m=1:length(si)
0541   ind = zeros(size(d.array.frame.features));
0542   ind(si(m):ei(m)) = 1;
0543   ind = logical(ind);
0544   
0545   dc = framecut(d, ind);
0546   time(m,1)    = mean(dc.array.frame.utc);
0547   Tamb(m,1)    = mean(dc.array.weather.airTemperature + 273.15);
0548   TambErr(m,1) = std(dc.array.weather.airTemperature + 273.15);
0549 
0550   % next to calculate the slope, let us just get the data where we are
0551   % actually scanning in elevation % scan down at 0.5deg/s
0552   del = deriv(dc.antenna0.servo.el);
0553   ind = del < -0.4/100;
0554   elvals = dc.antenna0.servo.el(ind);
0555   Ivals  = dc.antenna0.receiver.data(ind, [1 8]);
0556   Iflags = dc.flags.fast(ind, [1 3]);
0557   Ivals(Iflags>0) = nan;
0558 
0559   % Now we can find the corresponding source values
0560   aa = sourceCorrespondance(unique(dc.antenna0.tracker.source));
0561   f  = find( cal.source == aa);
0562   if(isempty(f))
0563     % means it is a skydip at end of previous sched or beginning of next
0564     % look for closest in time
0565     ff = find( abs(cal.time - time(m)) == min(abs(cal.time - time(m))));
0566     time(m,1) = cal.time(ff);
0567     isCoupled = 0;
0568   else
0569     % of those two, find the closest match in time
0570     ff = find( abs(cal.time(f) - time(m)) == min(abs(cal.time(f) - time(m))));
0571     time(m,1) = cal.time(f(ff));
0572     isCoupled = 1;
0573   end
0574   
0575     
0576   % set up the data we need
0577   Ivals(:,1)  = (Ivals(:,1) - cal.blank_noiseoff(ff,1))./(cal.src_noiseoff(ff,1)  - ...
0578       cal.blank_noiseoff(ff,1));
0579   Ivals(:,2)  = (Ivals(:,2) - cal.blank_noiseoff(ff,2))./(cal.src_noiseoff(ff,2)  - ...
0580       cal.blank_noiseoff(ff,2));
0581   x1          = 1./sin(cal.elev(ff));
0582   Tsrc        = cal.Tsrc(ff);
0583   
0584   % next we fit for a slope of the line
0585   mx = []; 
0586   b  = [];
0587   thisflag = 0;
0588   for mm=1:2
0589     elfits = 1./sind(elvals(Iflags(:,mm)==0));
0590     Ifit   = Ivals(:,mm);
0591     Ifit(isnan(Ifit)) = [];
0592     [mx b] = linfit(elfits, Ifit);
0593     
0594     % and now, we calculate Tau
0595 %    tau(m,mm) = (Tcmb - Tamb(m,1))./(2*x1*Tsrc) + sqrt( ((Tamb(m,1)-Tcmb)./Tsrc).^2 + ...
0596 %    4*mx*x1)./(2*x1);
0597     tau(m,mm) = -1./(2*x1) + sqrt( 1 + ...
0598     4*x1*mx*Tsrc/(Tamb(m,1)-Tcmb))./(2*x1);
0599     
0600     
0601     thisflag = 0;
0602     % flag if the x-range is not that large.
0603     if( (max(elfits) - min(elfits)) < 0.12 )
0604       thisflag = 1;
0605     end
0606     if (min(elvals) > 60)
0607       thisflag = 1;
0608     end
0609     
0610     % flag if the fit residuals are bad too
0611     fitline = mx.*elfits + b;
0612     resid   = Ifit - fitline;
0613     residrat= sqrt(nanvar(resid))*10;
0614     if(isCoupled)
0615       % when the skydip is associated with the source, the residuals should be
0616       % closer to zero.
0617       if((residrat > 0.5))
0618     thisflag = 1;
0619       end
0620     else
0621       if(any(residrat > 2))
0622     thisflag = 1;
0623       end
0624     end
0625     flag(m,mm) = thisflag;
0626   end
0627   % need to calculate the errors as well
0628 end
0629 tauErr = zeros(size(tau));
0630 flag = flag | isnan(tau);
0631 tauVals = [time tau tauErr flag];
0632 
0633 
0634 
0635 return;
0636 
0637 
0638 
0639 function cal2 = calcSourceTemp(cal2, pointoff)
0640 
0641 % calculates the source temperature given a pointing offset
0642 
0643 % constants
0644 apEff = 0.55;  % 55% theoretical
0645 k_b   = 1.38e-23;
0646 
0647 % first we need to find a pointing offset that matches each our data set
0648 for m=1:length(cal2.source)
0649   f = find(pointoff(:,2) == cal2.source(m));
0650   timediff = (pointoff(f,1) - cal2.time(m))*24*60;
0651   
0652   ff = find(abs(timediff) == min(abs(timediff)));
0653   if(isempty(ff))
0654     display('calcTnoise_v2::calcSourceTemp:: No good pointing cross');
0655     thisCorrection = nan;
0656   elseif(abs(timediff(ff))  > 20)
0657     display('calcTnoise_v2::calcSourceTemp:: Pointing Cross not commensurate with cal observation');
0658     thisCorrection = nan;
0659   else
0660     thisCorrection = pointoff(f(ff), 5);
0661   end
0662   
0663   Tsrc(m,1) = cal2.sourceFlux(m).*(pi*(6.1/2).^2*apEff)./(2*k_b*1e26)*thisCorrection;
0664 end
0665 
0666 cal2.Tsrc = Tsrc;
0667 cal2.Tsrc_o = cal2.sourceFlux.*(pi*(6.1/2).^2*apEff)./(2*k_b*1e26);
0668 
0669 return;
0670 
0671 function cal2 = matchOpacity(cal2, opacityInfo)
0672 
0673 % flag the bad values
0674 f = find(opacityInfo(:,6));
0675 opacityInfo(f,:) = [];
0676 
0677 if(isempty(opacityInfo))
0678   display('calcTnoise_v2:: matchOpacity:: You have no good skydip values');
0679   cal2.slope = nan(size(cal2.ynoise));
0680   cal2.Tamb  = nan(size(cal2.time));
0681   cal2.TambErr = nan(size(cal2.time));
0682   return;
0683 end
0684 
0685 for m=1:length(cal2.source)
0686   t = (opacityInfo(:,1) - cal2.time(m))*60*24;
0687   f = find(abs(t) == min(abs(t)));
0688   if(abs(t) > 30)
0689     display('calcTnoise_v2:: matchOpacity:: Sky Dip not commensurate with cal observation');
0690     cal2.slope(m,1:2)    = nan;
0691     cal2.Tamb(m,1)       = nan;
0692     cal2.TambErr(m,1)    = nan;
0693   else
0694     cal2.slope(m,1:2)    = opacityInfo(f, [2 3]);
0695     cal2.Tamb(m,1)       = opacityInfo(f, 4);
0696     cal2.TambErr(m,1)    = opacityInfo(f, 5);
0697   end
0698 end
0699 
0700 return;
0701