interpolateAlpha

 function [Ai,Gi] = interpolateAlpha(utc,t,A,G,selection)
   Apply the passed alpha values
 or
 function [Ai,Gi] = interpolateAlpha(utc,selection)
   Use the database values

0001 function [Ai,Gi,tsel,asel,gsel,Afilt,Gfilt] = interpolateAlpha(utc,varargin)
0002 % function [Ai,Gi] = interpolateAlpha(utc,t,A,G,selection)
0003 %   Apply the passed alpha values
0004 % or
0005 % function [Ai,Gi] = interpolateAlpha(utc,selection)
0006 %   Use the database values
0007 %
0008 
0009 if nargin == 2
0010     selection = varargin{1};
0011     mode = 'DATABASE';
0012 elseif nargin == 5
0013     selection = varargin{4};
0014     mode = 'CURRENT';
0015 else
0016     disp('interpolateAlpha:: Invalid number of arguments passed to interpolateAlpha!')
0017     return
0018 end
0019 
0020 if strcmp(mode,'CURRENT')
0021     % use the alpha values passed by the function call:
0022     t = varargin{1};
0023     A = varargin{2};
0024     G = varargin{3};
0025     % the number of columns in A and G are:
0026     NA = size(A,2); 
0027     NG = size(G,2);
0028     % Matrices of interpolated alpha and gain values:
0029     Ai = zeros(length(utc),NA);
0030     Gi = zeros(length(utc),NG);
0031     % The order of the polynomial to fit to the alpha and gain terms
0032     nOrder = 1;
0033     for k=1:NA
0034         Ai(:,k) = polyval(  polyfit(t-t(1),A(:,k),nOrder)  ,utc-t(1));
0035     end
0036     for k=1:NG   
0037         Gi(:,k) = polyval(  polyfit(t-t(1),G(:,k),nOrder)  ,utc-t(1));
0038     end
0039     tsel = zeros(size(t));
0040     asel = zeros(size(t));
0041     gsel = zeros(size(t));
0042     Afilt = zeros(size(t));
0043     Gfilt = zeros(size(t));
0044 else
0045     % use the values contained in the alpha database:
0046     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0047     %  Read in the database of alpha values
0048     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0049     
0050     % Read in the data for the relevant time. If it is within 2 hours of a
0051     % month change then read in the adjacent month data too.
0052     x = readAlphaDatabase([utc(1)-2/24;utc;utc(end)+2/24]);
0053     times = x{1};
0054     Iflagged = x{17};
0055     if strcmp(selection,'FILTERED')
0056         A = x{9};
0057         G = x{8};
0058     elseif strcmp(selection,'POLONLY')
0059         GPO = x{6};
0060         GDD = x{3};
0061         G = [GPO(:,1) GDD(:,2:4) GPO(:,2)];
0062         A = x{2}(:,2:4);
0063     elseif strcmp(selection,'CLASSIC')
0064         A = x{2};
0065         G = x{3};
0066     end
0067     
0068     % Get rid of flagged values:
0069     times = times(~Iflagged);
0070     A = A(~Iflagged,:);
0071     G = G(~Iflagged,:);
0072 
0073     if ( utc(1) > times(end) || utc(end) > times(end))
0074         display('interpolateAlpha:: -------WARNING WARNING WARNING -------');
0075         display('interpolateAlpha:: Alpha database not up to date');
0076         display('interpolateAlpha:: Using last valid value from database');
0077         display('interpolateAlpha:: Re-run dataset at a future date when database is complete');
0078 
0079         asel = A(end,:);
0080         tsel = times(end,:);
0081         gsel = G(end,:);
0082 
0083         Afilt = zeros(length(tsel),size(A,2));
0084         Gfilt = zeros(length(tsel),size(G,2));
0085 
0086         Ai = repmat(asel, [length(utc) 1]);
0087         Gi = repmat(gsel, [length(utc) 1]);
0088 
0089         return
0090     end  
0091 
0092     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0093     %  Identify the range of alpha values to use in the interpolation
0094     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0095 
0096     % Find the alpha value with lie within the time window
0097     I = find(times-utc(1) >= 0);
0098     Ifirst = I(1) - 1;
0099     if(Ifirst == 0)
0100       Ifirst = 1;
0101     end
0102     I = find(times-utc(end) >= 0);
0103     Ilast = I(1) + 1;
0104     if(Ilast > length(times))
0105       Ilast = length(times);
0106     end
0107     
0108     Irange = Ifirst:Ilast;
0109     tsel = times(Irange);
0110     asel = A(Irange,:);
0111     gsel = G(Irange,:);
0112 
0113     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0114     %  Set up the filtering parameters
0115     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0116     % Determine the window width and polynomial order it from the cadence
0117     % of the noise events in the data
0118     
0119     % The smoothing window width must be 20 minutes wide for the total
0120     % intensity channels, and 10 minutes wide for the polarization vectors
0121 
0122     Nevents = Ilast-Ifirst+1; % number of events
0123     duration = (utc(end)-utc(1))*24*60; % number of minutes spanned by the data
0124     freq = Nevents/duration; % number of noise diode events per minute
0125 
0126     % Therefore, to get the targeted window width:
0127     NA = zeros(1,size(A,2)); % apply a moving window average to the data
0128     NG = zeros(1,size(G,2));
0129     Ml = round(freq*20);
0130     Ms = round(freq*10);
0131 
0132     MG = [Ml Ms Ms Ms Ml];
0133     if strcmp(selection,'POLONLY')
0134         % PolOnly A vector has [P1 P2 P3]
0135         MA = [Ms Ms Ms];
0136     elseif strcmp(selection,'FILTERED')
0137         % Filtered A vector has [P1 P2 P3]
0138         MA = [Ms Ms Ms];
0139     elseif strcmp(selection,'CLASSIC')
0140         % Classic A vector has [I1 P1 P2 P3 I2]
0141         MA = [Ml Ms Ms Ms Ml];
0142     end
0143 
0144     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0145     %  Filter the alpha and gain values
0146     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0147     % Matrices to hold the filtered data points:
0148     Afilt = filterVector(times,tsel,A,MA,NA);
0149     Gfilt = filterVector(times,tsel,G,MG,NG);
0150     
0151     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0152     %  Interpolate linearly between the filtered values
0153     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0154     if(length(tsel)==2)
0155       AA = mean(Afilt);
0156       GG = mean(Gfilt);
0157       Ai = repmat(AA, [length(utc) 1]);
0158       Gi = repmat(GG, [length(utc) 1]);
0159     else
0160       Ai = interp1(tsel,Afilt,utc,'linear');
0161       Gi = interp1(tsel,Gfilt,utc,'linear');
0162     end
0163   end
0164 
0165 end
0166 
0167 function Vfilt = filterVector(times,tsel,V,M,N)
0168 % this function filters the data in vector V. For each sample k in column
0169 % m, fit a polynomial of order N(m) to the M(m) samples around sample k.
0170 % Replace the sample by the fitted polynomial at the sample time.
0171 
0172 Vfilt = zeros(length(tsel),size(V,2));
0173 for k=1:length(tsel)
0174     % Fit an N'th degree polynomial to a window of M points around the
0175     % current noise diode event
0176     I = find(times == tsel(k)); % index of position in times vector corresonding to the current point
0177 
0178     % Perform the filtering for each column of sample k
0179     for m=1:size(V,2)
0180         % Select the window, and check to make sure that none of the data
0181         % points are too far away in time
0182         fitI = (I-floor(M(m)/2)):(I+floor(M(m)/2)); % the samples to fit over
0183         % make sure that the fit indices do not exceed the actual times
0184         % available
0185         fitI = fitI(fitI <= length(times));
0186         fitI = fitI(fitI > 0); % added to prevent negative indexing MI 5/5/13
0187         x = times(fitI);
0188         % If the points are more than 2 hours away, then don't use them:
0189         fitI = fitI(abs((x-tsel(k))*24)<=2);
0190         x = x(abs((x-tsel(k))*24)<=2);
0191 
0192         % Fit a polynomial to the t-centred data, then assign the value of
0193         % the polynomial at t=0 to the matrix of filtered values
0194         P = polyfit(x-tsel(k),V(fitI,m),N(m));
0195         Vfilt(k,m) = P(end);
0196     end
0197 end
0198 
0199 end