time   calls  line
                  1 function [A,G,tsel,asel,gsel,Afilt,Gfilt] = interpolateAlpha_DD(t,varargin)
                  2 % function [A,G] = interpolateAlpha_DD(t,[N, M])
                  3 %
                  4 % This function reads in the database of alpha values, and interpolates the
                  5 % values onto the times given by t. It performs two steps: the first is to
                  6 % smooth the alpha values over the range considered, and the second in to
                  7 % interpolate linearly between the smoothed alpha values.
                  8 %
                  9 % The smoothing is accomplished by applying a Savitsky-Golay filter. This
                 10 % is a generalization of a moving window average. At each point, fit a
                 11 % polynomial of order N to the M surrounding data points. The smoothed data
                 12 % point is then the value of the fitted polynomial.
                 13 %
                 14 % Inputs:
                 15 %   t - times in MJD to return alpha and gain values for
                 16 %   N - (optional) order of the polynomial to use in the smoothing.
                 17 %       2-element vector. First entry is the total intensity diode 
                 18 %       channels (columns 1,2,7,8)
                 19 %       Second entry is the polarisation diode channels (columns 3,4,5,6)
                 20 %       N = 0: moving window average
                 21 %       N = 1: linear fit
                 22 %       N = 2: quadratic fit
                 23 %       etc.
                 24 %   M - (optional) number of noise diode events to fit over. 2-element
                 25 %   vector like the N.
                 26 % If you don't specify N and M, then the routine uses a moving window
                 27 % average (N=0) and works out the window widths necessary to get 
                 28 % smoothing windows 30 minutes wide in total intensity and 10 minutes
                 29 % wide in polarization channels.
                 30 %
                 31 % Outputs:
                 32 %   A - 5 column matrix of alpha values corresponding to times t
                 33 %   G - 5 column matrix of gain values corresponding to times t
                 34 %
                 35 % ogk - December 2010
                 36 %
                 37 
             1   38 if ~(size(varargin,2) == 0 || size(varargin,2) == 2) 
                 39     error('AlphaInterpolation:WrongNumberOfArgs',...
                 40         'Incorrect number of arguments. Correct usage: [A,G] = interpolateAlpha(t,[N,M])')
                 41 end
                 42 
                 43 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 44 %  Read in the database of alpha values
                 45 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  2.67       1   46 [times,alpha,gain,r,TsysTND,temps,horz,equa] = readAlphaDatabase_DD(); 
                 47 
                 48 
             1   49 if ( t(1) > times(end) || t(end) > times(end)) 
                 50   display('-------WARNING WARNING WARNING -------');
                 51   display('Alpha database not up to date');
                 52   display('Using last valid value from database');
                 53   display('Re-run dataset at a future date when database is complete');
                 54 
                 55   asel = alpha(end,:);
                 56   tsel = times(end,:);
                 57   gsel = gain(end,:);
                 58   
                 59   Afilt = zeros(length(tsel),5);
                 60   Gfilt = zeros(length(tsel),5);
                 61   
                 62   A = repmat(asel, [length(t) 1]);
                 63   G = repmat(gsel, [length(t) 1]);
                 64   
                 65   return;
                 66 end  
                 67   
                 68 
                 69 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 70 %  Identify the range of alpha values to use in the interpolation
                 71 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 72 
                 73 % Find the alpha value with lie within the time window
             1   74 I = find(times-t(1) >= 0); 
             1   75 Ifirst = I(1) - 1; 
  0.01       1   76 I = find(times-t(end) >= 0); 
             1   77 Ilast = I(1) + 1; 
                 78 
             1   79 Irange = Ifirst:Ilast; 
             1   80 tsel = times(Irange); 
             1   81 asel = alpha(Irange,:); 
             1   82 gsel = gain(Irange,:); 
                 83 
                 84 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 85 %  Set up the filtering parameters
                 86 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 87 % Either use the supplied polynomial order and window width, or determine
                 88 % it from the cadence of the noise events in the data
             1   89 if size(varargin,2) == 2 
                 90     NI = varargin{1}(1);
                 91     NP = varargin{1}(2);
                 92     N = [NI NI NP NP NP NP NI NI];
                 93     MI = varargin{2}(1);
                 94     MP = varargin{2}(2);
                 95     M = [MI MI MP MP MP MP MI MI];
             1   96 else 
                 97     % The smoothing window width must be 20 minutes wide for the 1 and 5
                 98     % vectors, and 10 minutes wide for the 2,3,4 vectors, of alpha and
                 99     % gain
                100     
             1  101     Nevents = Ilast-Ifirst+1; % number of events 
             1  102     duration = (t(end)-t(1))*24*60; % number of minutes spanned by the data 
             1  103     freq = Nevents/duration; % number of noise diode events per minute 
                104     
                105     % Therefore, to get the targeted window width:
             1  106     N = zeros(1,5); % moving window average 
             1  107     Ml = round(freq*20); 
             1  108     Ms = round(freq*10); 
             1  109     M = [Ml Ml Ms Ms Ms Ms Ml Ml]; 
                110 %     N = [2 2 3 3 3 3 2 2]; % the order of the polynomial to fit
                111 %     M = [10 10 3 3 3 3 10 10]; % the number of noise diode events to fit to
             1  112 end 
                113 
                114 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                115 %  Filter the alpha and gain values
                116 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                117 % Matrices to hold the filtered data points:
             1  118 Afilt = zeros(length(tsel),5); 
             1  119 Gfilt = zeros(length(tsel),5); 
             1  120 for k=1:length(tsel) 
                121     % Fit an N'th degree polynomial to a window of M points around the
                122     % current noise diode event
  0.10      24  123     I = find(times == tsel(k)); % index of position in times vector corresonding to the current point 
                124     
                125     % Perform the fit for the alpha and gain values
            24  126     for m=1:5 
                127         % Select the window, and check to make sure that none of the data
                128         % points are too far away in time
           120  129         fitI = (I-floor(M(m)/2)):(I+floor(M(m)/2)); % the samples to fit over 
           120  130         x = times(fitI); 
                131         % If the points are more than 2 hours away, then don't use them:
           120  132         fitI = fitI(abs((x-tsel(k))*24)<=2); 
           120  133         x = x(abs((x-tsel(k))*24)<=2); 
                134         
                135         % Fit a polynomial to the t-centred data, then assign the value of
                136         % the polynomial at t=0 to the matrices of filtered values
  0.11     120  137         Pa = polyfit(x-tsel(k),alpha(fitI,m),N(m)); 
  0.11     120  138         Pg = polyfit(x-tsel(k),gain(fitI,m),N(m)); 
           120  139         Afilt(k,m) = Pa(end); 
           120  140         Gfilt(k,m) = Pg(end); 
           120  141     end 
                142     
            24  143 end 
                144 
                145 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                146 %  Interpolate linearly between the filtered values
                147 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  0.40       1  148 A = interp1(tsel,Afilt,t,'linear'); 
  0.40       1  149 G = interp1(tsel,Gfilt,t,'linear'); 
                150 
                151 
             1  152 return;