%% distfl2.m calculates the distance a point is along the given 
%% field line. Based on distfl.m, it uses L-shell rather than th0.
%%

function dflm=distfl2(latmatm,lsh,r0)

% latmatm = magnetic latitude of point of interest
% lsh = L-shell point of interest on
% r0 = base radial height of field line
% dflm = distance from base (height r0) to point of interest

th0=lshpar(lsh,-ones(size(lsh)),r0);
th=abs(latmatm);
th0=abs(th0);
y=sqrt(3)*sin(th);
y0=sqrt(3)*sin(th0);
s1=log(y+sqrt(1+y.^2))+y.*sqrt(1+y.^2);
sm=log(y0+sqrt(1+y0.^2))+y0.*sqrt(1+y0.^2);

dflm=(r0./(2*sqrt(3)*cos(th0).^2)).*(sm-s1);


