calcTauNormCal

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

  function tauVals = calcTauNormCal(d)
   
    calculates the opacity by normalizing power to calibrator observations,
    as in Section 4.1 of Stephen's abscal memo.

  sjcm

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

0001 function tauVals = calcTauNormCal(d)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 %  function tauVals = calcTauNormCal(d)
0006 %
0007 %    calculates the opacity by normalizing power to calibrator observations,
0008 %    as in Section 4.1 of Stephen's abscal memo.
0009 %
0010 %  sjcm
0011 %
0012 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0013 
0014 % select the data (finding the start and stop points)
0015 [sdip edip] = findStartStop(d.index.skydip.fast);
0016 
0017 [dd indRx]  = cutObs(d, 'calibrator', 'only');
0018 [ssource esource] = findStartStop(indRx);
0019 
0020 if(isempty(ssource))
0021   display('calcTauNormCal:: No suitable calibrator observations');
0022   display('calcTauNormCal:: Can not use this method for calculating tau');
0023   tauVals = [];
0024   return;
0025 end
0026 
0027 % first we get the power measurements for each of the source observations,
0028 % along with their flux and time.
0029 for m=1:length(ssource)
0030   ind = zeros(size(indRx));
0031   ind(ssource(m):esource(m)) = 1;
0032   ind = logical(ind);
0033   dc  = framecut(d, ind);
0034   % check if it's only one source
0035   if(length(unique(dc.antenna0.tracker.source))~=1)
0036     source.power(m,:) = [0 0];
0037     source.num(m,1)   = 99;
0038     source.flux(m,1)  = 0;
0039     source.time(m,1)  = 0;
0040     source.elev(m,1)  = 0;
0041     source.az(m,1)    = 0;
0042   else
0043     source.power(m,:) = mean(dc.antenna0.receiver.data(:, [1 9]));
0044     [source.num(m,1) source.flux(m,1)] = ...
0045     sourceCorrespondance(unique(dc.antenna0.tracker.source));
0046     source.time(m,1)  = mean(dc.antenna0.receiver.utc);
0047     source.elev(m,1)  = mean(dc.antenna0.servo.el);
0048     source.az(m,1)    = mean(dc.antenna0.servo.az);
0049   end
0050 end
0051 
0052 % throw out the ones that are zero.
0053 ind = source.time==0;
0054 source = structcut(source, ~ind);
0055 
0056 % next we split up the skydip data and find the source that is closest to
0057 % each skydip
0058 for m=1:length(sdip)
0059   ind = zeros(size(indRx));
0060   ind(sdip(m):edip(m)) = 1;
0061   ind = logical(ind);
0062   dc  = framecut(d, ind);  
0063   
0064   % elevation of source is the first entry
0065   elevDip = dc.antenna0.servo.el(1);
0066   timeDip = dc.antenna0.servo.utc(1);
0067   azDip   = dc.antenna0.servo.az(1);
0068   
0069   % so we look for the corresponding source that is closest in elevation and
0070   % azimuth
0071   elDiff  = source.elev-elevDip;
0072   indElev = abs(elDiff) < 5;
0073   azDiff  = source.az - azDip;
0074   indAz   = (abs(azDiff) > 3 & abs(azDiff)<10) | (abs(azDiff-360) > 3 & ...
0075       abs(azDiff-360)<10) | (abs(azDiff+360) > 3 & abs(azDiff+360)<10); 
0076   indCoord= indElev & indAz;
0077   
0078   
0079   % of those criteria, we look for the one that was closest in time.
0080   timeDiff = source.time - timeDip;
0081   f  = find(indCoord);
0082   ff = find(abs(timeDiff(f)) == min(abs(timeDiff(f))));
0083   f  = f(ff);
0084   
0085   % check that the time is not way off (15 minutes)
0086   if(abs(timeDiff(f)*24*60) > 15)  
0087     f = [];
0088   end
0089   
0090   if(~isempty(f))
0091     % we have a good sky dip an corresponding power measurement
0092     skydip.data{m,1} = dc;
0093     skydip.power(m,:) = source.power(f,:);
0094     skydip.flux(m,:)  = source.flux(f,:);
0095     skydip.el0(m,:)   = source.elev(f,:);
0096   else
0097     skydip.data{m,1} = [];
0098     skydip.power(m,:) = [0 0];
0099     skydip.flux(m,:)  = 0;
0100     skydip.el0(m,:)   = 0;
0101   end
0102 end
0103 
0104 % throw out values where the power is zero
0105 ind = skydip.flux == 0;
0106 skydip = structcut(skydip, ~ind);
0107 if(isempty(skydip.flux))
0108   display('calcTauNormCal:: No suitable data for this method');
0109   tauVals = [];
0110   return;
0111 end
0112 
0113 
0114 % load apEff and k_b, Tcmb, etc etc
0115 telescope_constants
0116 
0117 keyboard;
0118 % normalize each skydip and fit a line.
0119 for m=1:length(skydip.flux)
0120   
0121   % calculate constants needed first
0122   beta = 1./sind(skydip.el0(m));
0123   Tsrc = skydip.flux(m).*(areaEff)./(2*k_b*1e26);
0124   Tamb = mean(skydip.data{m}.array.weather.airTemperature) + 273.15;
0125   
0126   % get the data to fit
0127   y = skydip.data{m}.antenna0.receiver.data(:, [1 9]);
0128   x = 1./sind(skydip.data{m}.antenna0.servo.el);
0129   
0130   % normalize the data
0131   norm = repmat(skydip.power(m,:), [size(y,1) 1]);
0132   y = y./norm;
0133 
0134   % calculate some constants
0135   
0136   for mm=1:2
0137     thisX = x;
0138     thisY = y(:,mm);
0139     thisX(isnan(thisY)) = [];
0140     thisY(isnan(thisY)) = [];
0141     
0142     if(~isempty(thisY))
0143       [slope(m,mm) intercept(m,mm)] = linfit(thisX, thisY);
0144     else
0145       slope(m,mm) = nan;
0146       intercept(m,mm) = nan;
0147     end
0148   end
0149                        
0150 
0151   tau0 = -Tsrc./(Tamb - Tcmb*(intercept - 1)./(slope) - beta*Tsrc);
0152 end
0153 
0154