function [a,phi,t,fs] = toylineamp(varargin)
% LINEAMP
% [a,phi,t,ff,fs] = lineamp(Channel, T0, f, NSec, NSeg)
% Channel string name of channel to process 
% T0        Start time (GPS secs)
% f         line frequency in Hz. May be a vector.
% NSec      number seconds per segment
% NSeg      number of segments
%
% a 	line amplitude. Rows correspond to line frequencies. 
% phi	line phase. Rows correspond to line frequencies.
% t     time stamps for the beginning of the averaging segment
% fs    sample rate
%
% Example: [a,phi,t,fs]=toylineamp('H1:LSC-AS_Q',754562512,[41.25 151.3,973.3],15,40);
%
% 

if nargin ~= 5
  error('[a,phi,t,fs] = toylineamp(Channel, T0, f, NSec, NSeg)');
end

Channel = varargin{1};
T0 = varargin{2};
f = varargin{3};
NSec = varargin{4};
NSeg = varargin{5};


% pre-allocate
a = zeros([length(f),NSeg]);
phi = a;
t = T0 + (0:(NSeg-1))*NSec;

% Set up chanvector for filelist operation (not necessary with LDR queries)
[gpsTimes,frameFiles,frameDurs]=dir2framelist('./samrds/S3/L3/LHO/H-RDS_R_L3-7545');

% Create Channel Structure
dataStruct = chanstruct(Channel);

f = f(:);
for j = 1:NSeg
  
  % get data
  GPS = t(j);
  [data,fs] = chanvector(dataStruct, GPS, NSec, gpsTimes, frameFiles, frameDurs);
  
  % Create window function
  window = hann(length(data))';
  
  % local oscillator has zero phase at GPS t=0
  if (j == 1)
    k0 = 0:(fs*NSec-1);
    k0 = f(:)*k0/fs;
    k0 = exp(2*pi*i*k0);
  end
  for l = 1:length(f)
    k = k0(l,:)*exp(2*pi*i*GPS*f(l));
    data = data.*window;
    % Create phasor
    z(l,j) = sum(k.*data);
  end

end

z=z*2/(fs*NSec); % Scale integral

% Create amplitude and phase
a = abs(z);
phi = angle(z);

% Plot Results

  xl = sprintf('Time (s) from GPS %d',T0);
      
  tt = zeros([length(f),length(t)]);
  lgnd = cell([1,length(f)]);
  for k = 1:length(f)
    tt(k,:) = t;
    lgnd{k} = sprintf('%8.2f Hz',f(k));
  end
  tt = tt - T0 + NSec/2;
  
  subplot(2,1,1);
  semilogy(tt',a');
  grid;
  xlabel(xl);
  ylabel('Amplitude (counts)');
  legend(lgnd);

  subplot(2,1,2);
  plot(tt',phi');
  xlabel(xl);
  ylabel('Phase (radians)');
  legend(lgnd);
  title(sprintf('%s %ds from %d',Channel(1:3),NSec*NSeg,T0));
      
  print(gcf,'-depsc2',sprintf('Line-%s-%d-%d',Channel(1:2),T0,NSec));


ff = f;

return