noiseWrapper

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

  function d = noiseWrapper(d, plotparams, parm, field);

    reduceData function that calls functions to calculate tnoise if there
    are calibrator observations.  if there are no calibrator observations,
    it loads the noise diode data from the archive.

    sjcm

   08/02/2012 - sjcm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function d = noiseWrapper(d, plotparams, parm, field,thedate_day_portion)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 %  function d = noiseWrapper(d, plotparams, parm, field);
0006 %
0007 %    reduceData function that calls functions to calculate tnoise if there
0008 %    are calibrator observations.  if there are no calibrator observations,
0009 %    it loads the noise diode data from the archive.
0010 %
0011 %    sjcm
0012 %
0013 %   08/02/2012 - sjcm
0014 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0015 
0016 % first we calculate the noise temperatures.
0017 
0018 % this calculation requires four measurements
0019 %  Sky obs, no noise
0020 %  Sky obs, noise
0021 %  Calibrator source, no noise
0022 %  Calibrator source, noise
0023 % if any are missing, the return value is an empty set, otherwise it's
0024 
0025 if(size(d.antenna0.receiver.data,2)~=8)
0026   display('You need to run alpha and stokes before this step');
0027   return;
0028 end
0029 
0030 [a, b] = logcal(d, 'pointing');
0031 if(b~=1)
0032   display('You need to calculate the pointing corrections first');
0033 end
0034 [a, b] = logcal(d, 'alpha');
0035 if(b~=1)
0036   display('You need to apply alpha');
0037 end
0038 [a, b] = logcal(d, 'stokes');
0039 if(b~=1)
0040   display('You need to calculate stokes');
0041 end
0042 clear a;
0043 
0044 
0045 noise_diode_values_offint_filename = ['noise_diode_values_offint_',thedate_day_portion,'.txt'];
0046 noise_diode_values_ptcross_filename = ['noise_diode_values_ptcross_',thedate_day_portion,'.txt'];
0047 calsource_values_filename = ['calsource_values_',thedate_day_portion,'.txt'];
0048 nonlin_values_filename = ['nonlin_values_',thedate_day_portion,'.txt'];  
0049 
0050 % calculate the tnoise
0051 display('in noise wrapper');
0052 tnoiseVals = calcTnoise_v3(d);
0053 method = 1;
0054 
0055 if(isempty(tnoiseVals))
0056   method = 0;
0057   d.correction.tsys.val = [];
0058   display('No System temperature calculated from your data set');
0059   display('Reading values from the archive.');
0060 
0061   if(issubfield(d, 'correction', 'pointing'))
0062     if(~isempty(d.correction.pointing))
0063       display('Will try to calculate it from the old method');
0064       tnoiseVals = calcTsys(d);
0065       method = 2;
0066     end
0067   else
0068     tnoiseVals = [];
0069   end
0070   if(isempty(tnoiseVals))
0071     method = 0;
0072     display('Obtaining Values from the Archive');
0073     [tnoiseVals, noiseInt, noiseIntErr] = getNoiseDiodeTemps(d);  
0074     d.correction.survey = 1;
0075   else
0076     d.correction.survey = 0;
0077   end
0078 else
0079   d.correction.survey = 0;
0080 end
0081 
0082 % next we flag
0083 parFlag = checkpar(parm, 'noise');
0084 if(parFlag)
0085   flagParams = [parm.noise.valuehi parm.noise.valuelo];
0086 else
0087   display('Using Default Noise diode flagging values');
0088   flagParams = [5 2];
0089 end
0090 
0091 if method == 1 
0092    [noiseflag] = flagNoiseV2(tnoiseVals, flagParams);
0093 else [noiseflag] = flagNoise(tnoiseVals, flagParams);
0094 end
0095 % the noiseflag are the bad values.
0096 f = find(noiseflag);
0097 display(sprintf('Tnoise Flagging: %i of %i measurements unreliable', length(f), ...
0098     length(noiseflag(:))));
0099 
0100 d.correction.tnoise = tnoiseVals;
0101 d.correction.tnoise.flags = noiseflag;
0102 
0103 % if all our data is bad, than this data is useless for calibration
0104 if(all(noiseflag(:))) 
0105   display('ATTENTION --- ATTENTION --- ATTENTION');
0106   display('No good tnoise measurements in your data');
0107   display('This data set will not be used for calibration');
0108   d.correction.tnoise.allbad = 1;
0109   % reading the noise values from the archive
0110   method = 0;
0111   [tnoiseVals, noiseInt, noiseIntErr] = getNoiseDiodeTemps(d);  
0112   [noiseflag] = flagNoise(tnoiseVals, flagParams);
0113   d.correction.tnoise = tnoiseVals;
0114   d.correction.tnoise.flags = noiseflag;
0115   d.correction.tnoise.allbad = 1;
0116 else
0117   d.correction.tnoise.allbad = 0;
0118 end
0119 
0120 % next we plot everything!
0121 setPlotDisplay(plotparams.plot);
0122 % all flagging has already happened in the tsys step, so we just plot.
0123 display('Plotting data');
0124 flags = d.correction.tnoise.flags;
0125 
0126 [d flags] = packd(d, flags, 'none', 'noise', plotparams, 'Noise Diode Temperature Plots', field);
0127 d.correction.tnoise.flags = flags;
0128 d.correction.tnoise.method = method;
0129 
0130 % last thing we do is interpolate the noise diode values onto our timestream
0131 % data -- this is already done for the survey observations.
0132 if(d.correction.survey || method==0)
0133   d.antenna0.receiver.noise = noiseInt; 
0134   d.antenna0.receiver.noiseErr = noiseIntErr; 
0135 else
0136   for m=1:2
0137     % let's flag also we only want tnd2
0138     y = d.correction.tnoise.Tnd2(:,m);
0139     yErr = d.correction.tnoise.TndErr2(:,m);
0140     y1 = y;
0141     y1Err = yErr;
0142     y1(d.correction.tnoise.flags(:,m)) = [];
0143     y1Err(d.correction.tnoise.flags(:,m)) = [];
0144     if(length(y1)>=1)
0145        [wmean error_wm] = weighted_mean(y1,y1Err);
0146        noiseInt(:,m) = repmat(wmean, length(d.antenna0.receiver.utc), 1);
0147        noiseIntErr(:,m) = repmat(error_wm, length(d.antenna0.receiver.utc), 1);
0148     elseif(isempty(y1))
0149       display('All values flagged for this channel');
0150       noiseInt(:,m) = nan(size(d.antenna0.receiver.utc)); 
0151       noiseIntErr(:,m) = nan(size(d.antenna0.receiver.utc)); 
0152     end
0153   end
0154   d.antenna0.receiver.noise = noiseInt; 
0155   d.antenna0.receiver.noiseErr = noiseIntErr;
0156 
0157   % next, let's save the nonlinearity factors
0158   d.correction.nonlin.time   = d.correction.tnoise.time;
0159   d.correction.nonlin.source = d.correction.tnoise.source;
0160   d.correction.nonlin.val    = d.correction.tnoise.nonlin;
0161 end
0162 
0163 % next, we save these noise value temperatures
0164 % what do we want to save is:
0165 if(d.correction.survey==0 & d.correction.tnoise.allbad==0)
0166   % we calculated the noise diode temperature here, and now want to save it
0167   % we want to save the time, the source, elevation, position angle, ...
0168   %  Tnd (2chan), Tnd err (2chan), and flags
0169   thisNoise = [d.correction.tnoise.time d.correction.tnoise.source ...
0170     d.correction.tnoise.elev d.correction.tnoise.pa ...
0171     d.correction.tnoise.Tnd d.correction.tnoise.TndErr ...
0172     d.correction.tnoise.flags ones(size(d.correction.tnoise.time)).*method];
0173 
0174   [home,installeddir]=where_am_i();
0175   
0176   if (any(isnan(thisNoise)))
0177       disp('WARN FOUND A NAN IN NOISE WRAPPER');     
0178   end;
0179 
0180   % make sure file is up-to-date
0181   parFlag = checkpar(parm, 'autosave');
0182   if(parFlag)
0183     if(parm.autosave.flag)
0184       %system('cvs update constants/noise_diode_values.txt');
0185       %allNoise  = load([home,'/',installeddir,'/constants/noise_diode_values_offint.txt']);
0186       %allNoise  = [allNoise; thisNoise];
0187       save([home,'/',installeddir,'/constants/',noise_diode_values_offint_filename], 'thisNoise', '-ascii', '-single','-append') 
0188       save([home,'/',installeddir,'/constants/noise_diode_values_offint.txt'], 'thisNoise', '-ascii', '-single','-append')
0189      
0190 %      system('cvs commit -m "update" constants/noise_diode_values.txt');
0191     end
0192   end
0193   
0194   if(method~=2)
0195     thisNoise = [d.correction.tnoise.time d.correction.tnoise.source ...
0196       d.correction.tnoise.elev d.correction.tnoise.pa ...
0197       d.correction.tnoise.Tnd2 d.correction.tnoise.TndErr2 ...
0198       d.correction.tnoise.flags ones(size(d.correction.tnoise.time)).*method];
0199     [home,installeddir]=where_am_i();
0200     % make sure file is up-to-date
0201     parFlag = checkpar(parm, 'autosave');
0202     if(parFlag)
0203       if(parm.autosave.flag)
0204     %      system('cvs update constants/noise_diode_values.txt');
0205     %allNoise  = load([home,'/',installeddir,'/constants/noise_diode_values_ptcross.txt']);
0206     %allNoise  = [allNoise; thisNoise];
0207     save([home,'/',installeddir,'/constants/noise_diode_values_ptcross.txt'], 'thisNoise', '-ascii', '-single','-append')
0208         save([home,'/',installeddir,'/constants/',noise_diode_values_ptcross_filename], 'thisNoise', '-ascii', '-single','-append')  
0209     %      system('cvs commit -m "update" constants/noise_diode_values_ptcross.txt');
0210       end
0211     end
0212   end
0213 end
0214 
0215 
0216 % and write the nonlin data out
0217 if(d.correction.survey==0 & d.correction.tnoise.allbad==0)
0218   thisNonlin = [d.correction.nonlin.time d.correction.nonlin.source ...
0219     d.correction.nonlin.val(:,1:2)];
0220   if(parFlag)
0221     if(parm.autosave.flag)
0222 %      system('cvs update constants/nonlin_values.txt');
0223       %allNonlin  = load([home,'/',installeddir,'/constants/nonlin_values.txt']);
0224       %allNonlin  = [allNonlin; thisNonlin];
0225       save([home,'/',installeddir,'/constants/nonlin_values.txt'], 'thisNonlin', '-ascii', '-single','-append')
0226       save([home,'/',installeddir,'/constants/',nonlin_values_filename], 'thisNonlin', '-ascii', '-single','-append')
0227 %      system('cvs commit -m "update" constants/nonlin_values.txt');
0228     end
0229   end
0230 end
0231 
0232 % and write out taua data
0233 if(d.correction.survey==0 & d.correction.tnoise.allbad==0 & method==1)
0234   % ACT - 8/10/2013 - write out data for all calibration sources
0235   % Save in file constants/calsource_values.txt
0236     %f = find(d.correction.tnoise.source == 2);
0237   %if(~isempty(f))
0238     % [ time, az, el, pa, I1, I2, Q, U];
0239     % still need to check if this is the case.
0240     thisTauA = [d.correction.tnoise.time  d.correction.tnoise.az ...
0241       d.correction.tnoise.elev d.correction.tnoise.pa ...
0242       d.correction.tnoise.src_noiseoff d.correction.tnoise.source];
0243 
0244     if(parFlag)
0245       if(parm.autosave.flag)
0246         disp('Saving cal source data')
0247     %system('cvs update constants/nonlin_values.txt');
0248     %allTauA  = load([home,'/',installeddir,'/constants/calsource_values.txt']);
0249         %allTauA  = [allTauA; thisTauA];
0250         save([home,'/',installeddir,'/constants/calsource_values.txt'], 'thisTauA', '-ascii', '-single','-append')
0251         save([home,'/',installeddir,'/constants/',calsource_values_filename], 'thisTauA', '-ascii', '-single','-append') 
0252           
0253      end
0254     end    
0255   %end
0256 end
0257 
0258 d = logcal(d, 'noise');
0259 
0260 return;
0261