%% dis2th.m determines the magnetic latitude of a of a point which is a
%% certain distance along a particular L-shell. This is based on
%% equation A5 given in Angerami and Thomas.

function latm=dis2th(dflm,latmb,r0)

% dflm = distance along field line of point of interest from base (km)
% latmb = magnetic latitude at base of field line (rads)
% r0 = radial distance of base (km) (scalar)
% latm = magnetic latitude of point of interest (rads) 

sz=size(dflm);                %size of input
poi=latmb<0;                  %terms with negative latitudes
latmb=abs(latmb);             %take absolute value of latitude
dflm=dflm(:);
latmb=latmb(:);

%Determine the constant terms from A4
c1=r0./(2*sqrt(3)*cos(latmb).^2);
y0=sqrt(3)*sin(latmb);
c2=log(y0+sqrt(1+y0.^2))+y0.*sqrt(1+y0.^2);
f=-(dflm./c1)+c2;

%Calculate lookup vector interpolate from f to magnetic latitude
d80=deg2rad(80);
d81=deg2rad(81);
latv=[linspace(0,d80,81) linspace(d81,pi/2,50)];       %latitude lookup vector
y=sqrt(3)*sin(latv);                                   %y value from A2
fv=log(y+sqrt(1+y.^2))+y.*sqrt(1+y.^2);                %f values at each latitude

%Check f values are not outside fv range
f(f<0)=0;
fmx=max(fv);
f(f>fmx)=fmx;

nop=length(dflm);                            %number of points of interest
nopd=100;                                    %nop divider
latm=zeros(nop,1);                           %magnetic latitude of points of interest
imx=ceil(nop/nopd);                          %number of required loops
for i=1:imx
 li=1+nopd*(i-1);                            %lower index
 if i~=imx
  ui=nopd*i;                                 %upper index
 else
  ui=nop;
 end 
 latmi=interp1(fv,latv,f(li:ui),'cubic');    %i'th interpolation lookup 
 latm(li:ui)=latmi;                          %record i'th results
end  

latm=reshape(latm,sz);                       %reshape to original size
latm(latm<1e-6)=0;                           %assume values less than 1e-6 are 0 (round off error)
latm(poi)=-latm(poi);                        %reset negative latitudes to negative values

