%% flpar2.m (field line parameters) determines the arrays describing the
%% co-ordinates of the given field lines. Based on flpar.m, it is designed to
%% output the required arrays when magnetic local time is used, rather
%% than magnetic longitude.
%% Note that for the calculation of the geographic and magnetic co-ordinates 
%% that the solar declination will be taken as 0 (ie equinox conditions) at
%% all times (ie no annual variation). This at most introduces an error of 
%% the order of at most +-4 degrees (summer and winter solstice) in the 
%% magnetic longitude ie 4/360 ~1% error which is enough to ignore. See
%% timml2lon.m for correct declination values.
%% 
%% One major difference between flpar and flpar2 is that flpar2 also writes to
%% an ascii file a copy of lonb, lonbc, latb and latbc to both the msis90 and
%% iri2001 subdirectories. This removes the need to have the IRI2001 and MSIS90
%% lookup programs continually writing these values to hard drive when they
%% are called. The values saved are those for points lsh<=lshcp and is done 
%% if spp==1.
%%
%% Note that at very large L-shells (say L>35, ~80 degrees)
%% that the results obtained will be in error, the more so
%% the closer to 90 degrees.
%%
%% 31/7/2001 - modified to call IRI2001.
%% 9/3/2002 - modify to call NRLMSIS2000 or MSIS90. 
%%

function [latmatm,dismat,dflm,lonb,lonbc,latb,latbc]=flpar2(poip,lsh,mlt,timmat,r0,nop,spp,sms)

% poip = field lines for which to complete full field line calculations
% lsh = L-shell of field lines 
% mlt = magnetic local time of field lines (rads)
% timmat = UT of field lines
% r0 = base radial distance (km)
% nop = number of points along each field line
% spp = save parameter (spp=1 save, anything else does not save)
% sms = switch to determine which version of MSIS to use (0 = MSIS90, 1 = NRLMSIS2000)(scalar)
%
% 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)
% lonb = r0 geographic longitude of points of interest (rads)
% lonbc = conjugate r0 geographic longitude of points of interest (rads)
% latb = geographic latitude of points of interest (rads)
% latbc = conjugate r0 geographic latitude of points of interest (rads)

%% Calculate L-shell parameters for complete field lines (latmatm, dismat and dflm)
lsh=lsh(:);
sv1=ones(length(lsh),1);
sv2=ones(1,nop);
latmb=lshpar(lsh,-1*sv1,r0*sv1);                       %magnetic latitude of base of field line
dflt=distfleq(latmb,lsh);                              %length of the field line from the equator to r0
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
dismat=lshpar(lsh*sv2,latmatm,-sv1*sv2);               %radial displacement of points of interest
dismat(:,1)=r0;                                        %reset base value to r0 exactly

%% Calculate geographic and geomagnetic longitudes and latitudes at r0, which vary
%% as function of UT. As noted above, assume solar declination is 0 at all times.
mlt=mlt(:);                                            %
timmat=timmat(:)';                                     %make sure timmat is a row vector
nopt=length(timmat);                                   %number of points in timmat
svt=ones(1,nopt);
lonmats=longch(pi-(pi/12)*timmat);                     %geographic longitude of subsolar point 
lonmatms=geo2mag3(lonmats,0*svt);                      %magnetic longitude of subsolar point
lonmatm=(pi/12)*(mlt-12)*svt+sv1*lonmatms;             %magnetic longitude of point of interest     
lonmatm=longch(lonmatm);                               %correct magnetic longitude to -pi -> pi

[lonb,latb]=mag2geo3(lonmatm,latmatm(:,1)*svt);        %geographic co-ordinates of field lines at r0
[lonbc,latbc]=mag2geo3(lonmatm,-latmatm(:,1)*svt);     %conjugate geographic co-ordinates of field lines at r0

%%Write the longitude and latitudes to the ascii files data2.txt if required
if spp==1
 %convert from rads to degrees
 lonb2=rad2deg(lonb);                 
 latb2=rad2deg(latb);                 
 lonb2c=rad2deg(lonbc);                 
 latb2c=rad2deg(latbc);      
 
 lonb21=lonb2(poip,:);
 lonb21c=lonb2c(poip,:);
 latb21=latb2(poip,:);
 latb21c=latb2c(poip,:);

 if sms==0
  cd msis90
 else
  cd msis2000
 end

 %Save parameters to disk
 fid=fopen('data2.txt','w');
 fprintf(fid,'%4.3f %4.3f\r\n',[[lonb21(:)' lonb21c(:)'];[latb21(:)' latb21c(:)']]);
 fclose(fid);

 cd ..
    
 %convert negative longitudes into positive values (add 2pi, required format of IRI2001)
 lonb3=longch2(lonb);
 lonb3c=longch2(lonbc);
 
 %convert from rads to degrees
 lonb3=rad2deg(lonb3);                 
 latb3=rad2deg(latb);                 
 lonb3c=rad2deg(lonb3c);                 
 latb3c=rad2deg(latbc);     
 
 lonb31=lonb3(poip,:);
 lonb31c=lonb3c(poip,:);
 latb31=latb3(poip,:);
 latb31c=latb3c(poip,:);
 
 cd iri2000
 
 %Save parameters to disk
 fid=fopen('data2.txt','w');
 fprintf(fid,'%4.3f %4.3f\r\n',[[lonb31(:)' lonb31c(:)'];[latb31(:)' latb31c(:)']]);
 fclose(fid);
 
 cd ..
 
end 