function [f1save,s1save] = srftables(clist)

%% copied from /asl/matlib/sconv/sconv2.m
%% input : clist = select 1 -- 2378
%% output : cells of wavenumber and srf (at 0.0025 cm-1 resolution)
%%
%% example
%% [f1save,s1save] = srftables([1 1400 2378]);
%% plot(f1save{1},s1save{1})
%% plot(f1save{2},s1save{2})
%% plot(f1save{3},s1save{3})

% disp('WARNING off, MATLAB:imagesci:hdf:removalWarningHDFSD')
warning('off', 'MATLAB:imagesci:hdf:removalWarningHDFSD')

fin = (1:890000)-1;
fin = 605 + 0.0025 * fin;

rin = ones(size(fin'));

sfile = '/asl/matlab2012/srftest/srftables_m140f_withfake_mar08.hdf';

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

% guarantee input vectors are columns
fin = fin(:);
clist = clist(:);

% check rin dimensions
[m,n] = size(rin);
if m ~= length(fin)
  error('number of rows of rin do not match length of fin')
end

% default to zero channel offsets
coff = zeros(length(clist), 1);

% read the SRF data
[alist, fattr] = h4sdread(sfile);

for i = 1 : length(alist)
  switch alist{i}{1}
    case 'chanid', chanid = double(alist{i}{2})';
    case 'freq',   freq   = double(alist{i}{2})';
    case 'fwgrid', fwgrid = double(alist{i}{2})';
    case 'srfval', srfval = double(alist{i}{2})';
    case 'width',  width  = double(alist{i}{2})';
  end
end
clear alist fattr

% use the subset of channels specified in clist
cind = interp1(chanid, 1:length(chanid), clist, 'nearest');
chanid = chanid(cind);
freq = freq(cind) + coff;
srfval = srfval(cind,:);
width = width(cind);
[nchan, nspts] = size(srfval);

% srfreq is an nchan x nspts array of frequency points for srfval
srfreq = (width * fwgrid) + freq * ones(1, nspts);

% set the output frequency list
fout = freq;

% check that fin is uniformly spaced and increasing
dfin = diff(fin);
dmax = max(dfin);
dmin = min(dfin);
if abs(dmax - dmin) > 10e-9
  error('fin must be uniformly spaced')
end

% get the step size of fin
dv = dfin(1);
if dv <= 0
  error('fin must be in increasing order')
end

% get the span of fin
nfin = length(fin);
f1 = fin(1);
f2 = fin(nfin);

% initialize output array
[m,n] = size(rin);
rout = zeros(nchan, n);

% loop on SRFs in cind
for j = 1 : nchan

  % get the frequency span of the current SRF
  v1 = srfreq(j, 1);
  v2 = srfreq(j, nspts);

  % if the SRF is outside fin, skip this channel
  if v2 <= f1 | f2 <= v1
    fprintf(1, 'sconv2(): WARNING -- SRF %d outside of input range\n', ...
	    chanid(j));
    continue
  end

  % if the SRF overlaps fin, lop it off to fit; 
  % giv a warning message if we lop more than dv
  if v1 < f1
    if v1 < f1 - dv
      fprintf(1, 'sconv2(): WARNING -- truncating LHS of SRF %d\n', ...
	      chanid(j));
    end
    v1 = f1;
  end
  if f2 < v2
    if f2 + dv < v2
      fprintf(1, 'sconv2(): WARNING -- truncating RHS of SRF %d\n', ...
	      chanid(j));
    end
    v2 = f2;
  end

  % find the indices of the current SRF in fin 
  v1ind = ceil((v1-f1)/dv) + 1;
  v2ind = floor((v2-f1)/dv) + 1;

  % interpolate the SRF to a subinterval of the fin grid
% s1 = interp1(srfreq(j,:), srfval(j,:), fin(v1ind:v2ind), 'spline');
  s1 = interp1(srfreq(j,:), srfval(j,:), fin(v1ind:v2ind), 'linear');

  % normalize the SRF
  s1 = s1 ./ sum(s1);

  s1save{j} = s1;
  f1save{j} = fin(v1ind:v2ind);
  
  % apply the SRF
  rout(j, :) = s1' * rin(v1ind:v2ind, :);
end

