tsysWrapper

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

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

    reduceData function that calls functions to calculate tsys

    sjcm

   3/30/2011 sjcm:  saving plots without display
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function d = tsysWrapper(d, plotparams, parm, field,thedate_day_portion);
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 %  function d = tsysWrapper(d, plotparams, parm, field);
0006 %
0007 %    reduceData function that calls functions to calculate tsys
0008 %
0009 %    sjcm
0010 %
0011 %   3/30/2011 sjcm:  saving plots without display
0012 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0013 
0014 % first we calculate the system temperatures.
0015 
0016 % this calculation requires four measurements
0017 %  Sky obs, no noise
0018 %  Sky obs, noise
0019 %  Calibrator source, no noise
0020 %  Calibrator source, noise
0021 % if any are missing, the return value is an empty set, otherwise it's
0022 
0023 
0024 [a, b] = logcal(d, 'pointing');
0025 if(b~=1)
0026   display('You need to calculate the pointing corrections first');
0027   return;
0028 end
0029 [a, b] = logcal(d, 'alpha');
0030 if(b~=1)
0031   display('You need to apply alpha');
0032   return;
0033 end
0034 [a, b] = logcal(d, 'stokes');
0035 if(b~=1)
0036   display('You need to calculate stokes');
0037   return;
0038 end
0039 [a, b] = logcal(d, 'noise');
0040 if(b~=1)
0041   display('You need to calculate the noise diode temps');
0042   return;
0043 end
0044 clear a;
0045 
0046 tsys_values_filename = ['tsys_values_',thedate_day_portion,'.txt'];
0047 
0048 % calculating the system temperature at this stage is pretty easy:  you get
0049 % the value from the alpha database, and multiply it by your noise diode
0050 % temperature
0051 x = readAlphaDatabase();
0052 times   = x{1};
0053 alpha   = x{9};
0054 gain    = x{8};
0055 r       = x{12};
0056 TsysTND = x{13};
0057 flags   = x{17};
0058 
0059 % get the noise diode temperatures on the same timescale
0060 f = find(times >= d.antenna0.receiver.utc(1) & times <= ...
0061     last(d.antenna0.receiver.utc));
0062 hr = 1;
0063 doSave = 1;
0064 while(length(f)<2)
0065   f = find(times >= d.antenna0.receiver.utc(1) - hr/24 & times <= ...
0066       last(d.antenna0.receiver.utc) + hr/24);
0067   hr = hr+1;
0068   doSave = 0;
0069 end
0070 eval(sprintf('display(''Found tsys within %d hour of track'');', hr-1));
0071 
0072 tnd   = interp1(d.antenna0.receiver.utc, d.antenna0.receiver.noise, times(f));
0073 Tsys  = TsysTND(f,:).*tnd;
0074 flags = repmat( flags(f,:), [1 2]);;
0075 gain  = gain(f,:);
0076 
0077 % next we flag and plot
0078 parFlag = checkpar(parm, 'tsys');
0079 if(parFlag)
0080   flagParams = [parm.tsys.valuehi parm.tsys.valuelo];
0081 else
0082   display('Using Default 1st stage RFI flagging parameters');
0083   flagParams = [100 10];
0084 end
0085 tsysFlag = Tsys > flagParams(1) & Tsys < flagParams(2);
0086 % set our variables and plot the data.
0087 d.correction.tsys.time = times(f);
0088 d.correction.tsys.val  = Tsys;
0089 d.correction.tsys.gain = gain;
0090 d.correction.tsys.flag = tsysFlag | flags;
0091 
0092 tsysflags = d.correction.tsys.flag;
0093 setPlotDisplay(plotparams.plot);
0094 display('Plotting data');
0095 [d tsysflags] = packd(d, tsysflags, 'none', 'tsys', plotparams, 'System Temperature Plots', field);
0096 d.correction.tsys.flag = tsysflags;
0097 
0098 % now we apply the flags, and interpolate the system temperature
0099 for m=1:2
0100   x = times(f);
0101   y = Tsys(:,m);
0102   x(tsysflags(:,m)) = [];
0103   y(tsysflags(:,m)) = [];
0104   if(length(x)==0)
0105     % no data to make plots of
0106     tsysInterp(:,m) = nan(size(d.antenna0.receiver.utc));
0107   elseif(length(x)==1)
0108     tsysInterp(:,m) = ones(size(d.antenna0.receiver.utc))*y;
0109   else
0110     tsysInterp(:,m) = interp1(x, y, d.antenna0.receiver.utc, 'linear', ...
0111     'extrap');
0112   end
0113 end
0114 indbad = tsysInterp > flagParams(1) | tsysInterp < flagParams(2);
0115 flag2set(:,1) = indbad(:,1);
0116 flag2set(:,3) = indbad(:,2);
0117 flag2set(:,2) = indbad(:,1) | indbad(:,2);
0118 d = setNewFlag(d, flag2set, 'tsys');
0119 d.antenna0.receiver.tsys = tsysInterp;
0120 
0121 % next we write things to disk -- all we want is time, Tsys, flag.
0122 thisTsys = [d.correction.tsys.time d.correction.tsys.val tsysflags];
0123 
0124 [home,installeddir]=where_am_i();
0125 
0126 % make sure file is up-to-date
0127 parFlag = checkpar(parm, 'autosave');
0128 if(parFlag)
0129   if(parm.autosave.flag & doSave)
0130     %system('cvs update constants/tsys_values.txt');
0131     %allTsys  = load([home,'/',installeddir,'/constants/tsys_values.txt']);
0132     %allTsys  = [allTsys; thisTsys];
0133     save([home,'/',installeddir,'/constants/tsys_values.txt'], 'thisTsys', '-ascii', '-single','-append')
0134     save([home,'/',installeddir,'/constants/',tsys_values_filename], 'thisTsys', '-ascii', '-single','-append') 
0135 
0136     %system('cvs commit -m "update" constants/tsys_values.txt');
0137   
0138   end
0139 end
0140 
0141 d = logcal(d, 'tsys');
0142 
0143 return;
0144 
0145