time   calls  line
                  1 function tsys = calcTsys(d)
                  2 
                  3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                  4 %
                  5 %  function tsysVals = calcTsys(d)
                  6 %
                  7 %    should calculate tsys according to page 49 of stephen's lab book
                  8 %    (combining equations 27 and 28 in calibration document)
                  9 %
                 10 %  sjcm
                 11 %
                 12 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 13 
                 14 % we know when teh noise source went on and off from 
                 15 %  d.correction.alpha.indices = [blankStart noiseOff endBlank noiseOn]
                 16 
                 17 % we need to know if we're observing any sources.
             1   18 sources = unique(d.antenna0.tracker.source); 
             1   19 for m=1:length(sources) 
             2   20   [sourceNum(m) sourceFlux(m)] = sourceCorrespondance(sources(m)); 
             2   21 end 
             1   22 if(all(sourceNum==99)) 
             1   23   display('No Calibrator Source Observations in this schedule'); 
             1   24   display('Can not calculate system temperatures'); 
             1   25   tsys = []; 
             1   26   return; 
                 27 end
                 28 
                 29 
                 30 % now we go through the noise diode events, calculate y-factors, and
                 31 % determine whether it corresponds to the offset of a source.
                 32 numNoiseObs = size(d.correction.alpha.indices,1);
                 33 cal.whatSource = zeros(numNoiseObs,1);
                 34 cal.onSource   = zeros(numNoiseObs,1);
                 35 cal.time       = zeros(numNoiseObs,1);
                 36 cal.tau        = zeros(numNoiseObs,1);
                 37 cal.sourceFlux = zeros(numNoiseObs,1);
                 38 indices = d.correction.alpha.indices;
                 39 for m=1:numNoiseObs
                 40   ind = zeros(size(d.antenna0.receiver.utc));
                 41   ind(indices(m,1):indices(m,3)) = 1;
                 42   ind = logical(ind);
                 43   dc = framecut(d, ind);
                 44   % we should have 4 intensity sky values, states on/off for chans 1/2
                 45   dataNoiseOn = d.antenna0.receiver.data((indices(m,4)+5):indices(m,2), [1 9]);
                 46   dataNoiseOff= ...
                 47       d.antenna0.receiver.data([indices(m,1):(indices(m,4)-5) ...
                 48 	(indices(m,2)+5):indices(m,3)], [1 9]);
                 49 
                 50   
                 51   % calculate the opacities from theory 
                 52   % WILL MODIFY LATER TO GET THE OPACITY FROM SKYDIPS
                 53   tau = calcOpacity(5, dc.array.weather.airTemperature+273.15, ...
                 54       dc.array.weather.relativeHumidity, ...
                 55       length(dc.array.weather.relativeHumidity), 1);
                 56   cal.tau(m,1) = mean(tau)./cos( mean(dc.antenna0.servo.el)*pi/180);
                 57   
                 58   % make sure this obs is on a calibrator
                 59   thisSource = unique(dc.antenna0.tracker.source);
                 60   if(length(thisSource)>1)
                 61     cal.whatSource(m,1) = 99;
                 62   else
                 63     [cal.whatSource(m,1) cal.sourceFlux(m,1)] = sourceCorrespondance(thisSource);
                 64   end
                 65   
                 66   if(cal.whatSource(m) ~= 99)
                 67     % we should now determine whether it's on or off the source.
                 68     offsets = unique(dc.antenna0.tracker.horiz_off(:,1));
                 69     if(offsets==0)
                 70       cal.onSource(m,1) = 1;
                 71     end
                 72     cal.time(m,1) = dc.antenna0.receiver.utc(1);
                 73     cal.PnoiseOn(m,:)  = mean(dataNoiseOn);
                 74     cal.PnoiseOff(m,:) = mean(dataNoiseOff);
                 75     cal.varNoiseOn(m,:) = var(dataNoiseOn);
                 76     cal.varNoiseOff(m,:) = var(dataNoiseOff);
                 77   else
                 78     cal.PnoiseOn(m,1:2) = nan;
                 79     cal.PnoiseOff(m,1:2) = nan;
                 80     cal.varNoiseOn(m,1:2) = nan;
                 81     cal.varNoiseOff(m,1:2) = nan;
                 82   end
                 83   
                 84 end
                 85 
                 86 % next we get rid of any observations that are not on a recognized
                 87 % calibrator.
                 88 ind = cal.whatSource == 99;
                 89 
                 90 cal = structcut(cal, ~ind);
                 91 
                 92 if(length(find(cal.onSource==1)) ~= length(find(cal.onSource==0)))
                 93   display('Number of noise events on and off source do not match');
                 94   display('Can not calculate the system temperature');
                 95   tsys = [];
                 96   return;
                 97 end
                 98 
                 99   
                100 % check if we actually have any source
                101 % 18-Mar-2011 (MAS) -- Modified, this crashed if structct() rendered cal no
                102 % longer a structure (which is possible!).
                103 if(~isstruct(cal) || length(cal.whatSource)<1)
                104   display('No Calibrator Source Observations in this schedule');
                105   display('Can not calculate system temperatures');
                106   tsys = [];
                107   return;
                108 end
                109 
                110 
                111 
                112 % ok.  now we start with the first one, find the corresponding match, and
                113 % calculate our y-factors.
                114 index = 1;
                115 for m=1:(length(cal.time)/2)
                116   
                117   % first we find the correspondance
                118   ind =  cal.whatSource==cal.whatSource(2*m-1) & ...   % match the source
                119       cal.onSource == ~cal.onSource(2*m-1);          % if it's on, match an off
                120   f = find(ind);
                121   % search for the closest one in time.
                122   ff = find( abs(cal.time(f) - cal.time(2*m-1)) == min(abs(cal.time(f) - ...
                123       cal.time(2*m-1))));
                124   matchVal = f(ff);
                125 
                126   % ok...so we have our match, [2*m-1, matchVal]
                127   if(cal.onSource(2*m-1))
                128     Psrc_noiseon    = cal.PnoiseOn(2*m-1,:);
                129     Psrc_noiseoff   = cal.PnoiseOff(2*m-1,:);
                130     Pblank_noiseon  = cal.PnoiseOn(matchVal,:);
                131     Pblank_noiseoff = cal.PnoiseOff(matchVal,:);
                132     Vsrc_noiseon    = cal.varNoiseOn(2*m-1,:);
                133     Vsrc_noiseoff   = cal.varNoiseOff(2*m-1,:);
                134     Vblank_noiseon  = cal.varNoiseOn(matchVal,:);
                135     Vblank_noiseoff = cal.varNoiseOff(matchVal,:);
                136   else
                137     Psrc_noiseon    = cal.PnoiseOn(matchVal,:);
                138     Psrc_noiseoff   = cal.PnoiseOff(matchVal,:);
                139     Pblank_noiseon  = cal.PnoiseOn(2*m-1,:);
                140     Pblank_noiseoff = cal.PnoiseOff(2*m-1,:);
                141     Vsrc_noiseon    = cal.varNoiseOn(matchVal,:);
                142     Vsrc_noiseoff   = cal.varNoiseOff(matchVal,:);
                143     Vblank_noiseon  = cal.varNoiseOn(2*m-1,:);
                144     Vblank_noiseoff = cal.varNoiseOff(2*m-1,:);
                145   end    
                146   
                147   % set up our y-values
                148   ynoise = Psrc_noiseoff./Psrc_noiseon;
                149   ycrazy = (Pblank_noiseon - Pblank_noiseoff)./(Psrc_noiseon - ...
                150       Pblank_noiseon);
                151   
                152   varynoise = Vsrc_noiseoff./(Psrc_noiseoff.^2) + Vsrc_noiseon./Psrc_noiseon.^2;
                153   varynoise = varynoise.*ynoise.^2;
                154   varycrazy = (Vblank_noiseon + Vblank_noiseoff)./(Pblank_noiseon - Pblank_noiseoff).^2 + (Vsrc_noiseon + ...
                155       Vblank_noiseoff)./(Psrc_noiseon - Pblank_noiseoff).^2;
                156   varycrazy = varycrazy.*ycrazy.^2;
                157   
                158   % calculate the Tsrc
                159   %apEff = getApEff(d.array.frame.utc(1));
                160   apEff = 0.55;  % 55% theoretical
                161   k_b   = 1.38e-23;
                162   Tsrc  = cal.sourceFlux(2*m-1).*(pi*(6.1/2).^2*apEff)./(2*k_b*1e26);
                163 
                164   % get the mean tau
                165   tauVal = mean(cal.tau([2*m-1 matchVal]));
                166   
                167   % and now we put it all together
                168 
                169   tsys.tnoise(index,:) = ycrazy.*Tsrc.*exp(-tauVal);
                170 
                171   Tcmb = 2.728*exp(-tauVal);
                172   Tatm = mean(dc.array.weather.airTemperature + 273.15).*(1-exp(-tauVal));
                173   Tsky = Tcmb + Tatm;
                174 
                175 %  tsys.rxtemp(index,:) = (2./(1-ynoise)).*(Tsky + ...
                176 %      Tsrc.*exp(-tauVal).*(ynoise-1) + ynoise*tsys.tnoise(index));
                177 %  tsys.val(index,:) = tsys.rxtemp(index,:) + Tcmb + Tatm;
                178 
                179    tsys.val(index,:) = ycrazy.*Tsrc*exp(-tauVal).*ynoise./(1-ynoise) - ...
                180        Tsrc*exp(-tauVal);
                181    tsys.trx(index,:) = tsys.val(index,:) - Tsky;
                182 
                183   tsys.time(index,1) = cal.time(2*m-1);
                184   tsys.source(index,1) = cal.whatSource(2*m-1);
                185   tsys.tnoise(index,:) = ycrazy*Tsrc*exp(-tauVal);
                186   tsys.Tsrc(index,:) = Tsrc;
                187   tsys.tau(index,:) = tauVal;
                188   
                189   tsys.ynoise(index,:) = ynoise;
                190   tsys.ycrazy(index,:) = ycrazy;
                191   
                192   g1 = (cal.PnoiseOn(2*m-1,:) - ...
                193       cal.PnoiseOff(2*m-1,:))./(k_b*tsys.tnoise(index,:).*1e9);
                194   g1dB = 10*log10(g1);
                195   g2 = (cal.PnoiseOn(matchVal,:) - ...
                196       cal.PnoiseOff(matchVal,:))./(k_b*tsys.tnoise(index,:).*1e9);
                197   g2dB = 10*log10(g2);
                198 
                199   thisGain = [cal.time(2*m-1) g1 cal.time(matchVal) g2];
                200   thisGaindB = [cal.time(2*m-1) g1dB cal.time(matchVal) g2dB];
                201   
                202   tsys.gain(index,:) = thisGain;
                203   tsys.gaindB(index,:) = thisGaindB;
                204   
                205   
                206   % and now, errors
                207   ddycrazy = Tsrc.*exp(tauVal).*ynoise./(1-ynoise);
                208   ddy      = ycrazy.*Tsrc*exp(tauVal)./(1-ynoise).^2;
                209   varTot   = ddycrazy.^2.*varycrazy + ddy.^2.*varynoise;
                210   tsys.err(index,:) = sqrt(varTot);
                211 
                212   
                213   
                214   
                215   
                216   index = index+1;
                217 end
                218 
                219 return;

Other subfunctions in this file are not included in this listing.