Description of calcTnoise

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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
