calcTsys

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

  function tsysVals = calcTsys(d)

    should calculate tsys according to page 63 of stephen''s lab book
    when the cold load is included.  Only main difference is with regards
    to the system temperature.

  sjcm

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

0001 function tsys = calcTsys(d)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 %  function tsysVals = calcTsys(d)
0006 %
0007 %    should calculate tsys according to page 63 of stephen''s lab book
0008 %    when the cold load is included.  Only main difference is with regards
0009 %    to the system temperature.
0010 %
0011 %  sjcm
0012 %
0013 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0014 
0015 % we need to know if we are observing any sources.
0016 sources = unique(d.antenna0.tracker.source);
0017 for m=1:length(sources)
0018   [sourceNum(m) sourceFlux(m)] = sourceCorrespondance(sources(m));
0019 end
0020 if(all(sourceNum==99))
0021   display('calcTsys:: No Calibrator Source Observations in this schedule');
0022   display('calcTsys:: Can not calculate system temperatures');
0023   tsys = [];
0024   return;
0025 end
0026 
0027 % now we go through the noise diode events, calculate y-factors, and
0028 % determine whether it corresponds to the offset of a source.
0029 numNoiseObs = size(d.correction.alpha.indices,1);
0030 cal.whatSource = zeros(numNoiseObs,1);
0031 cal.onSource   = zeros(numNoiseObs,1);
0032 cal.time       = zeros(numNoiseObs,1);
0033 cal.tau        = zeros(numNoiseObs,1);
0034 cal.sourceFlux = zeros(numNoiseObs,1);
0035 indices = d.correction.alpha.indices;
0036 
0037 cols = [1 8];
0038 
0039 for m=1:numNoiseObs
0040   ind = zeros(size(d.antenna0.receiver.utc));
0041   ind(indices(m,1):indices(m,3)) = 1;
0042   ind = logical(ind);
0043   dc = framecut(d, ind);
0044   % we should have 4 intensity sky values, states on/off for chans 1/2
0045   dataNoiseOn = d.antenna0.receiver.data((indices(m,4)+5):indices(m,2), cols);
0046   dataNoiseOff= d.antenna0.receiver.data([indices(m,1):(indices(m,4)-5) ...
0047     (indices(m,2)+5):indices(m,3)], cols);
0048 
0049   flagsOn  = d.flags.fast((indices(m,4)+5):indices(m,2), [1 3]);
0050   flagsOff = d.flags.fast([indices(m,1):(indices(m,4)-5) ...
0051     (indices(m,2)+5):indices(m,3)], [1 3]);
0052   
0053   dataNoiseOn(flagsOn>0) = nan;
0054   dataNoiseOff(flagsOff>0) = nan;
0055   
0056   
0057   % get Opacity from skydips
0058   tau = getTauVal(dc);
0059   tau.values(tau.flag>0) = nan;
0060   tau.values = nanmean(tau.values,2);
0061   tau = structcut(tau, ~isnan(tau.values));
0062   % find the closest val
0063   if(~isempty(tau.time))
0064     f = find( abs(tau.time - dc.array.frame.utc(1)) == min(abs(tau.time - ...
0065     dc.array.frame.utc(1))));
0066     tau2use = tau.values(f);
0067   else
0068     display('calcTsys:: Using default tau');
0069     tau2use = 0.01;
0070   end
0071   if(isempty(tau2use))
0072     display('calcTsys:: Using default tau');
0073     tau2use = 0.01;
0074   end
0075   cal.tau(m,1) = tau2use;
0076   
0077   % make sure this obs is on a calibrator
0078   thisSource = unique(dc.antenna0.tracker.source);
0079   if(length(thisSource)>1)
0080     cal.whatSource(m,1) = 99;
0081     cal.sourceFlux(m,1) = nan;
0082     cal.fluxErr(m,1)    = nan;
0083     cal.isoff(m,1)      = 0;
0084   else
0085     [cal.whatSource(m,1) cal.sourceFlux(m,1) cal.fluxErr(m,1)] = sourceCorrespondance(thisSource);
0086     if(isempty(strfind(thisSource{1}, 'ff')))
0087       cal.isoff(m,1) = 0;
0088     else
0089       cal.isoff(m,1) = 1;
0090     end
0091   end
0092   
0093   if(cal.whatSource(m) ~= 99)
0094     % we should now determine whether it is on or off the source.
0095     offsets = unique(dc.antenna0.tracker.horiz_off(:,1));
0096     if(offsets==0)
0097       cal.onSource(m,1) = 1;
0098     else
0099       cal.onSource(m,1) = 0;
0100     end
0101     cal.time(m,1) = dc.antenna0.receiver.utc(1);
0102     cal.PnoiseOn(m,:)  = nanmean(dataNoiseOn);
0103     cal.PnoiseOff(m,:) = nanmean(dataNoiseOff);
0104     cal.varNoiseOn(m,:) = nanvar(dataNoiseOn);
0105     cal.varNoiseOff(m,:) = nanvar(dataNoiseOff);
0106     cal.datanoiseon{m,1} = dataNoiseOn;
0107     cal.datanoiseoff{m,1} = dataNoiseOff;
0108   else
0109     cal.onSource(m,1) = 0;
0110     cal.PnoiseOn(m,1:2) = nan;
0111     cal.PnoiseOff(m,1:2) = nan;
0112     cal.varNoiseOn(m,1:2) = nan;
0113     cal.varNoiseOff(m,1:2) = nan;
0114     cal.datanoiseon{m,1} = nan;
0115     cal.datanoiseoff{m,1} = nan;
0116   end
0117 
0118   % get info for parallactic angle
0119   az  = mean(dc.antenna0.servo.az)*pi/180;
0120   el  = mean(dc.antenna0.servo.el)*pi/180;
0121   ra  = mean(dc.antenna0.tracker.equat_geoc(:,1))*15*pi/180;
0122   dec = mean(dc.antenna0.tracker.equat_geoc(:,2))*pi/180;
0123   lat = mean(dc.antenna0.tracker.siteActual(:,2))*pi/180;
0124   lst = mean(dc.antenna0.tracker.lst)*15*pi/180;
0125   ha  = (lst - ra);
0126   cosel_cospa = sin(lat)*cos(dec) - sin(dec)*cos(lat)*cos(ha);
0127   cosel_sinpa = sin(ha)*cos(lat);
0128   pa  = atan2(cosel_sinpa, cosel_cospa);
0129   cal.pa(m,1) = pa;
0130   cal.el(m,1) = el;
0131   cal.az(m,1) = az;
0132   % get info on the temperatures
0133   cal.Tamb(m,1) = mean(dc.array.weather.airTemperature);
0134 end
0135 % make sure that our off-source values are correction
0136 cal.onSource = ~cal.isoff & cal.onSource;
0137 
0138 
0139 % next we get rid of any observations that are not on a recognized
0140 % calibrator.
0141 ind = cal.whatSource == 99;
0142 
0143 cal = structcut(cal, ~ind);
0144 
0145 % check if we actually have any source
0146 % 18-Mar-2011 (MAS) -- Modified, this crashed if structct() rendered cal no
0147 % longer a structure (which is possible!).
0148 if(~isstruct(cal) || length(cal.whatSource)<1)
0149   display('calcTsys:: No Calibrator Source Observations in this schedule');
0150   display('calcTsys:: Can not calculate system temperatures');
0151   tsys = [];
0152   return;
0153 end
0154 
0155 
0156 % so now we should have a structure that goes:  on source - off source - on
0157 % source - off source PER SOURCE.
0158 %  in between the two cycles, we have pointing crosses.
0159 %  so let us calculate the correction for the pointing offset.
0160 % first let us find the associated correction for each source
0161 point = d.correction.pointing;
0162 pointoff = [];
0163 while(~isempty(point))
0164   dt = (point(:,1) - point(1,1))*24*60;
0165   ind = find(dt<6 & point(:,2) == point(1,2));
0166   angoff = spaceangle(point(ind,3), point(ind,4), ...
0167       point(ind,3)+point(ind,5), point(ind,4)+point(ind,6), 'deg');
0168   ptmean = nanmean(point(ind,:),1);
0169   pointoff = [pointoff; [ptmean([1 2 3 4]) nanmean(angoff) nanstd(angoff)]];
0170   point(ind,:) = [];
0171 end
0172 
0173 % next we calculate the rescale value for the bad pointing
0174 % FWHM = 2sqrt(2ln(2)) call it beta
0175 beta = 0.73/(2*sqrt(2*log(2)));
0176 
0177 correction = gauss([1 0 beta], pointoff(:,5));
0178 if(any(correction<0.5))
0179   warning('calcTsys:: ********** WARNING ************');
0180   warning('calcTsys:: The pointing is off by more than half a beam');
0181   warning('calcTsys:: Be suspect of your data');
0182 end
0183 pointoff(:,7) = correction;
0184 
0185 indBad = pointoff(:,6) > 0.012;
0186 pointoff(indBad,:) = [];
0187 
0188 % we do not rescale the power -- that would scale up everything else as
0189 % well.  We need to rescale only the
0190 % next we rescale the 'on' fluxes
0191 for m=1:length(cal.whatSource)
0192   % only correct when it is on source
0193   if(cal.onSource(m) == 1)
0194     isSource = (pointoff(:,2) == cal.whatSource(m));
0195     f = find(isSource);
0196     if(sum(isSource)>1)
0197       % look for the one that was taken in closest time
0198       isTime   = abs(pointoff(f,1) - cal.time(m));
0199       isTime   = isTime == min(isTime);
0200       ff = find(isTime);
0201       f = f(ff);
0202     end
0203     if(length(f)~=1)
0204       cal.hasPoint(m,1) = 0;
0205     else
0206       cal.hasPoint(m,1)      = 1;
0207       cal.pointscale(m,1)    = pointoff(f,7);
0208       cal.fullpoint(m,:)   = pointoff(f,:);
0209     end
0210   else
0211     cal.hasPoint(m,1) = 0;
0212     cal.pointscale(m,1) = 0;
0213     cal.fullpoint(m,1:7) = 0;
0214   end
0215 end
0216 
0217 % bandwidths.  have to change when we bump up the gain on LNA3 again.
0218 if(d.array.frame.utc > date2mjd(2011, 06, 01))
0219   bw = [763 696]*1e6;
0220 else
0221   bw = [1000 1000]*1e6;
0222 end
0223 
0224 % ok.  now we start with the first one, find the corresponding match, and
0225 % calculate our y-factors.
0226 index = 1;
0227 for m=1:(length(cal.time)/2)
0228 
0229   % only do it if we haev a pointing correction
0230   if(cal.hasPoint(2*m-1))
0231     % first we find the correspondance
0232     ind =  cal.whatSource==cal.whatSource(2*m-1) & ...   % match the source
0233     cal.onSource == ~cal.onSource(2*m-1);          % if it is on, match an off
0234     f = find(ind);
0235     % search for the closest one in time.
0236     ff = find( abs(cal.time(f) - cal.time(2*m-1)) == min(abs(cal.time(f) - ...
0237     cal.time(2*m-1))));
0238     matchVal = f(ff);
0239     
0240     % ok...so we have our match, [2*m-1, matchVal]
0241     if(cal.onSource(2*m-1))
0242       Psrc_noiseon    = cal.PnoiseOn(2*m-1,:);
0243       Psrc_noiseoff   = cal.PnoiseOff(2*m-1,:);
0244       Pblank_noiseon  = cal.PnoiseOn(matchVal,:);
0245       Pblank_noiseoff = cal.PnoiseOff(matchVal,:);
0246       Vsrc_noiseon    = cal.varNoiseOn(2*m-1,:);
0247       Vsrc_noiseoff   = cal.varNoiseOff(2*m-1,:);
0248       Vblank_noiseon  = cal.varNoiseOn(matchVal,:);
0249       Vblank_noiseoff = cal.varNoiseOff(matchVal,:);
0250     else
0251       Psrc_noiseon    = cal.PnoiseOn(matchVal,:);
0252       Psrc_noiseoff   = cal.PnoiseOff(matchVal,:);
0253       Pblank_noiseon  = cal.PnoiseOn(2*m-1,:);
0254       Pblank_noiseoff = cal.PnoiseOff(2*m-1,:);
0255       Vsrc_noiseon    = cal.varNoiseOn(matchVal,:);
0256       Vsrc_noiseoff   = cal.varNoiseOff(matchVal,:);
0257       Vblank_noiseon  = cal.varNoiseOn(2*m-1,:);
0258       Vblank_noiseoff = cal.varNoiseOff(2*m-1,:);
0259     end    
0260     Pon_Poff = Pblank_noiseon - Pblank_noiseoff;
0261     Psrc_Pon = Psrc_noiseoff - Pblank_noiseon;
0262     Psrc_Poff= Psrc_noiseoff - Pblank_noiseoff;
0263     Psrc_onoff= Psrc_noiseon - Psrc_noiseoff;
0264     
0265     nonlin_ratio = Psrc_onoff./Pon_Poff;
0266 
0267     sigmaSource = cal.fluxErr(2*m-1);
0268     
0269     % set up our y-values
0270     % given that our fluxes are rescaled, we cannot calcuate the y-values for
0271     % comparing on-source flux to on-source noise.
0272     ynoisesky   = Pblank_noiseoff./Pblank_noiseon;
0273     
0274     ycrazysky = (Pblank_noiseon - Pblank_noiseoff)./(Psrc_noiseoff - ...
0275     Pblank_noiseoff);
0276     
0277     
0278     % calculate the Tsrc
0279     apEff = 0.55;  % 55% theoretical
0280     k_b   = 1.38e-23;
0281     Tsrc  = cal.sourceFlux(2*m-1).*(pi*(6.1/2).^2*apEff)./(2*k_b*1e26)*(cal.pointscale(2*m-1));
0282     
0283     % get the mean tau
0284     tauVal = mean(cal.tau([2*m-1 matchVal]));
0285     tauVar = var(cal.tau([2*m-1 matchVal]));
0286 
0287     % get the CMB temp and amb temp
0288     Tcmb = 2.728;
0289 
0290     
0291     % check for the elevation differences, if they're not different, then
0292     % just use the previous way.  If they are, then use our new method.
0293     elvals = cal.el([2*m-1 matchVal]);
0294     diffel = abs(diff(elvals)*180/pi);
0295     srcVal = mean(cal.whatSource([2*m-1 matchVal]));
0296     isOff  = (diffel > 0.4) & (srcVal == 1 | srcVal == 2 | srcVal == 8);
0297     
0298     if(~isOff)
0299       tsys.tnoise(index,:)    = ycrazysky.*Tsrc.*exp(-tauVal/cos(mean(elvals)));
0300     else
0301       onSource = cal.onSource([2*m-1 matchVal]);
0302       fon = find(onSource);
0303       foff = find(~onSource);
0304       Tamb =  cal.Tamb([2*m-1 matchVal])+273.15;
0305       Tamb(:) = mean(Tamb);
0306       
0307       
0308       term1 = Tsrc.*exp(-tauVal/cos(elvals(fon)));
0309       term2 = Tcmb.*exp(-tauVal/cos(elvals(fon))) - Tcmb.*exp(-tauVal/cos(elvals(foff)));
0310       term3 = Tamb(fon).*exp(-tauVal/cos(elvals(fon))) - Tamb(foff).*exp(-tauVal/cos(elvals(foff)));
0311 
0312       denom = term1 + term2 - term3;
0313       
0314       tsys.tnoise(index,:) = ycrazysky.*denom;
0315       
0316     end
0317     
0318     
0319     tsys.time(index,1) = cal.time(2*m-1);
0320     tsys.source(index,1) = cal.whatSource(2*m-1);
0321     tsys.Tsrc(index,:) = Tsrc;
0322     tsys.tau(index,:) = tauVal;
0323     
0324     tsys.ynoise(index,:) = ynoisesky;
0325     tsys.ycrazy(index,:) = ycrazysky;
0326     tsys.nonlin(index,:) = nonlin_ratio;
0327 
0328     tsys.src_noiseoff(index,:) = Psrc_noiseoff;
0329     
0330     % and now, the gains.
0331     g1sky = (Pblank_noiseon - ...
0332     Pblank_noiseoff)./(k_b*tsys.tnoise(index,:).*bw); 
0333     
0334     g1skydB = 10*log10(g1sky);
0335     
0336     thisGain = [cal.time(2*m-1) cal.time(matchVal) g1sky];
0337     thisGaindB = [cal.time(2*m-1) cal.time(matchVal) g1skydB];
0338     
0339     tsys.gain(index,:) = thisGain;
0340     tsys.gaindB(index,:) = thisGaindB;
0341     
0342 
0343     tsys.pa(index,1) = cal.pa(2*m-1);
0344     tsys.az(index,1) = cal.fullpoint(2*m-1,3);
0345     tsys.el(index,1) = cal.el(2*m-1);
0346     
0347     % NOW FOR THE ERROR PROPAGATION  page 85-89 in Stephen's CBASS notebook #2
0348     thirdTerm  = tsys.tnoise(index,:).^2.*tauVar;
0349     secondTerm = tsys.tnoise(index,:).^2.*(sigmaSource.^2 + ...
0350     cal.fullpoint(2*m-1, 5).^2./beta.^4.*cal.fullpoint(2*m-1,6).^2);
0351     firstTerm  = tsys.tnoise(index,:).^2./Pon_Poff.^2.*(Vblank_noiseon + ...
0352     (Psrc_Pon./Psrc_Poff).^2.*Vblank_noiseoff + ...
0353     ycrazysky.^2.*Vsrc_noiseoff);
0354     sigmaTnd = firstTerm + secondTerm + thirdTerm;
0355 
0356     tsys.err1(index,:) = firstTerm;
0357     tsys.err2(index,:) = secondTerm;
0358     tsys.err3(index,:) = thirdTerm;
0359     tsys.err(index,:)  = sigmaTnd;
0360     tsys.errTau(index,1) = sqrt(tauVar);
0361     tsys.errSrc(index,1) = sigmaSource;
0362     index = index+1;
0363   else
0364     tsys.tnoise(index,:) = [nan nan];
0365     tsys.tnoise(index,:) = [nan nan];
0366     tsys.time(index,1) = cal.time(2*m-1);
0367     tsys.source(index,1) = cal.whatSource(2*m-1);
0368     tsys.el(index,1)     = cal.el(2*m-1);
0369     tsys.pa(index,1)     = cal.pa(2*m-1);
0370     tsys.az(index,1)     = cal.az(2*m-1);
0371     tsys.Tsrc(index,:) = nan;
0372     tsys.tau(index,:) = nan;
0373     tsys.ynoise(index,:) = [nan nan];
0374     tsys.ycrazy(index,:) = [nan nan];
0375     tsys.nonlin(index,:) = [nan nan];
0376     tsys.src_noiseoff(index,:) = [nan nan];
0377     tsys.gain(index,:) = [nan nan nan nan];
0378     tsys.gaindB(index,:) = [nan nan nan nan];
0379     tsys.err1(index,:) = [nan nan];
0380     tsys.err2(index,:) = [nan nan];
0381     tsys.err3(index,:) = [nan nan];
0382     tsys.err(index,:) = [nan nan];
0383     tsys.errTau(index,1) = nan;
0384     tsys.errSrc(index,1) = nan;
0385     index = index+1;
0386   end
0387 end
0388 
0389 
0390 % re-arrange the outputs
0391 
0392 tnoise.source = tsys.source;
0393 tnoise.time   = tsys.time;
0394 tnoise.ynoise = tsys.ynoise;
0395 tnoise.ycrazy = tsys.ycrazy;
0396 tnoise.nonlin = tsys.nonlin;
0397 tnoise.src_noiseoff = tsys.src_noiseoff;
0398 tnoise.elev   = tsys.el;
0399 tnoise.pa     = tsys.pa;
0400 tnoise.az     = tsys.az;
0401 tnoise.Tsrc   = tsys.Tsrc;
0402 tnoise.Tnd    = tsys.tnoise;
0403 tnoise.TndErr = tsys.err;
0404 
0405 tnoise.calcTau= nan(size(tsys.tnoise));
0406 tnoise.calcTauErr = nan(size(tsys.tnoise));
0407 tnoise.tauSrcDip  = [];
0408 
0409 tsys = tnoise;
0410 
0411 return;
0412 
0413   
0414