%+=====================================================================+
%
%     Programmers:  Carlos Roithmayr          			Feb 1997
%
%		    NASA Langley Research Center
%		    Spacecraft and Sensors Branch (CBC)
%		    757 864 6778
%		    c.m.roithmayr@larc.nasa.gov
%
%+---------------------------------------------------------------------+
%
%     Purpose:
%
%     Calculate values of the geomagnetic field at 12 points spaced
%     equally on a circle inclined 51.6 deg to Earth's equator, and
%     400 km above Earth's surface.  The field is calculated with
%     IGRF coefficients up to degree and order 10, for the year 1995.00.
%
%     The results reported below list position vector "repe" (km)
%     from Earth's center E* to a point P, expressed in a basis fixed
%     in the Earth: unit vectors e1 and e2 lie in the equatorial plane
%     with e1 in the plane containing the prime meridian, and e3 in the
%     direction of the north pole.
%     The magnetic field vector, "bepe" (Tesla), is also projected into
%     the e1-e2-e3 basis.
%+---------------------------------------------------------------------+

function bepe=get_XYZ_field(alt,t)

global R_mean

R_mean = 6371.2;              % Mean radius for International Geomagnetic
                              % Reference Field (6371.2 km)

[G,H] = IGRF95;               % IGRF coefficients for 1995

nmax = 10;                    % max degree of geopotential
mmax = 10;                    % max order  of geopotential

Kschmidt = schmidt(nmax,mmax);

R_E = 6378.139;               % radius of Earth, km
R_km = R_E + alt;         % radius of circular orbit, km

i_rad=90*pi/180;            % inclination of orbit plane
  Si=sin(i_rad);
  Ci=cos(i_rad);
arg_lat=pi/2-t/(48*60)*pi  ;        % values of argument of latitude starting at the equator, at an input time t seconds

k = 1;

St=sin(arg_lat(k));
Ct=cos(arg_lat(k));

% direction cosine matrix from E to LVLH,
% assuming longitude of ascending node = 0

D(1,1)=Ct;     D(1,2)=-St;     D(1,3)=0;
D(2,1)=Ci*St;  D(2,2)=Ci*Ct;   D(2,3)=-Si;
D(3,1)=Si*St;  D(3,2)=Si*Ct;   D(3,3)=Ci;

% position vector from E* to P, E-basis

repe = R_km*D(:,1)';
repe_array2(k,:) = repe;

[A,ctilde,stilde] = recursion(repe,nmax,mmax);
bepe(k,:) = bfield(repe,nmax,mmax,Kschmidt,A,ctilde,stilde,G,H);


format long
repe_array2;
bepe;

end