Description of getNoiseDiodeTemps

0001 function [noiseVals, noiseInterp, noiseInterpErr] = getNoiseDiodeTemps(d)
0002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0003 %
0004 %  function [noiseVals noiseInterp] = getNoiseDiodeTemps(d, primCal)
0005 %
0006 %
0007 %   function that reads an ASCII text file, re-orders everything in that
0008 %   file, calculates the flux for secondary calibrators, determines the
0009 %   noise diode temperature value for the times in your current schedule,
0010 %   and interpolates the output for you
0011 %
0012 %   sjcm
0013 %
0014 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0015 
0016 if(nargin<2)
0017   primCal = [1,2,8];  % casa
0018   % other prims are 1: Cyga, and 2: taua
0019 end
0020 
0021 filename = 'constants/noise_diode_values_ptcross.txt';
0022 
0023 % file is in form:
0024 %[time source elev pa Tnd1 Tnd2 Tnd1err Tnd2err flag1 flag2 method];
0025 % only want to use method 1!
0026 
0027 tnoise = load(filename);
0028 
0029 % first let us cut on method
0030 indmethod = tnoise(:,11) == 1;
0031 tnoise = tnoise(indmethod,:);
0032 tnoise(:,11)  = [];
0033 
0034 
0035 % next let us re-order everything
0036 tDate  = tnoise(:,1);
0037 % first put things in chronological order
0038 [tDateChron, sortOrder] = sort(tDate);
0039 tChron     = tnoise(sortOrder,:);
0040 
0041 % redefine the variables in chronologiacl order
0042 tDate     = tDateChron;
0043 tnoise    = tChron;
0044 
0045 % let us get rid of the repeats - always trust the lastest one
0046 reps   = find(diff(tDate)==0);
0047 tDate(reps,:) = [];
0048 tnoise(reps,:) = [];
0049 
0050 % redefine variables for ease to work with
0051 sources   = tnoise(:,2);
0052 times     = tDate;
0053 noiseVals = tnoise(:,5:10);  % Tnd1 Tnd2 Tnd1err Tnd2err flag1 flag2
0054 tnd       = noiseVals(:,1:2);
0055 tndErr    = noiseVals(:, 3:4);
0056 flags     = noiseVals(:,5:6)  | isnan(tnd);
0057 
0058 % first we want to flag the noise diode temperature and interpolate their
0059 % values for the calibrators.
0060 
0061 indPrim1 = find(sources == primCal(1));
0062 indPrim2 = find(sources == primCal(2));
0063 indPrim3 = find(sources == primCal(3));
0064 indPrim = cat(1, indPrim1, indPrim2, indPrim3);
0065 sources = sources(indPrim);
0066 times =times(indPrim);
0067 noiseVals = noiseVals(indPrim, 1:6);
0068 tnd = tnd(indPrim, 1:2);
0069 tndErr = tndErr(indPrim, 1:2);
0070 flags =flags(indPrim, 1:2);
0071  
0072 % next...we find values which are near our data set:
0073 if(isstruct(d))
0074   ind = times > d.array.frame.utc(1) - 2/24  & times < ...
0075       last(d.array.frame.utc) + 2/24;
0076   finc= find(ind);
0077   % if nothing within a few hours, keep getting more!
0078   hr = 73;
0079   while(isempty(finc))
0080     %display(' getNoiseDiodeTemps:: WARNING WARNING WARNING');
0081     %eval(sprintf('display(''getNoiseDiodeTemps:: No noise diode within %d hours of track'')', (hr-1)));
0082     ind =  times > d.array.frame.utc(1) - hr/24  & times < ...
0083     last(d.array.frame.utc) + hr/24;
0084     finc = find(ind);
0085     hr   = hr+1;
0086   end
0087   eval(sprintf('display(''getNoiseDiodeTemps:: Nearest noise diode is within %d hours of track'')', (hr-1)));
0088 else
0089   ind = times > d(1) - 2/24  & times < ...
0090       last(d) + 2/24;
0091   finc= find(ind);
0092   % if nothing within a few hours, keep getting more!
0093   hr = 73;
0094   while(isempty(finc))
0095     %display('getNoiseDiodeTemps:: WARNING WARNING WARNING');
0096     %eval(sprintf('display(''getNoiseDiodeTemps::No noise diode within %d hours of track'')', (hr-1)));
0097     ind =  times > d(1) - hr/24  & times < ...
0098     last(d) + hr/24;
0099     finc = find(ind);
0100     hr   = hr+1;
0101   end
0102   eval(sprintf('display(''getNoiseDiodeTemps:: Nearest noise diode is within %d hours of track'')', (hr-1)));
0103   if ( (hr-1) > 23)
0104     disp('getNoiseDiodeTemps:: WARN - Nearest noise diode event is greater than 24 hours away from track')
0105   end 
0106 end
0107 
0108   
0109   
0110 if(~isempty(finc))
0111   tnd.time = times(finc);
0112   tnd.source = sources(finc);
0113   tnd.elev   = tnoise(finc,3);
0114   tnd.pa     = tnoise(finc,4);
0115   tnd.Tnd    = noiseVals(finc,1:2);
0116   tnd.TndErr = tndErr(finc, 1:2);
0117   tnd.flags  = flags(finc, 1:2);
0118 else
0119   error('getNoiseDiodeTemps:: You have no noise diode info in your archive');
0120 end
0121 % lastly, we interpolate our primary calibrator noise diode values onto our
0122 % data.
0123 
0124 acceptable1 = find(tnd.flags(:,1) == 0);
0125 acceptable2 = find(tnd.flags(:,2) == 0);
0126 using1 = tnd.Tnd(acceptable1, 1);
0127 using2 = tnd.Tnd(acceptable2, 2);
0128 usingErr1 = tnd.TndErr(acceptable1, 1);
0129 usingErr2 = tnd.TndErr(acceptable2, 2);
0130 
0131 [wmean1 error_wm1] = weighted_mean(using1,usingErr1);
0132 [wmean2 error_wm2] = weighted_mean(using2,usingErr2);
0133 if(isnan(wmean1))
0134   wmean1 = 3;
0135   error_wm1 = inf;
0136 end
0137 if(isnan(wmean2))
0138   wmean2 = 3;
0139   error_wm2 = inf;
0140 end
0141 
0142 
0143 
0144 if(isstruct(d))
0145   noiseInterp(:,1) = repmat(wmean1, length(d.antenna0.receiver.utc), 1);
0146   noiseInterpErr(:,1) = repmat(error_wm1, length(d.antenna0.receiver.utc), 1);
0147   noiseInterp(:,2) = repmat(wmean2, length(d.antenna0.receiver.utc), 1);
0148   noiseInterpErr(:,2) = repmat(error_wm2, length(d.antenna0.receiver.utc), 1);
0149 else
0150   noiseInterp = [wmean1 wmean2];
0151   noiseInterpErr = [error_wm1 error_wm2];
0152 end
0153 
0154   
0155 noiseVals = tnd;
0156 
0157 timeM = mean(tnd.time);
0158   %
0159   % OXFORD - can you uncommment and save these values please
0160   % JL Commented out - this file location is now invalid (moved to mammoth)
0161   %fid = fopen('//mnt/usersVol1/cbassvisitor1/cbass_analysis/MelsTests/TestingNDs.txt','a');
0162   %fprintf(fid, '%10.4f %6.4f %5.4f %6.4f %5.4f \n', timeM, wmean1, error_wm1, wmean2, error_wm2);
0163   %fclose(fid);
0164 
0165 return;
0166 
0167 
0168 % plot if we want to.
0169 
0170 
0171 %marks = {'r.', 'g.', 'k.', 'b.', 'c.', 'y.', 'm.'};
0172 %full_flux(full_flags>0) = nan;
0173 
0174 % let us s make a plot
0175 %unCals = unique(sources);
0176 %for m=1:length(unCals)
0177 %  f = find(sources == unCals(m));
0178 %  plot(times(f), full_flux(f,1), marks{m});
0179 %  hold on
0180 %  thisName = getSourceName(unCals(m));
0181 %  fullname{m} = thisName;
0182 %end
0183 %legend(fullname);
0184 
0185 
0186
getNoiseDiodeTemps

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
