9.2 C Program File: RF_GetReference.c
/***********************************************************************
*
*
* FILE NAME: RF_GetReference.c
*
*
* PURPOSE: This function determines if a requested reference
* vector is in need of update. It uses analytic routines
* to update vectors and these updates are reflected in the
* reference.h include file.
*
*
* FILE REFERENCES:
*
* Name IO Description
* ------------ -- ---------------------------------
* none
*
*
* EXTERNAL VARIABLES:
*
* Source : debug.h
*
* Name Type IO Description
* ------------- ------- -- -----------------------------
* debug_file_handle FILE* I File handle for debug file
* name
* debug_level int[9] I Debug level array
*
* Source : HD_reference.h
*
* Name Type IO Description
* ------------- ------- -- -----------------------------
* ephem_file_lu long I FORTRAN logical unit number
* for the ephemeris file
* ephem_method char I Method for computing
* ephemeris information:
* F = Use ephemeris file
* A = Compute analytically
* using Keplerian
* elements
* keplerian double[6] I Keplerian orbital elements at
* the epoch time
* (orbital_t_epoch):
* [1] Semimajor axis [km]
* [2] Eccentricity
* [3] Inclination [rad]
* [4] Right ascension of
* the ascending node
* [rad]
* [5] Argument of perigee
* [rad]
* [6] Mean anomaly [rad]
* m_order long I Order of magnetic field
* maxit long I Maximum number of iterations
* to converge the true
* anomaly
* MU_E double I Earth gravitational constant
* [km^3/sec^2]
* NUMPTS int I Number of points used by the
* EPHEMRD interpolator
* orbital_t_epoch double I Base epoch time of the
* orbital elements [sec]
* THREEB double I Gravitational constant of
* perturbations [Km^2]
* ttol double I Tolerance in the calculations
* of the true anomaly [rad]
* t_b_ref double IO Time of last calculated Earth
* magnetic field vector [sec]
* t_e_ref double IO Time of last calculated s/c
* to Earth unit vector [sec]
* t_m_ref double IO Time of last calculated s/c
* to Moon unit vector [sec]
* t_o_ref double IO Time of last calculated orbit
* normal unit vector [sec]
* t_rv_ref double IO Time of last calculated s/c
* position and velocity
* vectors[sec]
* t_s_ref double IO Time of last calculated s/c
* to Sun unit vector [sec]
* e_pos double[3] O S/C to Earth unit vector
* m_pos double[3] O S/C to Moon unit vector
* mag_field double[3] O Earth magnetic field vector
* [mG]
* mag_field_unit double[3] O Earth magnetic field unit
* vector
* orbit_normal double[3] O Orbit normal unit vector
* s_c_pos double[3] O S/C position vector [km]
* s_c_vel double[3] O S/C velocity vector [km/sec]
* s_pos double[3] O S/C to Sun unit vector
*
*
*
* EXTERNAL REFERENCES:
*
* Name Description
* ------------------ ----------------------------------------------
* c_ephemrd Retrieves vectors from an ephemeris file and
* interpolates them for a requested time
* c_calpvs Generates s/c position and velocity vectors
* using J2 effects
* c_sunlunp Generates Earth to Sun or Earth to Moon
* vectors
* c_emagfld Generates Earth magnetic field vectors
* c_nmlist Opens the magnetic field file for reading
* GetSun Compute s/c to Sun unit vector
* GetOrbitNormal Compute orbit normal vector
* GetEarth Compute s/c to Earth vector
* GetMoon Compute s/c to Moon unit vector
* SecsToCalendar Converts time from secornds to standard
* calendar format
* c_packst Converts time from standard calendar format to
* an unpacked array format
* c_calmjd Computes the modified Julian date of an
* unpacked array format time
* c_jgrenha Computes the Greenwich Hour Angle using
* analytical data
* c_unvec3 Unitizes a vector and computes its magnitude
*
*
* ABNORMAL TERMINATION CONDITIONS, ERROR AND WARNING MESSAGES:
* none
*
*
* ASSUMPTIONS, CONSTRAINTS, RESTRICTIONS: none
*
*
* NOTES:
* CALLED BY: InitReference, CalcNadirAngle, ConvertAttitude,
* ComputeAttitude, CompSunNad, CalcLambdaPhi
*
*
* REQUIREMENTS/FUNCTIONAL SPECIFICATIONS REFERENCES:
* FASTRAD Functional Specifications, Sections 4.3.1 - 4.3.6
*
*
* DEVELOPMENT HISTORY:
* Date Name Change Release Description
* ID
* -------- ------------- ------ ------- --------------------------
* 09-16-93 J. Programmer 1 Prolog and PDL
* 10-25-93 J. Programmer 1 Coded
* 11-16-93 J. Programmer 1 Controlled
* 12-02-93 J. Programmer 1 Integrated new RSL
* routines
* 12-20-93 J. Programmer 12 1 Created intermediate
* variables for #define
* arguments of calpvs
* in order to pass
* by address
* 02-15-94 J. Programmer 15 2 Corrected time errors
* using RSL routines
* 05-03-94 J. Programmer 3 Enhancements to RSL
* prototypes
* 05-10-94 J. Programmer 3 Added Earth magnetic
* field read capability
* 05-10-94 J. Programmer 3 Added ephemeris read
* capability
*
*
ALGORITHM
*
*
* DO CASE of reference type
*
* CASE 1 or 2, request is for s/c position or velocity vectors
*
* IF offset between request time and time of last calculated s/c
* position and velocity vectors exceeds wait time THEN
*
* COMPUTE elapsed seconds between epoch time and request time
*
* IF ephemeris method is for reading file THEN
* CALL c_ephemrd to read ephemeris file getting s/c position and
* velocity vectors
* ELSE (analytic computation)
* CALL c_calpvs to generate new s/c position and velocity vectors
* ENDIF
*
* SET new time of last calculated s/c position and velocity
* vectors to request time
*
* ENDIF
*
* IF reference type is for s/c position vector THEN
* SET return vector to s/c position vector
* ELSE
* SET return vector to s/c velocity vector
* ENDIF
*
* CASE 3, request is for s/c to Sun unit vector
*
* IF offset between request time and time of last calculated s/c to
* Sun unit vector exceeds wait time THEN
*
* CALL SecsToCalendar c_packst and c_calmjd to get modified Julian
* date
* CALL c_sunlunp to generate new Earth to Sun vector
* CALL GetSun to compute new s/c to Sun unit vector
*
* SET new time of last calculated s/c to Sun unit vector to
* request time
*
* ENDIF
* SET return vector to s/c to Sun unit vector
*
* CASE 4 or 5, request is for Earth magnetic field vector or Earth
* magnetic field unit vector
*
* IF offset between request time and time of last calculated Earth
* magnetic field vector exceeds wait time THEN
*
* CALL SecsToCalendar c_packst and c_calmjd to get modified Julian
* date
* CALL c_jgrenha to get the Greenwich Hour Angle
* CALL c_emagfld to generate new Earth magnetic field vector
* CALL c_unvec3 to SET Earth magnetic field unit vector
*
* SET new time of last calculated Earth magnetic field vector to
* request time
*
* ENDIF
*
* IF reference type is for Earth magnetic field vector THEN
* SET return vector to Earth magnetic field vector
* ELSE
* SET return vector to Earth magnetic field unit vector
* ENDIF
*
* CASE 6, request is for orbit normal unit vector
*
* IF offset between request time and time of last calculated orbit
* normal unit vector exceeds wait time THEN
*
* CALL GetOrbitNormal to generate new orbit normal unit vector
*
* SET new time of last calculated orbit normal unit vector to
* request time
*
* ENDIF
*
* SET return vector to orbit normal unit vector
*
* CASE 7, request is for s/c to Moon unit vector
*
* IF offset between request time and time of last calculated s/c to
* Moon unit vector exceeds wait time THEN
*
* CALL SecsToCalendar c_packst and c_calmjd to get modified Julian
* date
* CALL c_sunlunp to generate new Earth to Moon vector
* CALL GetMoon to compute new s/c to Moon unit vector
*
* SET new time of last calculated s/c to Moon unit vector to request
* time
*
* ENDIF
*
* SET return vector to s/c to Moon unit vector
*
* CASE 8, request is for s/c to Earth unit vector
*
* IF offset between request time and time of last calculated s/c to
* Earth unit vector exceeds wait time THEN
*
* CALL GetEarth to compute new s/c to Earth unit vector
*
* SET new time of last calculated s/c to Earth unit vector to
* request time
*
* ENDIF
*
* SET return vector to s/c to Earth unit vector
*
* END CASE
*
* RETURN
*
******************************************************************** */
/* Include global parameters*/
#include "HD_debug.h"
#include "HD_reference.h"
/* Declare Prototypes */
void c_ephemrd (long , long , long , double , double *,
double *, double *, double *, long *);
void c_calpvs (double , double , double *, double , double ,
long , double *, double *, long *);
void c_sunlunp (double , double , double *, double *);
void c_emagfl2 (long , double , double , double , double *,
long , double *, long *);
void c_nmlist (long , long * , char * , long *);
void c_packst (double , double *);
void c_calmjd (double *, double *);
void c_jgrenha (double , double , long , long , double *,
long *);
void c_unvec3 (double *, double *, double *);
void GetSun (double[3], double[3]);
void GetOrbitNormal(double[3]);
void GetEarth (double[3]);
void GetMoon (double[3], double[3]);
double SecsToCalendar(double);
/*******************************************************************
*
* FUNCTION NAME: GetReference
*
* ARGUMENT LIST:
*
* Argument Type IO Description
* ------------- -------- -- ---------------------------------
* ref_type int I Type of reference data requested
* = 1, S/C position vector
* = 2, S/C velocity vector
* = 3, S/C to Sun unit vector
* = 4, Earth magnetic field
* vector
* = 5, Earth magnetic field unit
* vector
* = 6, Orbit normal unit vector
* = 7, S/C to Moon unit vector
* = 8, S/C to Earth unit vector
* t_request double I Time of requested reference
* vector
* t_wait double I Wait time between reference
* vector calculations
* ref_vector double[3] O Requested reference vector
*
* RETURN VALUE: void
*
*******************************************************************/
void GetReference(int ref_type, double t_request, double t_wait,
double ref_vector[3])
{
/* LOCAL VARIABLES:
*
* Variable Type Description
* ------------- ------- -------------------------------------
* sun double[3] Earth to Sun vector [km] (from
* c_sunlunp)
* moon double[3] Earth to Moon vector [km] (from
* c_sunlunp)
* caldate double Epoch time in calendar format
* starray double[6] Epoch time in unpacked array format
* mjd double Modified Julian Date [days]
* gha double Greenwich Hour Angle [rad]
* a1diff double A.1 - UT1 time difference [sec]
* numselc long Number of secular terms of nutation
* to compute (1- 39, nominally 1)
* numterm long Number of nonsecular terms of
* nutation to compute (1-106,
* nominally 50)
* fdumm double Unused return value (from c_unvec3)
* ierr long Return code from RSL routines
* m double Variable for #defined MU_E
* t double Variable for #defined THREEB
* eptime double Elapsed seconds between epoch time
* and requested time [sec]
* dpos double Array of dummy position vectors used
* by ephemris read routine
* dvel double Array of dummy velocity vectors used
* by ephemris read routine
* loop_counter int Loop counter
* i int Loop counter
* j int Loop counter
*/
double int sun[3], moon[3], caldate, starray[6], mjd, gha, a1diff,
fdumm;
double int m, t;
double int eptime;
long int numselc, numterm;
long int ierr = -100;
long int two = 2;
long int four = 4;
long int zero = 0;
int int i,j;
char *mag_path = "/public/libraries/rsl/hpux/emag1990.dat";
static int loop_counter = 0
static double int dpos[3][100], dvel[3][100];
/* Initialize local parameters for RSL routines */
a1diff = 0.0;
numselc = 1;
numterm = 50;
if (debug_level[RF] > TRACE)
fprintf(debug_file_handle,"ENTER GetReference\n");
if (debug_level[RF] > INPUT)
{
fprintf(debug_file_handle,"\tINPUT\n");
switch (ref_type)
{
case 1:
fprintf(debug_file_handle,
"\t\treference type (ref_type = 1) S/C position vector\n");
break;
case 2:
fprintf(debug_file_handle,
"\t\treference type (ref_type = 2) S/C velocity vector\n");
break;
case 3:
fprintf(debug_file_handle,
"\t\treference type (ref_type = 3) S/C to Sun unit
vector\n");
break;
case 4:
fprintf(debug_file_handle,
"\t\treference type (ref_type = 4) Earth mag field
vector\n");
break;
case 5:
fprintf(debug_file_handle,
"\t\treference type (ref_type = 5) Earth mag field unit
vector\n");
break;
case 6:
fprintf(debug_file_handle,
"\t\treference type (ref_type = 6) Orbit normal unit
vector\n");
break;
case 7:
fprintf(debug_file_handle,
"\t\treference type (ref_type = 7) S/C to Moon unit
vector\n");
break;
case 8:
fprintf(debug_file_handle,
"\t\treference type (ref_type = 8) S/C to Earth unit
vector\n");
break;
}
fprintf(debug_file_handle,
"\t\trequest time [sec] (t_request) = %lf\n",t_request);
fprintf(debug_file_handle,
"\t\twait time [sec] (t_wait) = %lf\n",t_wait);
}
/* Begin Case of reference type */
switch (ref_type)
{
/* Perform case for either s/c position or velocity vector request
* using the RSL routine c_calpvs */
case 1:
case 2:
if (debug_level[RF] > INPUT)
{
fprintf(debug_file_handle,
"\t\tlast pos and vel vector time [sec] (t_rv_ref) =
%lf\n",
t_rv_ref);
fprintf(debug_file_handle,
"\t\tephemeris read method flag (ephem_method) = %c\n",
ephem_method);
}
if ((t_request - t_rv_ref) > t_wait)
{
eptime = t_request - orbital_t_epoch;
if (debug_level[RF] > INTERMEDIATE)
{
fprintf(debug_file_handle, "\tINTERMEDIATE\n");
fprintf(debug_file_handle,
"\t\tRequest time [secs from reference]
(eptime) = %lf\n",eptime);
}
if (ephem_method == 'F')
{
if (loop_counter == 0)
{
for (i=0; i<100; i++)
for (j=0; j<3; j++)
{
dpos[j][i] = 0.0;
dvel[j][i] = 0.0;
}
loop_counter++;
}
c_ephemrd(ephem_file_lu,four,zero,eptime,
dpos,dvel, s_c_pos,s_c_vel,&ierr);
if (ierr)
if (debug_level[RF] > TRACE)
fprintf(debug_file_handle,
"**** Error code from c_ephemrd = %ld\n",ierr);
}
else
{
m = MU_E;
t = THREEB;
c_calpvs(eptime,m,keplerian,t,ttol,maxit,s_c_pos,s_c_vel,&ierr);
if (ierr)
if (debug_level[RF] > TRACE)
fprintf(debug_file_handle,
"**** Error code from c_calpvs = %ld\n",ierr);
if (debug_level[RF] > INTERMEDIATE)
{
fprintf(debug_file_handle,
"\t\tEarth gravitational constant [km^3/sec^2]
(MU_E) = %lf\n",MU_E);
fprintf(debug_file_handle,
"\t\tGrav. constant [Km^2]
(THREEB) = %lf\n",THREEB);
fprintf(debug_file_handle,
"\t\ttolerance of true anomaly [rad]
(ttol) = %lf\n",ttol);
fprintf(debug_file_handle,
"\t\tmax iters of true anomaly
(maxit) = %d\n",maxit);
fprintf(debug_file_handle,
"\t\ttime of request [sec from epoch]
(eptime) = %lf\n",eptime);
fprintf(debug_file_handle,
"\t\tsemi major axis [km]
(keplerian[1]) = %lf\n",keplerian[0]);
fprintf(debug_file_handle,
"\t\teccentricity
(keplerian[2]) = %lf\n",keplerian[1]);
fprintf(debug_file_handle,
"\t\tinclination [rad]
(keplerian[3]) = %lf\n",keplerian[2]);
fprintf(debug_file_handle,
"\t\tra of asc node [rad]
(keplerian[4]) = %lf\n",keplerian[3]);
fprintf(debug_file_handle,
"\t\targ of perigee [rad]
(keplerian[5]) = %lf\n",keplerian[4]);
fprintf(debug_file_handle,
"\t\tmean anomaly [rad]
(keplerian[6]) = %lf\n",keplerian[5]);
}
}
t_rv_ref = t_request;
if (debug_level[RF] > INTERMEDIATE)
{
fprintf(debug_file_handle,
"\t\ts/c position vector [km] (s_c_pos) =
%lf,%lf,%lf\n",
s_c_pos[0],s_c_pos[1],s_c_pos[2]);
fprintf(debug_file_handle,
"\t\ts/c velocity vector [km] (s_c_vel) =
%lf,%lf,%lf\n",
s_c_vel[0],s_c_vel[1],s_c_vel[2]);
}
}
if (ref_type == 1)
for (i=0 ; i<3 ; i++)
ref_vector[i] = s_c_pos[i];
else
for (i=0 ; i<3 ; i++)
ref_vector[i] = s_c_vel[i];
break;
/* Perform case for s/c to Sun unit vector request using the RSL
* routine c_sunlunp */
case 3:
if (debug_level[RF] > INPUT)
fprintf(debug_file_handle,
"\t\tlast sun vector time [sec] (t_s_ref) = %lf\n",t_s_ref);
if ((t_request - t_s_ref) > t_wait)
{
caldate = SecsToCalendar(t_request);
c_packst (caldate,starray);
c_calmjd (starray,&mjd);
c_sunlunp(mjd,t_request,sun,moon);
GetSun (sun,s_pos);
t_s_ref = t_request;
if (debug_level[RF] > INTERMEDIATE)
{
fprintf(debug_file_handle, "\tINTERMEDIATE\n");
fprintf(debug_file_handle,
"\t\tModified Julian Date [days] (mjd) = %lf\n", mjd);
fprintf(debug_file_handle,
"\t\ttime of request [sec] (use t_request see above)
\n");
}
}
for (i=0 ; i<3 ; i++)
ref_vector[i] = s_pos[i];
break;
/* Perform case for Earth magnetic field vector or Earth magnetic
* field unit vector using RSL routines c_emagfld and c_unvec3 */
case 4:
case 5:
if (debug_level[RF] > INPUT)
fprintf(debug_file_handle,
"\t\tlast Earth mag field vector time [sec]
(t_b_ref) = %lf\n", t_b_ref);
if ((t_request - t_b_ref) > t_wait)
{
caldate = SecsToCalendar(t_request);
c_packst (caldate,starray);
c_calmjd (starray,&mjd);
c_jgrenha(mjd,a1diff,numselc,numterm,&gha,&ierr);
if (ierr)
if (debug_level[RF] > TRACE)
fprintf(debug_file_handle,
"**** Error code from c_jgrenha = %ld\n",ierr);
c_nmlist(1,&two,mag_path,&ierr);
if (ierr)
if (debug_level[RF] > TRACE)
fprintf(debug_file_handle,
"**** Error code from c_nmlist = %ld\n",ierr);
c_emagfl2(two,mjd,t_request,gha,s_c_pos,m_order,mag_field,&ierr);
if (ierr)
if (debug_level[RF] > TRACE)
fprintf(debug_file_handle,
"**** Error code from c_emagfl2 = %ld\n",ierr);
c_unvec3 (mag_field,mag_field_unit,&fdumm);
t_b_ref = t_request;
if (debug_level[RF] > INTERMEDIATE)
{
fprintf(debug_file_handle,"\tINTERMEDIATE\n");
fprintf(debug_file_handle,
"\t\tModified Julian Date [days] (mjd) = %lf\n", mjd);
fprintf(debug_file_handle,
"\t\ttime difference [sec] (a1diff) = %lf\n", a1diff);
fprintf(debug_file_handle,
"\t\tnutation number (numselc) = %d\n", numselc);
fprintf(debug_file_handle,
"\t\tnutation number (numterm) = %d\n", numterm);
fprintf(debug_file_handle,
"\t\tGreenwich Hour Angle [rad] (gha) = %lf\n", gha);
fprintf(debug_file_handle,
"\t\torder of magnetic field (m_order) = %d\n",
m_order);
fprintf(debug_file_handle,
"\t\ts/c position vector [km] (s_c_pos) =
%lf,%lf,%lf\n",
s_c_pos[0],s_c_pos[1],s_c_pos[2]);
fprintf(debug_file_handle,
"\t\ttime of request [sec] (use t_request see above)
\n");
}
}
if (ref_type == 4)
for (i=0 ; i<3 ; i++)
ref_vector[i] = mag_field[i];
else
for (i=0 ; i<3 ; i++)
ref_vector[i] = mag_field_unit[i];
break;
/* Perform case for orbit normal unit vector request */
case 6:
/* Debug : Intermediate */
if (debug_level[RF] > INPUT)
fprintf(debug_file_handle,
"\t\tlast normal unit vector time [sec] (t_o_ref) = %lf\n",
t_o_ref);
if ((t_request - t_o_ref) > t_wait)
{
GetOrbitNormal(orbit_normal);
t_o_ref = t_request;
}
for (i=0 ; i<3 ; i++)
ref_vector[i] = orbit_normal[i];
break;
/* Perform case for s/c to Moon unit vector request using the RSL
* routine c_sunlunp */
case 7:
if (debug_level[RF] > INPUT)
fprintf(debug_file_handle,
"\t\tlast moon vector time [sec] (t_m_ref) = %lf\n",t_m_ref);
if ((t_request - t_m_ref) > t_wait)
{
caldate = SecsToCalendar(t_request);
c_packst (caldate,starray);
c_calmjd (starray,&mjd);
c_sunlunp(mjd,t_request,sun,moon);
GetMoon (moon,m_pos);
t_m_ref = t_request;
if (debug_level[RF] > INTERMEDIATE)
{
fprintf(debug_file_handle,"\tINTERMEDIATE\n");
fprintf(debug_file_handle,
"\t\tModified Julian Date [days] (mjd) = %lf\n", mjd);
fprintf(debug_file_handle,
"\t\ttime of request [sec] (use t_request see above) \n");
}
}
for (i=0 ; i<3 ; i++)
ref_vector[i] = m_pos[i];
break;
/* Perform case for s/c to Earth unit vector request */
case 8:
if (debug_level[RF] > INPUT)
fprintf(debug_file_handle,
"\t\tlast Earth vector time [sec] (t_e_ref) = %lf\n",t_e_ref);
if ((t_request - t_e_ref) > t_wait)
{
GetEarth(e_pos);
t_e_ref = t_request;
}
for (i=0 ; i<3 ; i++)
ref_vector[i] = e_pos[i];
break;
} /* end switch */
if (debug_level[RF] > OUTPUT)
{
fprintf(debug_file_handle,"\tOUTPUT\n");
fprintf(debug_file_handle,
"\t\trequested reference vector (ref_vector) = %lf,%lf,%lf\n",
ref_vector[0],ref_vector[1],ref_vector[2]);
}
if (debug_level[RF] > TRACE)
fprintf(debug_file_handle,"EXIT GetReference\n\n");
return;
} /* end */