function out = get_LSDnSF_Cl36(in,datadir)

% Interpolates into precalculated scaling factor grids to obtain LSDn
% scaling factors for specified Rc, pressure, and solar parameter vectors. 
%
% This code ingests an input structure with the following information:
%   cell array sfdata.t contains time vectors
%   cell array sfdata.dipRc contains Rc at corresponding times
%   cell array sfdata.S contains solar parameter
%   cell array nuclide is nuclide code
%   vector pressure is atmospheric pressure
%
% Adds a cell array containing corresponding vectors of scaling factor
% and returns the structure. 
%
% This needs to be modified to include correct scaling grids for K, Ca,
% etc. and so on. 

% Load grid file

load ([datadir 'LSDn2015_grid']);

% Load data files. 

load([datadir 'LSDn2015_N36K']); 
load([datadir 'LSDn2015_N36Ca']); 
load([datadir 'LSDn2015_N36Ti']); 
load([datadir 'LSDn2015_N36Fe']);
load([datadir 'LSDn2015_N36eth']);

% Now do interpolation.    

for a = 1:length(in.tmax)

    % Make pressure vector. This is planning ahead for time-dependent
    % pressure changes. 
    if iscell(in.pressure)
        % Case in.pressure is cell array with time-dependent p
        temp_p = in.pressure{a};
    else
        % Case in.pressure is vector with non-time-dependent p
        % In this case is scalar for each sample, must be expanded to
        % vector that matches time mesh
        temp_p = zeros(size(in.tmax{a})) + in.pressure(a);
    end
    
     % Catch excessive Rc.
    bigRc = find(in.Rc{a} > 20);
    if ~isempty(bigRc)
        in.Rc{a}(bigRc) = 20 + zeros(size(bigRc));
    end
    
    % Do interpolation...with linear approximation for solar variability
    % assign vector to cell in field in.LSDn
    in.Ksp_LSDn{a} = interp2(sfgrid.Rc,sfgrid.p,sfN36K.b,in.Rc{a},temp_p) + in.S{a}.*interp2(sfgrid.Rc,sfgrid.p,sfN36K.m,in.Rc{a},temp_p);
    in.Casp_LSDn{a} = interp2(sfgrid.Rc,sfgrid.p,sfN36Ca.b,in.Rc{a},temp_p) + in.S{a}.*interp2(sfgrid.Rc,sfgrid.p,sfN36Ca.m,in.Rc{a},temp_p);
    in.Fesp_LSDn{a} = interp2(sfgrid.Rc,sfgrid.p,sfN36Fe.b,in.Rc{a},temp_p) + in.S{a}.*interp2(sfgrid.Rc,sfgrid.p,sfN36Fe.m,in.Rc{a},temp_p);
    in.Tisp_LSDn{a} = interp2(sfgrid.Rc,sfgrid.p,sfN36Ti.b,in.Rc{a},temp_p) + in.S{a}.*interp2(sfgrid.Rc,sfgrid.p,sfN36Ti.m,in.Rc{a},temp_p);
    in.Pf0sp_LSDn{a} = interp2(sfgrid.Rc,sfgrid.p,sfN36eth.b,in.Rc{a},temp_p) + in.S{a}.*interp2(sfgrid.Rc,sfgrid.p,sfN36eth.m,in.Rc{a},temp_p);
    
end

% Return

out = in;
    
    
    
    
    
    
    
    