Solid Angles
// In this program, we first determine the sub-observer points, positions and
// solid angles of all the planets at the current time. Then, to show more
// functions, we define the physical parameters for Pan (the moon of Saturn)
// and then find out some stuff about it at the current time.
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <aephem.h>
int main(int argc, char *argv[]) {
int i;
double jd_ut, jd_ut1, jd_tt, ra, dec, dist, sa, w, n[3], p[3], e[3], v_e[3];
double j[3], light_t, lat, lon, f, a, b;
struct ae_physical_t pan;
// Get the current Julian Date in the required time definitions.
jd_ut = ae_ctime_to_jd(time(NULL));
jd_ut1 = jd_ut + ae_dut1(jd_ut) * AE_D_PER_S;
jd_tt = jd_ut1 + ae_delta_t(jd_ut1) * AE_D_PER_S;
printf("The current Julian Date in UT1 is %0.5f.nn", jd_ut1);
// Get the planets' solid angles.
printf("Here are the planets' data (all geocentric):n");
printf(" -------------------------------------------------------------------"
printf(" Planet R. Ascension Declination Sub-lat.(°) "
printf(" -------------------------------------------------------------------"
for (i = 0; i <= AE_SS_NEPTUNE; i++) {
if (i == AE_SS_EMBARY)
continue; // Don't do the E-M barycentre!
ae_phys_planet[i], &lat, &lon, NULL);
&dec, NULL);
printf(" %-8s " ae_hms_sfmt " " ae_dms_sfmt " %6.2f "
"%7.2fn", ae_ss_name[i], ae_hms_arg(ra), ae_dms_arg(dec), lat,
sa * AE_RTS * AE_RTS);
// Now do some heavy-lifting on Pan. Define its physical parameters.
pan.r_mean = 17.2;
pan.r_eq = 10.4;
pan.r_pole = 9.1;
pan.pole_ra = 40.6;
pan.pole_ra_t = -0.036;
pan.ra_sin_term = ae_phys_no_term;
pan.pole_dec = 83.5;
pan.pole_dec_t = -0.004;
pan.dec_cos_term = ae_phys_no_term;
pan.w = 48.8;
pan.w_d = 626.044;
pan.w_d_sq = 0;
pan.w_sin_term = ae_phys_no_term;
// Get the unit vector pointing from us to Pan. We'll use Saturn as a proxy
// since we don't have ephemerides for Pan itself. We could use
// ae_geocentric_from_orbit() here, but we also need to know the light-travel
// times. First, our coordinates.
ae_kepler(jd_tt, &ae_orb_earth, e);
ae_v_orbit(jd_tt, &ae_orb_earth, v_e);
// Get the light-corrected coordinates of Saturn.
ae_kepler(jd_tt, &ae_orb_saturn, p);
light_t = 0;
for (i = 0; i < 2; i++) {
light_t = ae_light_t(e, p, 0);
ae_kepler(jd_tt - light_t, &ae_orb_saturn, p);
// Now get the geocentric position of Saturn. (We won't be precise and reduce
// to topocentric coordinates. After all, we are pretending that Pan is at
// Saturn's exact position.)
ae_geocentric(jd_tt, p, e, v_e, &ra, &dec, &dist);
// Get the unit vector pointing from us to Pan. We pass a distance of -1 so
// that the unit vector is pointing _towards_ us, not away from us.
ae_polar_to_rect(ra, dec, -1.0, j);
// Get the coordinates of Pan's pole and prime meridean, accounting for
// light-travel time.
ae_phys_pole(jd_ut1 - light_t, jd_tt, &pan, n, &w);
// Get the flattening of the spheroid.
f = ae_flattening(&pan);
// Get the subobserver point.
ae_subobs_point(j, n, w, f, ae_is_retrograde(&pan), &lat, &lon);
// Get the semi-minor axis of the projected disc.
a = pan.r_eq;
b = ae_disc_semiminor(a, pan.r_pole, lat);
// Get its solid angle.
sa = ae_disc_solid_angle(a, b, dist);
// Now print everything out.
printf("Here are some data on Saturn's moon Pan:n");
printf(" Ra/dec: " ae_hms_fmt ", " ae_dms_fmt "n",
ae_hms_arg(ra), ae_dms_arg(dec));
printf(" Distance: %.3f AUn", dist);
printf(" Light-travel time: %.2f minn", light_t * 24.0 * 60.0);
printf(" North-pole coords: [%.2f %.2f %.2f]n", n[0], n[1], n[2]);
printf(" Sub-obs. lon/lat: %.1f, %.1f (deg)n", lon, lat);
printf(" Projected axes: %.1f, %.1f (km)n", a, b);
printf(" Solid angle: %.2g arcsec^2n", sa * AE_RTS * AE_RTS);
return 0;