calculateTau

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

  function tau = calculateTau(d)

  calculates the opacity values from our track given a noise diode.

  sjcm

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

0001 function tau = calculateTau(d)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 %  function tau = calculateTau(d)
0006 %
0007 %  calculates the opacity values from our track given a noise diode.
0008 %
0009 %  sjcm
0010 %
0011 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0012 
0013 Tcmb = 2.728;
0014 
0015 % first we see if there is any skydip data
0016 if(length(find(d.index.skydip.fast))==0)
0017   display('calculateTau:: no sky dips in your data');
0018   tau = [];
0019   return;
0020 end
0021 
0022 % select the data (finding the start and stop points)
0023 [si ei] = findStartStop(d.index.skydip.fast);
0024 
0025 
0026 allFlag = [];
0027 % Calculation from Skydips using the noise diode temperature
0028 for m=1:length(si)
0029   ind = zeros(size(d.index.skydip.fast));
0030   ind(si(m):ei(m)) = 1;
0031   ind = logical(ind);
0032   dc = framecut(d, ind);
0033   time(m,1)    = mean(dc.array.frame.utc);
0034   Tamb(m,1)    = mean(dc.array.weather.airTemperature + 273.15);
0035   TambErr(m,1) = std(dc.array.weather.airTemperature + 273.15);
0036   
0037   % next to calculate the slope, let us just get the data where we are
0038   % actually scanning in elevation % scan down at 0.5deg/s
0039   del = deriv(dc.antenna0.servo.el);
0040   ind = del < -0.4/100;
0041   elvals = dc.antenna0.servo.el(ind);
0042   Ivals  = dc.antenna0.receiver.data(ind, [1 8]);
0043   Ivals  = Ivals.*2.*dc.antenna0.receiver.noise(ind,:);
0044 
0045   thisflag = [ 0 0 ];
0046   indflag = dc.flags.fast(ind,[1 3]);
0047   for mm=1:2
0048     thisInd = ~indflag(:,mm);
0049     elfits = 1./sind(elvals(thisInd));
0050     Ifit   = Ivals(thisInd,mm);
0051     if(length(Ifit) < 10)
0052       mx = nan;
0053       b = nan;
0054       thisflag(mm) = 1;
0055     else
0056       % next we fit for a slope of the line
0057       [mx b] = linfit(elfits, Ifit);
0058     end
0059     slopes(m,mm) = mx;
0060     if( (max(elfits) - min(elfits)) < 0.12)
0061       thisflag(mm) = 1;
0062     end
0063     if(min(elfits) > 60)
0064       thisflag(mm) = 1;
0065     end
0066     % flag if the fit residuals are bad too
0067     fitline = mx.*elfits + b;
0068     resid   = Ifit - fitline;
0069     residrat= abs(rms(resid)./mean(Ivals)*100);
0070     if(residrat>0.5)
0071       thisflag(mm) = 1;
0072     end
0073     tauVal(m,mm) = slopes(m,mm)./(Tamb(m) - Tcmb);
0074   end
0075   allFlags(m,:) = thisflag;
0076   el(m,1) = mean(dc.antenna0.servo.el);
0077   az(m,1) = mean(dc.antenna0.servo.az);
0078 end
0079 
0080 tauErr   = zeros(size(tauVal));
0081 tau.time = time;
0082 tau.Tamb = Tamb;
0083 tau.values  = tauVal;
0084 tau.elev = el;
0085 tau.az   = az;
0086 tau.error = tauErr;
0087 tau.flag  = allFlags;
0088 tau.source = zeros(size(el));
0089 
0090 return;
0091