%% Similar to flpar.m, flpar3.m determines the arrays describing the
%% co-ordinates of the given field lines in the polar regions. It is
%% different in that it requires the user to input a single maximum
%% radial distance as an L-shell. This specifies the upper bound for
%% all the field lines in the context that all field lines will have 
%% a length equal to the length of the specificed field line. If they
%% are on a lower L-shell, then they have the maximum for their L-shell.
%%

function [lonmat,latmat,latmatm,dismat,dflm]=flpar3(latmb,lonlsh,r0,lshrt,nop)

% latmb = base magnetic latitude of field lines (rads)
% lonlsh = magnetic longitude of field lines (rads)
% r0 = base radial distance (km) (scalar)
% lshrt = upper radial L-shell (scalar)
% nop = number of points per field line (scalar)
%
% lonmat = geographic longitude of points of interest (rads)
% latmat = geographic latitude of points of interest (rads)
% latmatm = magnetic latitude of points of interest (rads)
% dismat = distance from center of Earth of points of interest (km)
% dflm = distance along the field line from base (km)

%convert input to column vectors
latmb=latmb(:);
lonlsh=lonlsh(:);

sv1=ones(length(latmb),1);
sv2=ones(1,nop);
svm=sv1*sv2;

latmxb=lshpar(lshrt,-1,r0);                       %magnetic latitude of upper radial L-shell
dflmx=distfleq(latmxb,lshrt);                     %length of upper radial L-shell
lsh=lshpar(-sv1,latmb,r0*sv1);                    %L-shell of field lines
dflt=distfleq(latmb,lsh);                         %length of the field line from the equator to r0
dflt(dflt>dflmx)=dflmx;                           %reset field lines with lengths greater than the upper radial value to the radial values length
dv=logspace(-2,0,(nop-1));                        %fractional vector to break field line up into
dv=[dv-dv(1) 1];                                  %final version of corrected field line fractional vector
dflm=dflt*dv;                                     %required distant along field line of points of interest
latmatm=dis2th(dflm,latmb*sv2,r0);                %magnetic latitude of points of interest below 80 degrees
dismat=lshpar(lsh*sv2,latmatm,-svm);              %radial displacement of points of interest
dismat(:,1)=r0;                                   %reset base value to r0 exactly
poi=abs(latmb)>deg2rad(88.5);                     %check for terms above 88.5 degree
if any(poi)
 dismat(poi,:)=r0+ones(sum(poi),1)*dflmx*dv;      %replace dismat values at 90 degrees with base radial value plus dflm values
end 
[lonmat,latmat]=mag2geo3(lonlsh*sv2,latmatm);     %geographic co-ordinates of points of interest

