/* Program geom -- Compute solar geometric quantities This is a demo program written for Ashitey Trebi-Ollenu. The program demonstrates some simple geometric computations using CSPICE. The following items are computed: 1) Sun position as seen from a specified surface point on the earth, optionally corrected for light time and stellar aberration. The sun state is defined relative to the fixed surface point DSS-16. 2) Sun azimuth and elevation derived from (1). 3) Geometric sub-solar point on the earth. 4) Geometric sub-solar point on Mars. 5) Geometric sun position as seen from a surface point on Mars, where the surface point is defined in the metakernel. 6) Sun azimuth and elevation derived from (5). All sub-solar points are computed using the intercept of the sun-target body line and the target's surface. This program requires the following kernels: - Leapseconds. - Planetary SPK file. This provides the ephemeris of the earth and sun. - SPK file giving the position of DSS-16 in ITRF93 coordinates. - IAU PCK file. This provides radii of the earth, which are required for computation of geodetic coordinates. - High-precision earth PCK relating the terrestrial frame TRF93 to an inertial frame. This is needed because the DSS-16 ephemeris is given relative to ITRF93. - Frame kernel defining the topocentric frame centered at DSS-16. All of the kernels are loaded via the metakernel geom.metaker. The code has not been subjected to rigorous testing. Version 1.0.0 07-SEP-2000 (NJB) */ #include #include "SpiceUsr.h" int main() { /* Prototypes */ SpiceBoolean badkpv_c ( ConstSpiceChar *caller, ConstSpiceChar *name, ConstSpiceChar *comp, SpiceInt size, SpiceInt divby, SpiceChar type ); /* Local constants */ #define CORRLEN 10 #define TIMELEN 33 #define LINELEN 81 #define METAKER "geom.metaker" #define OBSERVER "DSS-16" #define TARGET "sun" #define REF "DSS-16_TOPO" #define MARS_PT_COORDS "MARS_SURFACE_PT_COORDS" /* Local variables */ SpiceBoolean found; SpiceChar abcorr [ CORRLEN ]; SpiceChar utc [ TIMELEN ]; SpiceDouble az; SpiceDouble el; SpiceDouble earthRadii [3]; SpiceDouble et; SpiceDouble f; SpiceDouble gpAlt; SpiceDouble gpLat; SpiceDouble gpLon; SpiceDouble lat; SpiceDouble lon; SpiceDouble lt; SpiceDouble marsF; SpiceDouble marsPt [3]; SpiceDouble marsPtCoords [3]; SpiceDouble marsRe; SpiceDouble marsRp; SpiceDouble marsSunAz; SpiceDouble marsSunEl; SpiceDouble marsSunGp [3]; SpiceDouble marsSunGpAlt; SpiceDouble marsSunGpLat; SpiceDouble marsSunGpLon; SpiceDouble marsSunPos [3]; SpiceDouble marsRadii [3]; SpiceDouble marsZ [3] = { 0.0, 0.0, 1.0 }; SpiceDouble range; SpiceDouble re; SpiceDouble rp; SpiceDouble sunGp [3]; SpiceDouble sunPos [3]; SpiceDouble topoMat [3][3]; SpiceDouble topoX [3]; SpiceDouble topoY [3]; SpiceDouble topoZ [3]; SpiceInt n; /* Check in: let the SPICE tracing system know which module is at the top of the call tree. */ chkin_c ( "geom" ); /* Load the SPICE kernels we'll need. */ furnsh_c ( METAKER ); /* Prompt the user for the UTC time of interest. Convert the input UTC time to ephemeris time. */ prompt_c ( "Enter UTC epoch > ", TIMELEN, utc ); str2et_c ( utc, &et ); /* Ask the user which aberration correction to use. See the source file spkpos_c.c for an explanation of aberration corrections. NOTE: to get a state vector, i.e., position and velocity, use spkezr_c. Then look up the position of the sun relative to our surface location, DSS-16. Get the position in the topocentric frame centered at the station. In addition to the position vector, the one-way light time is returned as a bonus. */ prompt_c ( "Enter aberration correction > ", CORRLEN, abcorr ); spkpos_c ( TARGET, et, REF, abcorr, OBSERVER, sunPos, < ); printf ( "The position of %s relative to %s\n" "in reference frame %s\n" "using aberration correction %s\n" "at UTC epoch %s is:\n\n" "x (km): %25.17e\n" "y (km): %25.17e\n" "z (km): %25.17e\n", TARGET, OBSERVER, REF, abcorr, utc, sunPos[0], sunPos[1], sunPos[2] ); /* Now find the azimuth and elevation of the sun. Converting the sun's position to latitudinal coordinates gets us most of the way there. The resulting angles are given in radians. */ reclat_c ( sunPos, &range, &lon, &lat ); el = lat; az = -lon; /* The x-axis of the topocentric frame points north, and the "longitude" we just found is measured in the counterclockwise direction from this axis. So, negating longitude yields the azimuth angle. Convert the angles to degrees on output. */ printf ( "\n" "Azimuth of the sun (deg): %13.6f\n" "Elevation of the Sun (deg): %13.6f\n", az * dpr_c(), el * dpr_c() ); /* Now find the geographic position of the sun. Equivalently, we want to find the sub-solar point on the earth, without using aberration corrections. The observer is not relevant in this situation, but one is required by the function subsol_c, so we make the sun the observer. */ subsol_c ( "Intercept", "earth", et, "NONE", "sun", sunGp ); /* Find the geodetic coordinates of the sub-solar point. First, look up the earth's radii from the kernel pool. Find the flattening factor. */ bodvar_c ( 399, "RADII", &n, earthRadii ); re = earthRadii[0]; rp = earthRadii[1]; f = ( re - rp ) / re; recgeo_c ( sunGp, re, f, &gpLon, &gpLat, &gpAlt ); /* Convert the angles to degrees on output. */ printf ( "\n" "Geodetic longitude of the sub-solar point on earth (deg): %13.6f\n" "Geodetic longitude of the sub-solar point on earth (deg): %13.6f\n", gpLon * dpr_c(), gpLat * dpr_c() ); /* Now we'll do some computations using Mars as the "observing" body. First, we'll find the geographic position of the sun. */ subsol_c ( "Intercept", "mars", et, "NONE", "sun", marsSunGp ); /* Look up the radii of Mars. The bodvar_c shortcut can be used here. We'll need the flattening factor for Mars, so compute it here. */ bodvar_c ( 499, "RADII", &n, marsRadii ); marsRe = marsRadii[0]; marsRp = marsRadii[1]; marsF = ( marsRe - marsRp ) / marsRe; /* Find the geodetic coordinates of the sub-solar point. */ recgeo_c ( marsSunGp, marsRe, marsF, &marsSunGpLon, &marsSunGpLat, &marsSunGpAlt ); /* Convert the angles to degrees on output. */ printf ( "\n" "Geodetic longitude of the sub-solar point on mars (deg): %25.17f\n" "Geodetic latitude of the sub-solar point on mars (deg): %25.17f\n", marsSunGpLon * dpr_c(), marsSunGpLat * dpr_c() ); /* Finally, we'll compute the position of the sun, as well as the sun's azimuth and elevation, as seen from an arbitrary point on Mars. In this case, we can't treat the observer as an ephemeris object, so we'll have to manually carry out some of the work that spkezr_c was doing for us when the observer was DSS-16. We'll compute two versions of the sun's position vector: a geometric vector and a vector corrected for light time and stellar aberration. First step: look up the planetodetic coordinates of the surface point. These coordinates have been loaded via the metakernel. The routine badkpv_c checks to ensure that specified kernel variables are present in the kernel pool and have the expected dimension and data type. The combination of badkpv_c and gdpool_c provides a more general, but slightly more complex, interface than does bodvar_c. */ badkpv_c ( "geom", MARS_PT_COORDS, "=", 3, 1, 'N' ); /* If we got past the badkpv_c check, the desired variable is present in the kernel pool, so it's not necessary to check the found flag. */ gdpool_c ( MARS_PT_COORDS, 0, 3, &n, marsPtCoords, &found ); /* The coordinates we looked up are expected to be planetodetic longitude, latitude (both in degrees), and altitude (in km). Find the equivalent rectangular coordinates. Note that georec_c requires angular inputs in radians. */ georec_c ( rpd_c()*marsPtCoords[0], rpd_c()*marsPtCoords[1], marsPtCoords[2], marsRe, marsF, marsPt ); printf ( "\n" "Geodetic longitude of mars surface point (deg): %25.17f\n" "Geodetic latitude of mars surface point (deg): %25.17f\n", marsPtCoords[0], marsPtCoords[1] ); /* Look up the geometric position of the sun as seen from Mars' center. Represent this vector in Mars body-fixed coordinates. Subtract of the position of the surface point to obtain the position of the sun relative to the surface point. */ spkpos_c ( "sun", et, "IAU_MARS", "NONE", "mars", marsSunPos, < ); vsub_c ( marsSunPos, marsPt, marsSunPos ); /* Construct the basis vectors of a topocentric frame centered as marsPt. The function surfnm_c gives us the unit surface normal vector at this point. The +y direction is the unitized cross product of normal and the Mars' pole, so +y points west. The +x direction is the cross product of the unit +y vector and the +z vector, so +x points north. */ surfnm_c ( marsRadii[0], marsRadii[1], marsRadii[2], marsPt, topoZ ); ucrss_c ( topoZ, marsZ, topoY ); ucrss_c ( topoY, topoZ, topoX ); /* Assemble the transformation matrix that maps Mars body-fixed coordinates to topocentric coordinates. */ vequ_c ( topoX, topoMat[0] ); vequ_c ( topoY, topoMat[1] ); vequ_c ( topoZ, topoMat[2] ); /* Convert the geometric sun position vector into topocentric coordinates and find the azimuth and elevation of the sun. */ mxv_c ( topoMat, marsSunPos, marsSunPos ); printf ( "\n\nThe geometric position of %s relative to \n" "%s\n" "in topocentric coordinates \n" "using aberration correction %s\n" "at UTC epoch %s is:\n\n" "x (km): %25.17e\n" "y (km): %25.17e\n" "z (km): %25.17e\n", "the sun", "the specified surface point on mars", "NONE", utc, marsSunPos[0], marsSunPos[1], marsSunPos[2] ); reclat_c ( marsSunPos, &range, &lon, &lat ); el = lat; az = -lon; /* The x-axis of the topocentric frame points north, and the "longitude" we just found is measured in the counterclockwise direction from this axis. So, negating longitude yields the azimuth angle. Convert the angles to degrees on output. */ printf ( "\n" "Azimuth of the sun at mars surface point (deg): %25.17f\n" "Elevation of the sun at mars surface point (deg): %25.17f\n", az * dpr_c(), el * dpr_c() ); return ( 0 ); }