Home > reduc > calcTsysPolOnly.m

calcTsysPolOnly

PURPOSE ^

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

SYNOPSIS ^

function tsys = calcTsys(d)

DESCRIPTION ^

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

  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 teh system temperature.

  sjcm

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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 teh 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('calcTsysPolOnly:: No Calibrator Source Observations in this schedule');
0022   display('calcTsysPolOnly:: Can not calculate system temperatures');
0023   tsys = [];
0024   return;
0025 end
0026 
0027 
0028 % now we go through the noise diode events, calculate y-factors, and
0029 % determine whether it corresponds to the offset of a source.
0030 numNoiseObs = size(d.correction.alpha.indices,1);
0031 cal.whatSource = zeros(numNoiseObs,1);
0032 cal.onSource   = zeros(numNoiseObs,1);
0033 cal.time       = zeros(numNoiseObs,1);
0034 cal.tau        = zeros(numNoiseObs,1);
0035 cal.sourceFlux = zeros(numNoiseObs,1);
0036 indices = d.correction.alpha.indices;
0037 for m=1:numNoiseObs
0038   ind = zeros(size(d.antenna0.receiver.utc));
0039   ind(indices(m,1):indices(m,3)) = 1;
0040   ind = logical(ind);
0041   dc = framecut(d, ind);
0042   % we should have 4 intensity sky values, states on/off for chans 1/2
0043   switch(size(d.antenna0.receiver.data,2))
0044     case 10
0045       dataNoiseOn = d.antenna0.receiver.data((indices(m,4)+5):indices(m,2), [1 9]);
0046       dataNoiseOff= d.antenna0.receiver.data([indices(m,1):(indices(m,4)-5) ...
0047         (indices(m,2)+5):indices(m,3)], [1 9]);
0048       
0049     case 8
0050       dataNoiseOn = d.antenna0.receiver.data((indices(m,4)+5):indices(m,2), [1 8]);
0051       dataNoiseOff= d.antenna0.receiver.data([indices(m,1):(indices(m,4)-5) ...
0052         (indices(m,2)+5):indices(m,3)], [1 8]);
0053   end
0054   
0055   % calculate the opacities from theory
0056   % WILL MODIFY LATER TO GET THE OPACITY FROM SKYDIPS
0057   tau = calcOpacity(5, dc.array.weather.airTemperature+273.15, ...
0058       dc.array.weather.relativeHumidity, ...
0059       length(dc.array.weather.relativeHumidity), 1);
0060   cal.tau(m,1) = mean(tau)./cos( mean(dc.antenna0.servo.el)*pi/180);
0061   
0062   % make sure this obs is on a calibrator
0063   thisSource = unique(dc.antenna0.tracker.source);
0064   if(length(thisSource)>1)
0065     cal.whatSource(m,1) = 99;
0066   else
0067     [cal.whatSource(m,1) cal.sourceFlux(m,1)] = sourceCorrespondance(thisSource);
0068   end
0069   
0070   if(cal.whatSource(m) ~= 99)
0071     % we should now determine whether it is on or off the source.
0072     offsets = unique(dc.antenna0.tracker.horiz_off(:,1));
0073     if(offsets==0)
0074       cal.onSource(m,1) = 1;
0075     end
0076     cal.time(m,1) = dc.antenna0.receiver.utc(1);
0077     cal.PnoiseOn(m,:)  = mean(dataNoiseOn);
0078     cal.PnoiseOff(m,:) = mean(dataNoiseOff);
0079     cal.varNoiseOn(m,:) = var(dataNoiseOn);
0080     cal.varNoiseOff(m,:) = var(dataNoiseOff);
0081   else
0082     cal.PnoiseOn(m,1:2) = nan;
0083     cal.PnoiseOff(m,1:2) = nan;
0084     cal.varNoiseOn(m,1:2) = nan;
0085     cal.varNoiseOff(m,1:2) = nan;
0086   end
0087   
0088 end
0089 
0090 % next we get rid of any observations that are not on a recognized
0091 % calibrator.
0092 ind = cal.whatSource == 99;
0093 
0094 cal = structcut(cal, ~ind);
0095 
0096 
0097 % check if we actually have any source
0098 % 18-Mar-2011 (MAS) -- Modified, this crashed if structct() rendered cal no
0099 % longer a structure (which is possible!).
0100 if(~isstruct(cal) || length(cal.whatSource)<1)
0101   display('calcTsysPolOnly:: No Calibrator Source Observations in this schedule');
0102   display('calcTsysPolOnly:: Can not calculate system temperatures');
0103   tsys = [];
0104   return;
0105 end
0106 
0107 % bandwidths.  have to change when we bump up the gain on LNA3 again.
0108 if(d.array.frame.utc > date2mjd(2011, 06, 01))
0109   bw = [763 696]*1e6;
0110 else
0111   bw = [1000 1000]*1e6;
0112 end
0113 
0114 
0115 % ok.  now we start with the first one, find the corresponding match, and
0116 % calculate our y-factors.
0117 index = 1;
0118 for m=1:(length(cal.time)/2)
0119   
0120   % first we find the correspondance
0121   ind =  cal.whatSource==cal.whatSource(2*m-1) & ...   % match the source
0122       cal.onSource == ~cal.onSource(2*m-1);          % if it is on, match an off
0123   f = find(ind);
0124   % search for the closest one in time.
0125   ff = find( abs(cal.time(f) - cal.time(2*m-1)) == min(abs(cal.time(f) - ...
0126       cal.time(2*m-1))));
0127   matchVal = f(ff);
0128 
0129   % ok...so we have our match, [2*m-1, matchVal]
0130   if(cal.onSource(2*m-1))
0131     Psrc_noiseon    = cal.PnoiseOn(2*m-1,:);
0132     Psrc_noiseoff   = cal.PnoiseOff(2*m-1,:);
0133     Pblank_noiseon  = cal.PnoiseOn(matchVal,:);
0134     Pblank_noiseoff = cal.PnoiseOff(matchVal,:);
0135     Vsrc_noiseon    = cal.varNoiseOn(2*m-1,:);
0136     Vsrc_noiseoff   = cal.varNoiseOff(2*m-1,:);
0137     Vblank_noiseon  = cal.varNoiseOn(matchVal,:);
0138     Vblank_noiseoff = cal.varNoiseOff(matchVal,:);
0139   else
0140     Psrc_noiseon    = cal.PnoiseOn(matchVal,:);
0141     Psrc_noiseoff   = cal.PnoiseOff(matchVal,:);
0142     Pblank_noiseon  = cal.PnoiseOn(2*m-1,:);
0143     Pblank_noiseoff = cal.PnoiseOff(2*m-1,:);
0144     Vsrc_noiseon    = cal.varNoiseOn(matchVal,:);
0145     Vsrc_noiseoff   = cal.varNoiseOff(matchVal,:);
0146     Vblank_noiseon  = cal.varNoiseOn(2*m-1,:);
0147     Vblank_noiseoff = cal.varNoiseOff(2*m-1,:);
0148   end    
0149   
0150   % set up our y-values
0151   %  there are many different ways, including when on/off source.  This will
0152   %  tell us about compression.
0153   ynoisesrc   = Psrc_noiseoff./Psrc_noiseon;
0154   ynoisesky   = Pblank_noiseoff./Pblank_noiseon;
0155 
0156   ycrazysrc = (Psrc_noiseon - Psrc_noiseoff)./(Psrc_noiseoff - ...
0157       Pblank_noiseoff);
0158   ycrazysky = (Pblank_noiseon - Pblank_noiseoff)./(Psrc_noiseoff - ...
0159       Pblank_noiseoff);
0160   
0161   % NOISE PROPAGATION TO COME LATER.
0162   %varynoise = Vsrc_noiseoff./(Psrc_noiseoff.^2) + Vsrc_noiseon./Psrc_noiseon.^2;
0163   %varynoise = varynoise.*ynoise.^2;
0164   %varycrazy = (Vblank_noiseon + Vblank_noiseoff)./(Pblank_noiseon - Pblank_noiseoff).^2 + (Vsrc_noiseon + ...
0165   %    Vblank_noiseoff)./(Psrc_noiseon - Pblank_noiseoff).^2;
0166   %varycrazy = varycrazy.*ycrazy.^2;
0167   
0168   % calculate the Tsrc
0169   %apEff = getApEff(d.array.frame.utc(1));
0170   apEff = 0.55;  % 55% theoretical
0171   k_b   = 1.38e-23;
0172   Tsrc  = cal.sourceFlux(2*m-1).*(pi*(6.1/2).^2*apEff)./(2*k_b*1e26);
0173 
0174   % get the mean tau
0175   tauVal = mean(cal.tau([2*m-1 matchVal]));
0176   
0177   % and now we put it all together
0178 
0179   tsys.tnoise.src(index,:) = ycrazysrc.*Tsrc.*exp(-tauVal);
0180   tsys.tnoise.sky(index,:)    = ycrazysky.*Tsrc.*exp(-tauVal);
0181   
0182   Tcmb = 2.728*exp(-tauVal);
0183   Tatm = mean(dc.array.weather.airTemperature + 273.15).*(1-exp(-tauVal));
0184   Tload = mean(dc.antenna0.thermal.coldLoad);
0185   
0186   % cannot actually calculate the system temperature from this data anymore.
0187   % but we can calculate Tground.
0188   tsys.tground.src(index,:) = Tload - Tcmb - Tatm + ...
0189       tsys.tnoise.src(index,:)./(ynoisesrc-1);
0190   tsys.tground.sky(index,:) = Tload - Tcmb - Tatm + ...
0191       tsys.tnoise.sky(index,:)./(ynoisesky-1);
0192   
0193   tsys.time(index,1) = cal.time(2*m-1);
0194   tsys.source(index,1) = cal.whatSource(2*m-1);
0195   tsys.Tsrc(index,:) = Tsrc;
0196   tsys.tau(index,:) = tauVal;
0197   
0198   tsys.ynoise.src(index,:) = ynoisesrc;
0199   tsys.ycrazy.src(index,:) = ycrazysrc;
0200   tsys.ynoise.sky(index,:) = ynoisesky;
0201   tsys.ycrazy.sky(index,:) = ycrazysky;
0202   
0203   % and now, the gains.
0204   g1src = (Psrc_noiseon - ...
0205       Psrc_noiseoff)./(k_b*tsys.tnoise.src(index,:).*bw);
0206   g1sky = (Pblank_noiseon - ...
0207       Pblank_noiseoff)./(k_b*tsys.tnoise.sky(index,:).*bw); 
0208 
0209   g1srcdB = 10*log10(g1src);
0210   g1skydB = 10*log10(g1sky);
0211 
0212   thisGain = [cal.time(2*m-1) g1src cal.time(matchVal) g1sky];
0213   thisGaindB = [cal.time(2*m-1) g1srcdB cal.time(matchVal) g1skydB];
0214   
0215   tsys.gain(index,:) = thisGain;
0216   tsys.gaindB(index,:) = thisGaindB;
0217   
0218   
0219   % and now, errors
0220   %ddycrazy = Tsrc.*exp(tauVal).*ynoise./(1-ynoise);
0221   %ddy      = ycrazy.*Tsrc*exp(tauVal)./(1-ynoise).^2;
0222   %varTot   = ddycrazy.^2.*varycrazy + ddy.^2.*varynoise;
0223   %tsys.err(index,:) = sqrt(varTot);
0224 
0225   index = index+1;
0226 end
0227 
0228 
0229 % for confusing other functions
0230 tsys.val = nan(size(tsys.tnoise.sky));
0231 
0232 
0233 return;
0234 
0235   
0236
Generated on Sun 14-Jun-2015 17:12:45 by m2html © 2005