Implementing calculation of mean elements of the planets based on Sim…

…on et al. - Numerical expressions for precession formulae and mean elements for the Moon and the planets

include/predict/predict.h.in

@@ -378,6 +378,11 @@ double predict_sun_declination(predict_julian_date_t time);
  **/
 double predict_sun_gha(predict_julian_date_t time);
 
+// TODO: documentation
+void predict_observe_venus(const predict_observer_t *observer, predict_julian_date_t time, struct predict_observation *obs);
+
+void predict_observe_jupiter(const predict_observer_t *observer, predict_julian_date_t time, struct predict_observation *obs);
+
 /** 
  * Find next acquisition of signal (AOS) of satellite (when the satellite rises above the horizon). Ignores previous AOS of current pass if the satellite is in range at the start time. 
  *

src/CMakeLists.txt

@@ -1,6 +1,6 @@
-add_library(predict SHARED orbit.c sun.c observer.c sdp4.c sgp4.c refraction.c unsorted.c julian_date.c version.c moon.c)
+add_library(predict SHARED orbit.c sun.c observer.c sdp4.c sgp4.c refraction.c unsorted.c julian_date.c version.c moon.c time.c coordinates.c body.c)
 
-add_library(predict_static STATIC orbit.c sun.c observer.c sdp4.c sgp4.c refraction.c unsorted.c julian_date.c moon.c)
+add_library(predict_static STATIC orbit.c sun.c observer.c sdp4.c sgp4.c refraction.c unsorted.c julian_date.c moon.c time.c coordinates.c body.c)
 
 set_target_properties(predict_static PROPERTIES OUTPUT_NAME predict)
 

src/body.c

@@ -0,0 +1,119 @@
+#include <math.h>
+#include <stdio.h> // debug
+#include "body.h"
+#include "coordinates.h"
+#include "time.h"
+
+// Mean elements of the planets referred to the mean dynamical ecliptic and equinox of date
+// Coefficients from Listings 5.9.*
+// Simon et al. - Numerical expressions for precession formulae and mean elements for the Moon and the planets
+// a, lambda, e, pi, i, Omega
+// Memory Allocation for 7 coefficients per variable
+// If definition is shorter, intrinsic zeros for the rest of the coefficients
+const predict_body_keplarian_elements_t KEPLERIAN_ELEMENTS_VENUS = {
+  {0.7233298200},
+  {181.97980085, 2106691666.31989, 111.65021, 0.05368, -0.23516, -0.00179, 0.0020},
+  {0.0067719164, -0.0004776521, 98127e-10, 4639e-10, 123e-10, -3e-10},
+  {131.56370300, 50477.47081, -387.42545, -20.44048, -0.95948, 0.00044, 0.00020},
+  {3.39466189, 36.13261, -0.31523, -0.02525, 0.00085, -0.00008},
+  {76.67992019, 32437.57636, 146.22586, -0.33446, 0.23007, -0.00088, 0.00009}
+};
+
+const predict_body_keplarian_elements_t KEPLERIAN_ELEMENTS_EARTH = {
+  {1.0000010178},
+  {100.46645683, 1296027711.03429, 109.15809, 0.07207, -0.23530, -0.00180, 0.00020},
+  {0.0167086342, -0.0004203654, -0.0000126734, 1444e-10, -2e-10, 3e-10},
+  {102.93734808, 61900.55290, 164.47797, -0.06365, -0.12090, 0.00298, 0.00020},
+  {0.},
+  {0.}
+};
+
+const predict_body_keplarian_elements_t KEPLERIAN_ELEMENTS_JUPITER = {
+  {5.2026032092, 19132e-10, -39e-10, -60e-10, -10e-10, 1e-10},
+  {34.35151874, 109306899.89453, 80.38700, 0.13327, -0.18850, 0.00411, -0.00014},
+  {0.0484979255, 0.0016322542, -0.0000471366, -20063e-10, 1018e-10, -21e-10, 1e-10},
+  {14.33120687, 58054.86625, 370.95016, -16.07110, 0.51186, -0.02268, 0.00004},
+  {1.30326698, -197.87442, 1.67744, -0.00838, -0.00735, 0.00085, 0.00004},
+  {100.46440702, 36755.18747, 145.13295, 1.45556, -0.59609, -0.04324, 0.00175}
+};
+
+// TODO: add additional bodies if necessary
+
+
+const predict_body_keplarian_elements_t *KEPLERIAN_ELEMENTS[10] = {
+  NULL,
+  NULL, //&KEPLERIAN_ELEMENTS_MERCURY,
+  &KEPLERIAN_ELEMENTS_VENUS,
+  &KEPLERIAN_ELEMENTS_EARTH,
+  NULL, //&KEPLERIAN_ELEMENTS_MARS,
+  &KEPLERIAN_ELEMENTS_JUPITER,
+  NULL, //&KEPLERIAN_ELEMENTS_SATURN,
+  NULL, //&KEPLERIAN_ELEMENTS_URANUS,
+  NULL, //&KEPLERIAN_ELEMENTS_NEPTUNE,
+  NULL, //&KEPLERIAN_ELEMENTS_PLUTO
+};
+
+void predict_body_calc_pos(predict_julian_ephimeris_day_t JDE, predict_body_t body, predict_pos_t *pos) {
+  double t = predict_julian_ephimeris_day_to_centuries(JDE)/10.; // millenia past J2000.0
+  // calculate Keplerian elements - constants
+  double a =      KEPLERIAN_ELEMENTS[body]->a[0];
+  double lambda = KEPLERIAN_ELEMENTS[body]->lambda[0] * 3600;
+  double e =      KEPLERIAN_ELEMENTS[body]->e[0];
+  double pi =     KEPLERIAN_ELEMENTS[body]->pi[0] * 3600;
+  double i =      KEPLERIAN_ELEMENTS[body]->i[0] * 3600;
+  double Omega =  KEPLERIAN_ELEMENTS[body]->Omega[0] * 3600;
+
+  // calculate Keplerian elements - time series
+  int n;
+  double t_pow;
+  for (n=1; n<N_COEFFS; n++) {
+    if (n > 1) {
+      t_pow = pow(t,n);
+    } else {
+      t_pow = t;
+    }
+    a       += KEPLERIAN_ELEMENTS[body]->a[n]       * t_pow;
+    lambda  += KEPLERIAN_ELEMENTS[body]->lambda[n]  * t_pow; // time coefficients in arcsec
+    e       += KEPLERIAN_ELEMENTS[body]->e[n]       * t_pow;
+    pi      += KEPLERIAN_ELEMENTS[body]->pi[n]      * t_pow; // time coefficients in arcsec
+    i       += KEPLERIAN_ELEMENTS[body]->i[n]       * t_pow; // time coefficients in arcsec
+    Omega   += KEPLERIAN_ELEMENTS[body]->Omega[n]   * t_pow; // time coefficients in arcsec
+  }
+
+  // Keplerian elements from arcsec to radians
+  lambda  = lambda / 3600 / 180*M_PI;
+  pi      = pi     / 3600 / 180*M_PI;
+  i       = i      / 3600 / 180*M_PI;
+  Omega   = Omega  / 3600 / 180*M_PI;
+
+  double omega = pi - Omega; // argument of perihelion
+  double M = lambda - pi; // mean anomaly
+  // for 3000 BC to 3000 AD + bT^2 + c * cos(f*T) + s * sin(f*T)
+
+  // Solving Kepler's Equation
+  // Meeus (30.7)
+  double E = M;
+  double delta_E = 1;
+  double delta_M = 0;
+
+  while (fabs(delta_E) > 1e-10) { // 1e-6 tolerance sufficient for mean ephimerides stuff
+    delta_M = M + e * sin(E) - E;
+    delta_E = delta_M / (1 - e * cos(E));
+    E += delta_E;
+  }
+
+  double x_prime = a * (cos(E) - e);
+  double y_prime = a * sqrt(1-pow(e,2)) * sin(E);
+
+  double x_ecl = (cos(omega)*cos(Omega) - sin(omega)*sin(Omega)*cos(i)) * x_prime + (-sin(omega)*cos(Omega) - cos(omega)*sin(Omega)*cos(i))* y_prime;
+  double y_ecl = (cos(omega)*sin(Omega) + sin(omega)*cos(Omega)* cos(i)) * x_prime + (-sin(omega)*sin(Omega) + cos(omega)*cos(Omega)*cos(i))* y_prime;
+  double z_ecl = (sin(omega) * sin(i)) * x_prime + (cos(omega) * sin(i)) * y_prime;
+
+  pos->coords[0] = x_ecl;
+  pos->coords[1] = y_ecl;
+  pos->coords[2] = z_ecl;
+  pos->type = CS_CARTESIAN;
+  pos->origin = CS_SUN;
+  pos->plane = CS_ECLIPTIC;
+
+}

src/body.h

@@ -0,0 +1,33 @@
+#ifndef _BODY_H_
+#define _BODY_H_
+
+#include "coordinates.h"
+#include "time.h"
+
+typedef enum {
+  SUN       = 0,
+  MERCURY   = 1,
+  VENUS     = 2,
+  EARTH     = 3,
+  MARS      = 4,
+  JUPITER   = 5,
+  SATURN    = 6,
+  URANUS    = 7,
+  NEPTUNE   = 8,
+  PLUTO     = 9
+} predict_body_t;
+
+#define N_COEFFS 7
+
+typedef struct {
+  const double a[N_COEFFS];       // semi-major axis [au]
+  const double lambda[N_COEFFS];  // mean longitude []
+  const double e[N_COEFFS];       // eccentricity
+  const double pi[N_COEFFS];      // longitude of the perihelion
+  const double i[N_COEFFS];       // inclination
+  const double Omega[N_COEFFS];   // longitude of the ascending node
+} predict_body_keplarian_elements_t;
+
+void predict_body_calc_pos(predict_julian_ephimeris_day_t JDE, predict_body_t body, predict_pos_t *pos);
+
+#endif

src/libpredict.symver

@@ -27,6 +27,8 @@ VER_2.0 {
 		predict_sun_ra;
 		predict_sun_declination;
 		predict_sun_gha;
+		predict_observe_venus;
+		predict_observe_jupiter;
 		predict_next_aos;
 		predict_next_los;
 		predict_at_max_elevation;

src/observer.c

@@ -2,8 +2,11 @@
 #include "unsorted.h"
 #include <stdlib.h>
 #include <string.h>
+#include "body.h"
+#include "coordinates.h"
 #include "defs.h"
 #include "sun.h"
+#include "time.h"
 
 void observer_calculate(const predict_observer_t *observer, double time, const double pos[3], const double vel[3], struct predict_observation *result);
 
@@ -130,19 +133,99 @@ void observer_calculate(const predict_observer_t *observer, double time, const d
 	// Elevation rate
 	double x_dot = (top_z_dot*range_length - range_rate_length*top_z) / (range_length * range_length);
 	double el_dot = x_dot / sqrt( 1 - x*x );
-	
+
 	result->azimuth = az;
 	result->azimuth_rate = az_dot;
 	result->elevation = el;
 	result->elevation_rate = el_dot;
 	result->range = range_length;
-	result->range_rate = range_rate_length; 
+	result->range_rate = range_rate_length;
 	result->range_x = range[0];
 	result->range_y = range[1];
 	result->range_z = range[2];
 
 }
 
+
+/**
+ * Calculate az,el coordinates of Venus with respect to an observer
+ *
+ * \param observer Observer in Geographic Coordinate System (WGS84) (Latitude, Longitude, Altitude)
+ * \param time Observation time (libpredict julian date)
+ * \param obs Observation in Horizontal Coordinate System (HCS) (Azimuth, Elevation)
+ */
+void predict_observe_venus(const predict_observer_t *observer, predict_time_t time, struct predict_observation *obs)
+{
+	predict_julian_day_t JD = predict_time_to_julian_day(time);
+	predict_julian_ephimeris_day_t JDE = predict_julian_day_to_ephimeris_day(JD);
+
+	/* calculate ecliptic coordinates of Venus and Earth (heliocentric) */
+	predict_pos_t venus, earth;
+	predict_body_calc_pos(JDE, VENUS, &venus);
+	predict_body_calc_pos(JDE, EARTH, &earth); // XXX: Earth Moon Barycenter
+
+	/* calculate ecliptic coordinates of Venus (geocentric) */
+	predict_helio_to_geo(&venus, &earth);
+
+	/* ecliptic - rectangular to spheric */
+	predict_cart_to_sph(&venus);
+
+	// correcting for effect of light time
+	double time_offset = 0.0057755183 * venus.coords[0]; // in days
+	predict_body_calc_pos(JDE-time_offset, VENUS, &venus);
+	predict_helio_to_geo(&venus, &earth);
+	predict_cart_to_sph(&venus);
+
+	/* ecliptic to equatorial - J2000.0 */
+	predict_ecliptic_to_equatorial(JDE, &venus);
+
+	/* equatorial to horizontal */
+	predict_equatorial_to_horizontal(JD, &venus, observer);
+
+	obs->azimuth = venus.coords[2];
+	obs->elevation = venus.coords[1];
+}
+
+/**
+ * Calculate az,el coordinates of Jupiter with respect to an observer
+ *
+ * \param observer Observer in GCS (Latitude, Longitude, Altitude)
+ * \param time Observation time TODO: which reference
+ * \param obs Observation in HCS (Azimuth, Elevation)
+ */
+void predict_observe_jupiter(const predict_observer_t *observer, predict_time_t time, struct predict_observation *obs)
+{
+	predict_julian_day_t JD = predict_time_to_julian_day(time);
+	predict_julian_ephimeris_day_t JDE = predict_julian_day_to_ephimeris_day(JD);
+
+	/* calculate ecliptic coordinates of Jupiter and Earth (heliocentric) */
+	predict_pos_t jupiter, earth;
+	predict_body_calc_pos(JDE, JUPITER, &jupiter);
+	predict_body_calc_pos(JDE, EARTH, &earth);
+
+	/* calculate ecliptic coordinates of Jupiter (geocentric) */
+	predict_helio_to_geo(&jupiter, &earth);
+
+	/* ecliptic - rectangular to spheric */
+	predict_cart_to_sph(&jupiter);
+
+	// correcting for effect of light time
+	double time_offset = 0.0057755183 * jupiter.coords[0]; // in days
+
+	predict_body_calc_pos(JDE-time_offset, JUPITER, &jupiter);
+	predict_helio_to_geo(&jupiter, &earth);
+	predict_cart_to_sph(&jupiter);
+
+	/* ecliptic to equatorial - J2000.0 */
+	predict_ecliptic_to_equatorial(JDE, &jupiter);
+
+	/* equatorial to horizontal */
+	predict_equatorial_to_horizontal(JD, &jupiter, observer);
+
+	obs->azimuth = jupiter.coords[2];
+	obs->elevation = jupiter.coords[1];
+}
+
 struct predict_observation predict_next_aos(const predict_observer_t *observer, const predict_orbital_elements_t *orbital_elements, double start_utc)
 {
 	double curr_time = start_utc;