%% intfl.m calculates the value of the inputted parameter at one point
%% for each inputted field lines using interpolation. The specific
%% format that is required is that each row corresponds to a seperate
%% field line (m field lines) each with n points. The user needs to
%% input the distance along each field line to calculate the inputted
%% parameter (which can be different for field line) which has size
%% m by 1. Linear interpolation is used, but if the switch 'intp' is 
%% set to 0, then the natural log of the inputted parameter is taken
%% before interpolation and then the exp taken to return the correct
%% value. 
%% NOTE: for points at greater distance along a field line than the maximum
%%       distance in the lookup table, then the maximum distance value is
%%       returned.

function vout=intfl(dflm,vmat,dflmp,intp)

% dflm = distance along field line matrix (size m x n).
% vmat = parameter values to be interpolated (size m x n). 
% dflmp = distance along each field line at which to conduct interpolation (size m x 1).
% intp = switch value (0 = take natural log) (scalar).
% vout = interpolated value for each field line (size m x n)

sv=ones(1,size(dflm,2));
poi=sum((dflmp*sv)>dflm,2);                                 %index terms
poic=poi==0;                                                %check for zero terms
poi(poic)=1;                                                %reset zero terms (will be corrected below)
if intp==0
 vmat=log(vmat);                                            %take log of input parameters if required
end 
ly=length(dflmp);                                           %number of input terms
indl=(poi-1)*ly+(1:ly)';                                    %lower matrix index terms
indu=indl+ly;                                               %upper matrix index terms
inmx=prod(size(dflm));                                      %maximum index size
poim=indu>inmx;                                             %upper terms greater than maximum
indl(poim)=indl(poim)-ly;                                   %reduce poic lower index terms by one (will correct below)
indu(poim)=indu(poim)-ly;                                   %reduce poic upper index terms by one (will correct below)
grad=(vmat(indu)-vmat(indl))./(dflm(indu)-dflm(indl));      %gradient term
vout=vmat(indl)+grad.*(dflmp-dflm(indl));                   %required parameter values at points of interest
vout(poic)=vmat(poic,1);                                    %reset values for points at base of field line
vout(poim)=vmat(poim,size(dflm,2));                         %reset terms above maximum dflm to maximum dflm vmat value
if intp==0
 vout=exp(vout);                                            %take exp of output values if required
end 
