stare_track

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  function stare_track(d)

   this function should analyize a stare track of a source.  In particular,
   it will take your data (d), do some flagging, and then make a plot of
   the noise as a function of time (in intervalMinutes), compared to
   theoretical predictions (that the noise goes down as 1/sqrt(time),
   normalized to the first data point).

   this track does a noise observation at the beginning, and then just
   stares at the  source.  in the future, we should put an observation of
   the noise source at regular time intervals.
  
  sjcm

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

0001 function stare_track(d, intervalMinutes)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 %  function stare_track(d)
0006 %
0007 %   this function should analyize a stare track of a source.  In particular,
0008 %   it will take your data (d), do some flagging, and then make a plot of
0009 %   the noise as a function of time (in intervalMinutes), compared to
0010 %   theoretical predictions (that the noise goes down as 1/sqrt(time),
0011 %   normalized to the first data point).
0012 %
0013 %   this track does a noise observation at the beginning, and then just
0014 %   stares at the  source.  in the future, we should put an observation of
0015 %   the noise source at regular time intervals.
0016 %
0017 %  sjcm
0018 %
0019 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0020 
0021 % separate the noise source and the on-source data
0022 indNoise = bitsearch(d.antenna0.receiver.flags, 2, 'any');
0023 d_noise  = framecut(d, indNoise);
0024 
0025 % on-source data
0026 indSrc  = bitsearch(d.array.frame.features, 12, 'any');
0027 % interpolate it onto the size of data
0028 indSrc  = interp1(d.array.frame.utc, indSrc, d.antenna0.receiver.utc, 'nearest');
0029 indSrc(isnan(indSrc)) = 0;
0030 indSrc  = indSrc & ~indNoise;
0031 d_src   = framecut(d, indSrc);
0032 
0033 % first we check for crazy RFI (flag the data a bit)
0034 % rudimentary checking for RFI (Matthew will improve)
0035 % first check I'm doing is basically a box-car average/rms to find regions
0036 % where the values peak relative to the trend.  this is a decent zeroth pass
0037 % at identifying RFI
0038 d_src = boxRejection(d_src, 3, 90);
0039 % plot the data before and after the flagging
0040 figure(1)
0041 t_min   = (d_src.antenna0.receiver.utc - d_src.antenna0.receiver.utc(1))*24*60;
0042 setwinsize(gcf, 575, 750);
0043 subplot(2,1,1)
0044 a = d_src.antenna0.receiver.data(:,1);
0045 plot(t_min, a, 'r');
0046 hold on
0047 a(d_src.flags.rfiFlags) = nan;
0048 plot(t_min, a, 'b');
0049 hold off
0050 xlabel('Time (minutes)');
0051 ylabel('Intensity (backend counts)');
0052 title('Total Intensity 1');
0053 legend('Flagged', 'Un-Flagged');
0054 
0055 subplot(2,1,2)
0056 a = d_src.antenna0.receiver.data(:,6);
0057 plot(t_min, a, 'r');
0058 hold on
0059 a(d_src.flags.rfiFlags) = nan;
0060 plot(t_min, a, 'b');
0061 hold off
0062 xlabel('Time (minutes)');
0063 ylabel('Intensity (backend counts)');
0064 title('Total Intensity 2');
0065 legend('Flagged', 'Un-Flagged');
0066 
0067 % remove the first 1.5 hours, as it looks pretty bad.
0068 d_src.flags.rfiFlags(t_min<90) = 1;
0069 
0070 
0071 % next we integrate the channels, and make sure the noise is decreasing.
0072 
0073 % unflagged time on source
0074 data = d_src.antenna0.receiver.data(:,[1 6]);
0075 data(d_src.flags.rfiFlags,:) = [];
0076 % each data point is 0.01 seconds long
0077 
0078 numIntervals = floor(size(data,1)/(intervalMinutes*60/0.01));
0079 rmsVals = [];
0080 for m=1:numIntervals
0081   if(m==numIntervals)
0082     rmsVals(m,:) = sqrt(var(data));
0083   else
0084     rmsVals(m,:) = sqrt(var(data(round(1:m*intervalMinutes*60/0.01),:)));
0085   end
0086 end
0087 
0088 t_min = intervalMinutes:intervalMinutes:(size(data,1)/(intervalMinutes*60/0.01))
0089     
0090 
0091 
0092 figure(2)
0093 setwinsize(gcf, 575, 500);
0094 plot(t_min, rmsVals)
0095 xlabel('Integration time (minutes)');
0096 ylabel('Rms (backend units)');
0097 legend('Channel 1', 'Channel 2');
0098 title('Integration on fixed source');
0099 
0100 return;