calcTnoise

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

  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:: No Calibrator Source Observations in this schedule');
0030   display('calcTnoise:: 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:: No pointing corrections found for your sources');
0040   display('calcTnoise:: 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:: No skydips found in your data');
0050   display('calcTnoise:: 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:: No Calibrator Source Observations in this schedule');
0062   display('calcTnoise:: 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);
0069 % pointing offset
0070 cal2 = calcSourceTemp(cal2, pointoff2);
0071 
0072 % next, we match up the opacity info with each source as well.
0073 cal2 = matchOpacity(cal2, opacityInfo);
0074 
0075 if(~all(isnan(cal2.slope)))
0076   % next we can calculate tau from the skydips and srcs
0077   tau  =  calcTauFromSrcDips(cal2, d);
0078 end
0079 
0080   
0081 if(all(isnan(cal2.Tamb)))
0082   display('calcTnoise:: No unflagged Opacity values in your data');
0083   display('calcTnoise:: Can not calculate noise diode');
0084   tnoise = [];
0085   return;
0086 end
0087 
0088 % and now, we are equipped with slope, Tamb, Tcmb, ycrazy, elevation, Tsrc
0089 % (all in cal2), and we can solve for opacity and TND
0090 % solver sucks!  we can do analytic solution.
0091 %  1.  tauo = log(y*Tsrc/Tnd)*sin(el_src)
0092 %  2.  tauo = slope*Tnd/(Tamb-Tcmb)
0093 for m=1:length(cal2.source)
0094   for mm=1:2
0095     thiscal.ycrazy = cal2.ycrazy(m,mm);
0096     thiscal.slope  = cal2.slope(m,mm);
0097     thiscal.Tamb   = cal2.Tamb(m);
0098     thiscal.Tcmb   = Tcmb;
0099     thiscal.elev   = cal2.elev(m);
0100     thiscal.Tsrc   = cal2.Tsrc(m);
0101     x0 = [0.01 3];
0102 
0103     alpha = 2*thiscal.slope./(thiscal.Tamb - thiscal.Tcmb); 
0104     beta  = thiscal.ycrazy.*thiscal.Tsrc.*sin(thiscal.elev);
0105     gamma = sin(thiscal.elev);
0106     
0107     tnd  =  (-gamma + sqrt(gamma.^2+4*alpha*beta))./(2*alpha);
0108     tauo1  = tnd.*alpha;
0109     tauo2  = beta./tnd - gamma;
0110     tau01(m,mm) = tauo1;
0111     Tnd(m,mm)  = tnd;
0112   end
0113 end
0114 cal2.calcTau = tau01;
0115 cal2.Tnd     = Tnd;
0116 cal2.point   = pointoff2;
0117 cal2.opacity = opacityInfo;
0118 cal2.tauSrcDip = tau;
0119 
0120 
0121 % need to calculate the errors
0122 cal2.TndErr  = zeros(size(cal2.Tnd));
0123 cal2.calcTauErr = zeros(size(cal2.calcTau));
0124 
0125 tnoise = cal2;
0126 
0127 return;
0128 
0129 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0130 %
0131 %  SUB-FUNCTIONS
0132 %
0133 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0134 
0135   
0136 
0137 function [pointoff1 pointoff2] = calculatePointingCorrection(d)
0138 
0139 % just takes our data and puts it in the format we want for the calculations
0140 
0141 % FWHM = 2sqrt(2ln(2)) call it beta
0142 beta = 0.73/(2*sqrt(2*log(2)));
0143 
0144 % old way:
0145 point = d.correction.pointing;
0146 pointoff1 = [];
0147 while(~isempty(point))
0148   dt = (point(:,1) - point(1,1))*24*60;
0149   ind = find(dt<6 & point(:,2) == point(1,2));
0150   angoff = spaceangle(point(ind,3), point(ind,4), ...
0151       point(ind,3)+point(ind,5), point(ind,4)+point(ind,6), 'deg');
0152   ptmean = nanmean(point(ind,:),1);
0153   pointoff1 = [pointoff1; [ptmean([1 2]) nanmean(angoff) nanstd(angoff)]];
0154   point(ind,:) = [];
0155 end
0156 correction = gauss([1 0 beta], pointoff1(:,3));
0157 pointoff1(:,5) = correction;
0158 
0159 % new way
0160 xoff      = d.correction.allpoint(:, [7 8]);
0161 yoff      = d.correction.allpoint(:, [9 10]);
0162 badpoint  = logical(d.correction.allpoint(:, [5 6]));
0163 time      = d.correction.allpoint(:, 1);
0164 source    = d.correction.allpoint(:, 2);
0165 errx      = d.correction.allpoint(:, [11 12]);
0166 erry      = d.correction.allpoint(:, [13 14]);
0167 %
0168 xoff(badpoint) = nan;
0169 yoff(badpoint) = nan;
0170 errx(badpoint) = inf;
0171 erry(badpoint) = inf;
0172 [xoff vxoff]   = vwstat(xoff, errx.^2, 2);
0173 [yoff vyoff]   = vwstat(yoff, erry.^2, 2);
0174 totaloff       = pyth(xoff, yoff);
0175 totalerroff    = pyth(sqrt(vxoff), sqrt(vyoff));
0176 % next, take a mean of values for the same source, within a short amount
0177 % of time.
0178 pointoff2 = [];
0179 while(~isempty(time))
0180   dt = (time(:,1) - time(1,1))*24*60;
0181   ind = find(dt<6 & source(:) == source(1));
0182   [ptmean pterr] =  vwstat(totaloff(ind), totalerroff(ind));
0183   pointoff2 = [pointoff2; [mean(time(ind)) mean(source(ind)) ptmean ...
0184       pterr]];
0185   time(ind) = [];
0186   source(ind) = [];
0187   totaloff(ind) = [];
0188   totalerroff(ind) = [];
0189 end
0190 correction = gauss([1 0 beta], pointoff2(:,3));
0191 pointoff2(:,5) = correction;
0192 
0193 return;
0194 
0195   
0196 
0197 function opacityInfo = calculateOpacityInfo(d)
0198 
0199 % calculates what we need from the opacity, namely the slope and Tamb
0200 
0201 [si ei] = findStartStop(d.index.skydip.slow);
0202 
0203 for m=1:length(si)
0204   ind = zeros(size(d.array.frame.features));
0205   ind(si(m):ei(m)) = 1;
0206   ind = logical(ind);
0207   
0208   dc = framecut(d, ind);
0209   time(m,1)    = mean(dc.array.frame.utc);
0210   Tamb(m,1)    = mean(dc.array.weather.airTemperature + 273.15);
0211   TambErr(m,1) = std(dc.array.weather.airTemperature + 273.15);
0212 
0213   % next to calculate the slope, let's just get the data where we are
0214   % actually scanning in elevation % scan down at 0.5deg/s
0215   del = deriv(dc.antenna0.servo.el);
0216   ind = del < -0.4/100;
0217   elvals = dc.antenna0.servo.el(ind);
0218   Ivals  = dc.antenna0.receiver.data(ind, [1 8]);
0219 
0220 
0221   thisflag = [ 0 0 ];
0222   indflag = dc.flags.fast(ind,[1 3]);
0223   for mm=1:2
0224     thisInd = ~indflag(:,mm);
0225     elfits = 1./sind(elvals(thisInd));
0226     Ifit   = Ivals(thisInd,mm);
0227     if(length(Ifit) < 10)
0228       mx = nan;
0229       b = nan;
0230       thisflag(mm) = 1;
0231     else
0232       % fit for slope of the line
0233       [mx b] = linfit(elfits, Ifit);
0234     end
0235     slopes(m,mm) = mx;
0236     if( (max(elfits) - min(elfits)) < 0.12)
0237       thisflag(mm) = 1;
0238     end
0239     if(min(elfits) > 60)
0240       thisflag(mm) = 1;
0241     end
0242     % flag if the fit residuals are bad too
0243     fitline = mx.*elfits + b;
0244     resid   = Ifit - fitline;
0245     if(rms(resid)*100>1)
0246       residrat= abs(rms(resid)./mean(Ivals)*100);
0247       if(residrat>0.5)
0248     thisflag(mm) = 1;
0249       end
0250     end
0251   end
0252   flag(m,:) = thisflag;
0253 end
0254 
0255 opacityInfo = [time slopes Tamb TambErr flag];
0256 
0257 return;
0258 
0259 
0260   
0261 function cal = obtainPowerMeasurements(d)
0262 
0263 % splits our data up, and does the proper averaging
0264 numNoiseObs = size(d.correction.alpha.indices,1);
0265 cal.whatSource = zeros(numNoiseObs,1);
0266 cal.onSource   = zeros(numNoiseObs,1);
0267 cal.time       = zeros(numNoiseObs,1);
0268 cal.tau        = zeros(numNoiseObs,1);
0269 cal.sourceFlux = zeros(numNoiseObs,1);
0270 indices = d.correction.alpha.indices;
0271 cols = [1 8];
0272 
0273 for m=1:numNoiseObs
0274   ind = zeros(size(d.antenna0.receiver.utc));
0275   ind(indices(m,1):indices(m,3)) = 1;
0276   ind = logical(ind);
0277   dc = framecut(d, ind);
0278   % we should have 4 intensity sky values, states on/off for chans 1/2
0279   
0280   % apply the flags
0281   dataNoiseOn = d.antenna0.receiver.data((indices(m,4)+5):indices(m,2), cols);
0282   dataNoiseOff= d.antenna0.receiver.data([indices(m,1):(indices(m,4)-5) ...
0283     (indices(m,2)+5):indices(m,3)], cols);
0284 
0285   flagsOn  = d.flags.fast((indices(m,4)+5):indices(m,2), [1 3]);
0286   flagsOff = d.flags.fast([indices(m,1):(indices(m,4)-5) (indices(m,2)+5):indices(m,3)], [1 3]);
0287 
0288   dataNoiseOn(flagsOn>0) = nan;
0289   dataNoiseOff(flagsOff>0) = nan;
0290 
0291   % calculate the opacities from theory
0292   tau = calcOpacity(5, dc.array.weather.airTemperature+273.15, ...
0293       dc.array.weather.relativeHumidity, ...
0294       length(dc.array.weather.relativeHumidity), 1);
0295   cal.tau(m,1) = mean(tau);
0296   
0297   % make sure this obs is on a calibrator
0298   thisSource = unique(dc.antenna0.tracker.source);
0299   if(length(thisSource)>1)
0300     cal.whatSource(m,1) = 99;
0301     cal.sourceFlux(m,1) = nan;
0302     cal.fluxErr(m,1)    = nan;
0303     cal.isoff(m,1)      = 0;
0304   else
0305     [cal.whatSource(m,1) cal.sourceFlux(m,1) cal.fluxErr(m,1)] = ...
0306     sourceCorrespondance(thisSource, dc.array.frame.utc(1));
0307     if(isempty(strfind(thisSource{1}, 'ff')))
0308       cal.isoff(m,1) = 0;
0309     else
0310       cal.isoff(m,1) = 1;
0311     end
0312   end
0313   
0314   if(cal.whatSource(m) ~= 99)
0315     % we should now determine whether it's on or off the source.
0316     offsets = unique(dc.antenna0.tracker.horiz_off(:,1));
0317     if(offsets==0)
0318       cal.onSource(m,1) = 1;
0319     else
0320       cal.onSource(m,1) = 0;
0321     end
0322     cal.time(m,1) = dc.antenna0.receiver.utc(1);
0323     cal.PnoiseOn(m,:)  = nanmean(dataNoiseOn);
0324     cal.PnoiseOff(m,:) = nanmean(dataNoiseOff);
0325     cal.varNoiseOn(m,:) = nanvar(dataNoiseOn);
0326     cal.varNoiseOff(m,:) = nanvar(dataNoiseOff);
0327     cal.datanoiseon{m,1} = dataNoiseOn;
0328     cal.datanoiseoff{m,1} = dataNoiseOff;
0329   else
0330     cal.onSource(m,1) = 0;
0331     cal.PnoiseOn(m,1:2) = nan;
0332     cal.PnoiseOff(m,1:2) = nan;
0333     cal.varNoiseOn(m,1:2) = nan;
0334     cal.varNoiseOff(m,1:2) = nan;
0335     cal.datanoiseon{m,1} = nan;
0336     cal.datanoiseoff{m,1} = nan;
0337   end
0338 
0339   % get info for parallactic angle
0340   az  = mean(dc.antenna0.servo.az)*pi/180;
0341   el  = mean(dc.antenna0.servo.el)*pi/180;
0342   ra  = mean(dc.antenna0.tracker.equat_geoc(:,1))*15*pi/180;
0343   dec = mean(dc.antenna0.tracker.equat_geoc(:,2))*pi/180;
0344   lat = mean(dc.antenna0.tracker.siteActual(:,2))*pi/180;
0345   lst = mean(dc.antenna0.tracker.lst)*15*pi/180;
0346   ha  = (lst - ra);
0347   cosel_cospa = sin(lat)*cos(dec) - sin(dec)*cos(lat)*cos(ha);
0348   cosel_sinpa = sin(ha)*cos(lat);
0349   pa  = atan2(cosel_sinpa, cosel_cospa);
0350   cal.pa(m,1) = pa;
0351   cal.el(m,1) = el;
0352   cal.az(m,1) = az;
0353   
0354   % get info on the temperatures
0355   cal.Tamb(m,1) = mean(dc.array.weather.airTemperature);
0356 end
0357 % make sure that our off-source values are correction
0358 cal.onSource = ~cal.isoff & cal.onSource;
0359 
0360 % next we get rid of any observations that are not on a recognized
0361 % calibrator.
0362 ind = cal.whatSource == 99 | isnan(mean(cal.PnoiseOn,2)) | isnan(mean(cal.PnoiseOff,2));
0363 cal = structcut(cal, ~ind);
0364 
0365 
0366 return;
0367 
0368 function cal2 = calculateYFactors(cal)
0369 
0370 % calcualte the y-factors
0371 
0372 % ok.  now we start with the first one, find the corresponding match, and
0373 % calculate our y-factors.
0374 % y-factor for our measurements:
0375 % ( PblankND - PblankOFF) ./ (Psrc - PblanckOff)
0376 
0377 
0378 index = 1;
0379 for m=1:(length(cal.time)/2)
0380 
0381   % first we find the correspondance
0382   ind =  cal.whatSource==cal.whatSource(2*m-1) & ...   % match the source
0383       cal.onSource == ~cal.onSource(2*m-1);          % if it's on, match an off
0384   f = find(ind);
0385   % search for the closest one in time.
0386   ff = find( abs(cal.time(f) - cal.time(2*m-1)) == min(abs(cal.time(f) - ...
0387       cal.time(2*m-1))));
0388   matchVal = f(ff);
0389 
0390   cal2.source(index,1)     = cal.whatSource(2*m-1);
0391   cal2.time(index,1)       = cal.time(2*m-1);
0392   cal2.sourceFlux(index,1) = cal.sourceFlux(2*m-1);
0393   cal2.fluxErr(index,1)    = cal.fluxErr(2*m-1);
0394   cal2.theorTau(index,1)   = cal.tau(2*m-1);
0395   
0396   % ok...so we have our match, [2*m-1, matchVal]
0397   if(cal.onSource(2*m-1))
0398     Psrc_noiseon    = cal.PnoiseOn(2*m-1,:);
0399     Psrc_noiseoff   = cal.PnoiseOff(2*m-1,:);
0400     Pblank_noiseon  = cal.PnoiseOn(matchVal,:);
0401     Pblank_noiseoff = cal.PnoiseOff(matchVal,:);
0402     Vsrc_noiseon    = cal.varNoiseOn(2*m-1,:);
0403     Vsrc_noiseoff   = cal.varNoiseOff(2*m-1,:);
0404     Vblank_noiseon  = cal.varNoiseOn(matchVal,:);
0405     Vblank_noiseoff = cal.varNoiseOff(matchVal,:);
0406   else
0407     Psrc_noiseon    = cal.PnoiseOn(matchVal,:);
0408     Psrc_noiseoff   = cal.PnoiseOff(matchVal,:);
0409     Pblank_noiseon  = cal.PnoiseOn(2*m-1,:);
0410     Pblank_noiseoff = cal.PnoiseOff(2*m-1,:);
0411     Vsrc_noiseon    = cal.varNoiseOn(matchVal,:);
0412     Vsrc_noiseoff   = cal.varNoiseOff(matchVal,:);
0413     Vblank_noiseon  = cal.varNoiseOn(2*m-1,:);
0414     Vblank_noiseoff = cal.varNoiseOff(2*m-1,:);
0415   end    
0416   Pon_Poff   = Pblank_noiseon - Pblank_noiseoff;
0417   Psrc_Pon   = Psrc_noiseoff - Pblank_noiseon;
0418   Psrc_Poff  = Psrc_noiseoff - Pblank_noiseoff;
0419   Psrc_onoff = Psrc_noiseon  - Psrc_noiseoff;
0420 
0421   cal2.blank_noiseon(index,:)  = Pblank_noiseon;
0422   cal2.blank_noiseoff(index,:) = Pblank_noiseoff;
0423   cal2.src_noiseon(index,:)    = Psrc_noiseon;
0424   cal2.src_noiseoff(index,:)   = Psrc_noiseoff;
0425   cal2.var_blank_noiseon(index,:)  = Vblank_noiseon;
0426   cal2.var_blank_noiseoff(index,:) = Vblank_noiseoff;
0427   cal2.var_src_noiseon(index,:)    = Vsrc_noiseon;
0428   cal2.var_src_noiseoff(index,:)   = Vsrc_noiseoff;
0429 
0430   nonlin_ratio = Psrc_onoff./Pon_Poff;
0431   
0432   ynoisesky   = Pblank_noiseoff./Pblank_noiseon;
0433   
0434   ycrazysky = (Pblank_noiseon - Pblank_noiseoff)./(Psrc_noiseoff - ...
0435       Pblank_noiseoff);
0436 
0437   cal2.ynoise(index,:) = ynoisesky;
0438   cal2.ycrazy(index,:) = ycrazysky;
0439   cal2.elev(index,:)   = mean(cal.el([2*m-1 matchVal]));
0440   cal2.pa(index,:)     = mean(cal.pa([2*m-1 matchVal]));
0441   cal2.az(index,:)     = mean(cal.az([2*m-1 matchVal]));
0442   cal2.nonlin(index,:) = nonlin_ratio;
0443   index = index+1;
0444 end
0445 
0446 
0447 function tauVals = calcTauFromSrcDips(cal, d)
0448 
0449 % Uses Method 2 (notebook Pages 94-95) to calculate tau
0450 
0451 % ok.  we select the sky dips, and then look at
0452 % ( Pskdip - PblankOFF) ./ (Psrc - PblanckOff)
0453 Tcmb = 2.728;
0454 
0455 [si ei] = findStartStop(d.index.skydip.slow);
0456 
0457 for m=1:length(si)
0458   ind = zeros(size(d.array.frame.features));
0459   ind(si(m):ei(m)) = 1;
0460   ind = logical(ind);
0461   
0462   dc = framecut(d, ind);
0463   time(m,1)    = mean(dc.array.frame.utc);
0464   Tamb(m,1)    = mean(dc.array.weather.airTemperature + 273.15);
0465   TambErr(m,1) = std(dc.array.weather.airTemperature + 273.15);
0466 
0467   % next to calculate the slope, let's just get the data where we are
0468   % actually scanning in elevation % scan down at 0.5deg/s
0469   del = deriv(dc.antenna0.servo.el);
0470   ind = del < -0.4/100;
0471   elvals = dc.antenna0.servo.el(ind);
0472   Ivals  = dc.antenna0.receiver.data(ind, [1 8]);
0473   Iflags = dc.flags.fast(ind, [1 3]);
0474   Ivals(Iflags>0) = nan;
0475 
0476   % Now we can find the corresponding source values
0477   aa = sourceCorrespondance(unique(dc.antenna0.tracker.source));
0478   f  = find( cal.source == aa);
0479   if(isempty(f))
0480     % means it's a skydip at end of previous sched or beginning of next
0481     % look for closest in time
0482     ff = find( abs(cal.time - time(m)) == min(abs(cal.time - time(m))));
0483     time(m,1) = cal.time(ff);
0484     isCoupled = 0;
0485   else
0486     % of those two, find the closest match in time
0487     ff = find( abs(cal.time(f) - time(m)) == min(abs(cal.time(f) - time(m))));
0488     time(m,1) = cal.time(f(ff));
0489     isCoupled = 1;
0490   end
0491   
0492     
0493   % set up the data we need
0494   Ivals(:,1)  = (Ivals(:,1) - cal.blank_noiseoff(ff,1))./(cal.src_noiseoff(ff,1)  - ...
0495       cal.blank_noiseoff(ff,1));
0496   Ivals(:,2)  = (Ivals(:,2) - cal.blank_noiseoff(ff,2))./(cal.src_noiseoff(ff,2)  - ...
0497       cal.blank_noiseoff(ff,2));
0498   x1          = 1./sin(cal.elev(ff));
0499   Tsrc        = cal.Tsrc(ff);
0500   % next we fit for a slope of the line
0501   mx = []; 
0502   b  = [];
0503   thisflag = 0;
0504   for mm=1:2
0505     elfits = 1./sind(elvals(Iflags(:,mm)==0));
0506     Ifit   = Ivals(:,mm);
0507     Ifit(isnan(Ifit)) = [];
0508     if isempty(Ifit)  
0509         tau(m, mm) = NaN;
0510         flag(m,mm) = 1;
0511     else [mx b] = linfit(elfits, Ifit);
0512     % and now, we calculate Tau
0513     tau(m,mm) = (Tcmb - Tamb(m,1))./(2*x1*Tsrc) + sqrt( ((Tamb(m,1)-Tcmb)./Tsrc).^2 + ...
0514     4*mx*x1)./(2*x1);
0515   
0516     thisflag = 0;
0517     % flag if the x-range is not that large.
0518     if( (max(elfits) - min(elfits)) < 0.12 )
0519       thisflag = 1;
0520     end
0521     if (min(elvals) > 60)
0522       thisflag = 1;
0523     end
0524     
0525     % flag if the fit residuals are bad too
0526     fitline = mx.*elfits + b;
0527     resid   = Ifit - fitline;
0528     residrat= sqrt(nanvar(resid))*10;
0529     if(isCoupled)
0530       % when the skydip is associated with the source, the residuals should be
0531       % closer to zero.
0532       if((residrat > 0.5))
0533     thisflag = 1;
0534       end
0535     else
0536       if(any(residrat > 2))
0537     thisflag = 1;
0538       end
0539     end
0540     flag(m,mm) = thisflag;
0541   end
0542   % need to calculate the errors as well
0543   end
0544 end
0545 tauErr = zeros(size(tau));
0546 flag = flag | isnan(tau);
0547 tauVals = [time tau tauErr flag];
0548 
0549 return;
0550 
0551 
0552 
0553 function cal2 = calcSourceTemp(cal2, pointoff)
0554 
0555 % calculates the source temperature given a pointing offset
0556 
0557 % constants
0558 apEff = 0.55;  % 55% theoretical
0559 k_b   = 1.38e-23;
0560 
0561 % first we need to find a pointing offset that matches each our data set
0562 for m=1:length(cal2.source)
0563   f = find(pointoff(:,2) == cal2.source(m));
0564   timediff = (pointoff(f,1) - cal2.time(m))*24*60;
0565   
0566   ff = find(abs(timediff) == min(abs(timediff)));
0567   if(abs(timediff(ff)) > 20)
0568     display('calcTnoise::calcSourceTemp:: Pointing Cross not commensurate with cal observation');
0569     thisCorrection = nan;
0570   else
0571     thisCorrection = pointoff(f(ff), 5);
0572   end
0573   
0574   Tsrc(m,1) = cal2.sourceFlux(m).*(pi*(6.1/2).^2*apEff)./(2*k_b*1e26)*thisCorrection;
0575 end
0576 
0577 cal2.Tsrc = Tsrc;
0578 
0579 return;
0580 
0581 function cal2 = matchOpacity(cal2, opacityInfo)
0582 
0583 % flag the bad values
0584 f = find(opacityInfo(:,6));
0585 opacityInfo(f,:) = [];
0586 
0587 if(isempty(opacityInfo))
0588   display('calcTnoise::matchOpacity:: You have no good skydip values');
0589   cal2.slope = nan(size(cal2.ynoise));
0590   cal2.Tamb  = nan(size(cal2.time));
0591   cal2.TambErr = nan(size(cal2.time));
0592   return;
0593 end
0594 
0595 for m=1:length(cal2.source)
0596   t = (opacityInfo(:,1) - cal2.time(m))*60*24;
0597   f = find(abs(t) == min(abs(t)));
0598   if(abs(t) > 30)
0599     display('calcTnoise::matchOpacity:: Sky Dip not commensurate with cal observation');
0600     cal2.slope(m,1:2)    = nan;
0601     cal2.Tamb(m,1)       = nan;
0602     cal2.TambErr(m,1)    = nan;
0603   else
0604     cal2.slope(m,1:2)    = opacityInfo(f, [2 3]);
0605     cal2.Tamb(m,1)       = opacityInfo(f, 4);
0606     cal2.TambErr(m,1)    = opacityInfo(f, 5);
0607   end
0608 end
0609 
0610 return;
0611