gergely
/
swe-glib
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity
swe-glib/swe/src/swemplan.c

969 lines
27 KiB
C
Raw Blame History

/* SWISSEPH
$Header: /home/dieter/sweph/RCS/swemplan.c,v 1.74 2008/06/16 10:07:20 dieter Exp $
Moshier planet routines
modified for SWISSEPH by Dieter Koch
**************************************************************/
/* Copyright (C) 1997 - 2008 Astrodienst AG, Switzerland. All rights reserved.
License conditions
------------------
This file is part of Swiss Ephemeris.
Swiss Ephemeris is distributed with NO WARRANTY OF ANY KIND. No author
or distributor accepts any responsibility for the consequences of using it,
or for whether it serves any particular purpose or works at all, unless he
or she says so in writing.
Swiss Ephemeris is made available by its authors under a dual licensing
system. The software developer, who uses any part of Swiss Ephemeris
in his or her software, must choose between one of the two license models,
which are
a) GNU public license version 2 or later
b) Swiss Ephemeris Professional License
The choice must be made before the software developer distributes software
containing parts of Swiss Ephemeris to others, and before any public
service using the developed software is activated.
If the developer choses the GNU GPL software license, he or she must fulfill
the conditions of that license, which includes the obligation to place his
or her whole software project under the GNU GPL or a compatible license.
See http://www.gnu.org/licenses/old-licenses/gpl-2.0.html
If the developer choses the Swiss Ephemeris Professional license,
he must follow the instructions as found in http://www.astro.com/swisseph/
and purchase the Swiss Ephemeris Professional Edition from Astrodienst
and sign the corresponding license contract.
The License grants you the right to use, copy, modify and redistribute
Swiss Ephemeris, but only under certain conditions described in the License.
Among other things, the License requires that the copyright notices and
this notice be preserved on all copies.
Authors of the Swiss Ephemeris: Dieter Koch and Alois Treindl
The authors of Swiss Ephemeris have no control or influence over any of
the derived works, i.e. over software or services created by other
programmers which use Swiss Ephemeris functions.
The names of the authors or of the copyright holder (Astrodienst) must not
be used for promoting any software, product or service which uses or contains
the Swiss Ephemeris. This copyright notice is the ONLY place where the
names of the authors can legally appear, except in cases where they have
given special permission in writing.
The trademarks 'Swiss Ephemeris' and 'Swiss Ephemeris inside' may be used
for promoting such software, products or services.
*/
#include <string.h>
#include "swephexp.h"
#include "sweph.h"
#include "swephlib.h"
#include "swemptab.h"
#define TIMESCALE 3652500.0
#define mods3600(x) ((x) - 1.296e6 * floor ((x)/1.296e6))
#define FICT_GEO 1
#define KGAUSS_GEO 0.0000298122353216 /* Earth only */
/* #define KGAUSS_GEO 0.00002999502129737 Earth + Moon */
static void embofs_mosh(double J, double *xemb);
static int check_t_terms(double t, char *sinp, double *doutp);
static int read_elements_file(int32 ipl, double tjd,
double *tjd0, double *tequ,
double *mano, double *sema, double *ecce,
double *parg, double *node, double *incl,
char *pname, int32 *fict_ifl, char *serr);
static int pnoint2msh[] = {2, 2, 0, 1, 3, 4, 5, 6, 7, 8, };
/* From Simon et al (1994) */
static double freqs[] =
{
/* Arc sec per 10000 Julian years. */
53810162868.8982,
21066413643.3548,
12959774228.3429,
6890507749.3988,
1092566037.7991,
439960985.5372,
154248119.3933,
78655032.0744,
52272245.1795
};
static double phases[] =
{
/* Arc sec. */
252.25090552 * 3600.,
181.97980085 * 3600.,
100.46645683 * 3600.,
355.43299958 * 3600.,
34.35151874 * 3600.,
50.07744430 * 3600.,
314.05500511 * 3600.,
304.34866548 * 3600.,
860492.1546,
};
static struct plantbl *planets[] =
{
&mer404,
&ven404,
&ear404,
&mar404,
&jup404,
&sat404,
&ura404,
&nep404,
&plu404
};
static double FAR ss[9][24];
static double FAR cc[9][24];
static void sscc (int k, double arg, int n);
int swi_moshplan2 (double J, int iplm, double *pobj)
{
int i, j, k, m, k1, ip, np, nt;
signed char FAR *p;
double FAR *pl, *pb, *pr;
double su, cu, sv, cv, T;
double t, sl, sb, sr;
struct plantbl *plan = planets[iplm];
T = (J - J2000) / TIMESCALE;
/* Calculate sin( i*MM ), etc. for needed multiple angles. */
for (i = 0; i < 9; i++)
{
if ((j = plan->max_harmonic[i]) > 0)
{
sr = (mods3600 (freqs[i] * T) + phases[i]) * STR;
sscc (i, sr, j);
}
}
/* Point to start of table of arguments. */
p = plan->arg_tbl;
/* Point to tabulated cosine and sine amplitudes. */
pl = plan->lon_tbl;
pb = plan->lat_tbl;
pr = plan->rad_tbl;
sl = 0.0;
sb = 0.0;
sr = 0.0;
for (;;)
{
/* argument of sine and cosine */
/* Number of periodic arguments. */
np = *p++;
if (np < 0)
break;
if (np == 0)
{ /* It is a polynomial term. */
nt = *p++;
/* Longitude polynomial. */
cu = *pl++;
for (ip = 0; ip < nt; ip++)
{
cu = cu * T + *pl++;
}
sl += mods3600 (cu);
/* Latitude polynomial. */
cu = *pb++;
for (ip = 0; ip < nt; ip++)
{
cu = cu * T + *pb++;
}
sb += cu;
/* Radius polynomial. */
cu = *pr++;
for (ip = 0; ip < nt; ip++)
{
cu = cu * T + *pr++;
}
sr += cu;
continue;
}
k1 = 0;
cv = 0.0;
sv = 0.0;
for (ip = 0; ip < np; ip++)
{
/* What harmonic. */
j = *p++;
/* Which planet. */
m = *p++ - 1;
if (j)
{
k = j;
if (j < 0)
k = -k;
k -= 1;
su = ss[m][k]; /* sin(k*angle) */
if (j < 0)
su = -su;
cu = cc[m][k];
if (k1 == 0)
{ /* set first angle */
sv = su;
cv = cu;
k1 = 1;
}
else
{ /* combine angles */
t = su * cv + cu * sv;
cv = cu * cv - su * sv;
sv = t;
}
}
}
/* Highest power of T. */
nt = *p++;
/* Longitude. */
cu = *pl++;
su = *pl++;
for (ip = 0; ip < nt; ip++)
{
cu = cu * T + *pl++;
su = su * T + *pl++;
}
sl += cu * cv + su * sv;
/* Latitiude. */
cu = *pb++;
su = *pb++;
for (ip = 0; ip < nt; ip++)
{
cu = cu * T + *pb++;
su = su * T + *pb++;
}
sb += cu * cv + su * sv;
/* Radius. */
cu = *pr++;
su = *pr++;
for (ip = 0; ip < nt; ip++)
{
cu = cu * T + *pr++;
su = su * T + *pr++;
}
sr += cu * cv + su * sv;
}
pobj[0] = STR * sl;
pobj[1] = STR * sb;
pobj[2] = STR * plan->distance * sr + plan->distance;
return OK;
}
/* Moshier ephemeris.
* computes heliocentric cartesian equatorial coordinates of
* equinox 2000
* for earth and a planet
* tjd julian day
* ipli internal SWEPH planet number
* xp array of 6 doubles for planet's position and speed
* xe earth's
* serr error string
*/
int swi_moshplan(double tjd, int ipli, AS_BOOL do_save, double *xpret, double *xeret, char *serr)
{
int i;
int do_earth = FALSE;
double dx[3], x2[3], xxe[6], xxp[6];
double *xp, *xe;
double dt;
char s[AS_MAXCH];
int iplm = pnoint2msh[ipli];
struct plan_data *pdp = &swed.pldat[ipli];
struct plan_data *pedp = &swed.pldat[SEI_EARTH];
double seps2000 = swed.oec2000.seps;
double ceps2000 = swed.oec2000.ceps;
if (do_save) {
xp = pdp->x;
xe = pedp->x;
} else {
xp = xxp;
xe = xxe;
}
if (do_save || ipli == SEI_EARTH || xeret != NULL)
do_earth = TRUE;
/* tjd beyond ephemeris limits, give some margin for spped at edge */
if (tjd < MOSHPLEPH_START - 0.3 || tjd > MOSHPLEPH_END + 0.3) {
if (serr != NULL) {
sprintf(s, "jd %f outside Moshier planet range %.2f .. %.2f ",
tjd, MOSHPLEPH_START, MOSHPLEPH_END);
if (strlen(serr) + strlen(s) < AS_MAXCH)
strcat(serr, s);
}
return(ERR);
}
/* earth, for geocentric position */
if (do_earth) {
if (tjd == pedp->teval
&& pedp->iephe == SEFLG_MOSEPH) {
xe = pedp->x;
} else {
/* emb */
swi_moshplan2(tjd, pnoint2msh[SEI_EMB], xe); /* emb hel. ecl. 2000 polar */
swi_polcart(xe, xe); /* to cartesian */
swi_coortrf2(xe, xe, -seps2000, ceps2000);/* and equator 2000 */
embofs_mosh(tjd, xe); /* emb -> earth */
if (do_save) {
pedp->teval = tjd;
pedp->xflgs = -1;
pedp->iephe = SEFLG_MOSEPH;
}
/* one more position for speed. */
swi_moshplan2(tjd - PLAN_SPEED_INTV, pnoint2msh[SEI_EMB], x2);
swi_polcart(x2, x2);
swi_coortrf2(x2, x2, -seps2000, ceps2000);
embofs_mosh(tjd - PLAN_SPEED_INTV, x2);/**/
for (i = 0; i <= 2; i++)
dx[i] = (xe[i] - x2[i]) / PLAN_SPEED_INTV;
/* store speed */
for (i = 0; i <= 2; i++) {
xe[i+3] = dx[i];
}
}
if (xeret != NULL)
for (i = 0; i <= 5; i++)
xeret[i] = xe[i];
}
/* earth is the planet wanted */
if (ipli == SEI_EARTH) {
xp = xe;
} else {
/* other planet */
/* if planet has already been computed, return */
if (tjd == pdp->teval && pdp->iephe == SEFLG_MOSEPH) {
xp = pdp->x;
} else {
swi_moshplan2(tjd, iplm, xp);
swi_polcart(xp, xp);
swi_coortrf2(xp, xp, -seps2000, ceps2000);
if (do_save) {
pdp->teval = tjd;/**/
pdp->xflgs = -1;
pdp->iephe = SEFLG_MOSEPH;
}
/* one more position for speed.
* the following dt gives good speed for light-time correction
*/
#if 0
for (i = 0; i <= 2; i++)
dx[i] = xp[i] - pedp->x[i];
dt = LIGHTTIME_AUNIT * sqrt(square_sum(dx));
#endif
dt = PLAN_SPEED_INTV;
swi_moshplan2(tjd - dt, iplm, x2);
swi_polcart(x2, x2);
swi_coortrf2(x2, x2, -seps2000, ceps2000);
for (i = 0; i <= 2; i++)
dx[i] = (xp[i] - x2[i]) / dt;
/* store speed */
for (i = 0; i <= 2; i++) {
xp[i+3] = dx[i];
}
}
if (xpret != NULL)
for (i = 0; i <= 5; i++)
xpret[i] = xp[i];
}
return(OK);
}
/* Prepare lookup table of sin and cos ( i*Lj )
* for required multiple angles
*/
static void sscc (int k, double arg, int n)
{
double cu, su, cv, sv, s;
int i;
su = sin (arg);
cu = cos (arg);
ss[k][0] = su; /* sin(L) */
cc[k][0] = cu; /* cos(L) */
sv = 2.0 * su * cu;
cv = cu * cu - su * su;
ss[k][1] = sv; /* sin(2L) */
cc[k][1] = cv;
for (i = 2; i < n; i++)
{
s = su * cv + cu * sv;
cv = cu * cv - su * sv;
sv = s;
ss[k][i] = sv; /* sin( i+1 L ) */
cc[k][i] = cv;
}
}
/* Adjust position from Earth-Moon barycenter to Earth
*
* J = Julian day number
* xemb = rectangular equatorial coordinates of Earth
*/
static void embofs_mosh(double tjd, double *xemb)
{
double T, M, a, L, B, p;
double smp, cmp, s2mp, c2mp, s2d, c2d, sf, cf;
double s2f, sx, cx, xyz[6];
double seps = swed.oec.seps;
double ceps = swed.oec.ceps;
int i;
/* Short series for position of the Moon
*/
T = (tjd-J1900)/36525.0;
/* Mean anomaly of moon (MP) */
a = swe_degnorm(((1.44e-5*T + 0.009192)*T + 477198.8491)*T + 296.104608);
a *= DEGTORAD;
smp = sin(a);
cmp = cos(a);
s2mp = 2.0*smp*cmp; /* sin(2MP) */
c2mp = cmp*cmp - smp*smp; /* cos(2MP) */
/* Mean elongation of moon (D) */
a = swe_degnorm(((1.9e-6*T - 0.001436)*T + 445267.1142)*T + 350.737486);
a = 2.0 * DEGTORAD * a;
s2d = sin(a);
c2d = cos(a);
/* Mean distance of moon from its ascending node (F) */
a = swe_degnorm((( -3.e-7*T - 0.003211)*T + 483202.0251)*T + 11.250889);
a *= DEGTORAD;
sf = sin(a);
cf = cos(a);
s2f = 2.0*sf*cf; /* sin(2F) */
sx = s2d*cmp - c2d*smp; /* sin(2D - MP) */
cx = c2d*cmp + s2d*smp; /* cos(2D - MP) */
/* Mean longitude of moon (LP) */
L = ((1.9e-6*T - 0.001133)*T + 481267.8831)*T + 270.434164;
/* Mean anomaly of sun (M) */
M = swe_degnorm((( -3.3e-6*T - 1.50e-4)*T + 35999.0498)*T + 358.475833);
/* Ecliptic longitude of the moon */
L = L
+ 6.288750*smp
+ 1.274018*sx
+ 0.658309*s2d
+ 0.213616*s2mp
- 0.185596*sin( DEGTORAD * M )
- 0.114336*s2f;
/* Ecliptic latitude of the moon */
a = smp*cf;
sx = cmp*sf;
B = 5.128189*sf
+ 0.280606*(a+sx) /* sin(MP+F) */
+ 0.277693*(a-sx) /* sin(MP-F) */
+ 0.173238*(s2d*cf - c2d*sf); /* sin(2D-F) */
B *= DEGTORAD;
/* Parallax of the moon */
p = 0.950724
+0.051818*cmp
+0.009531*cx
+0.007843*c2d
+0.002824*c2mp;
p *= DEGTORAD;
/* Elongation of Moon from Sun
*/
L = swe_degnorm(L);
L *= DEGTORAD;
/* Distance in au */
a = 4.263523e-5/sin(p);
/* Convert to rectangular ecliptic coordinates */
xyz[0] = L;
xyz[1] = B;
xyz[2] = a;
swi_polcart(xyz, xyz);
/* Convert to equatorial */
swi_coortrf2(xyz, xyz, -seps, ceps);
/* Precess to equinox of J2000.0 */
swi_precess(xyz, tjd, 0, J_TO_J2000);/**/
/* now emb -> earth */
for (i = 0; i <= 2; i++)
xemb[i] -= xyz[i] / (EARTH_MOON_MRAT + 1.0);
}
/* orbital elements of planets that are computed from osculating elements
* epoch
* equinox
* mean anomaly,
* semi axis,
* eccentricity,
* argument of perihelion,
* ascending node
* inclination
*/
#define SE_NEELY /* use James Neely's revised elements
* of Uranian planets*/
static char *plan_fict_nam[SE_NFICT_ELEM] =
{"Cupido", "Hades", "Zeus", "Kronos",
"Apollon", "Admetos", "Vulkanus", "Poseidon",
"Isis-Transpluto", "Nibiru", "Harrington",
"Leverrier", "Adams",
"Lowell", "Pickering",};
char *swi_get_fict_name(int32 ipl, char *snam)
{
if (read_elements_file(ipl, 0, NULL, NULL,
NULL, NULL, NULL, NULL, NULL, NULL,
snam, NULL, NULL) == ERR)
strcpy(snam, "name not found");
return snam;
}
static double plan_oscu_elem[SE_NFICT_ELEM][8] = {
#ifdef SE_NEELY
{J1900, J1900, 163.7409, 40.99837, 0.00460, 171.4333, 129.8325, 1.0833},/* Cupido Neely */
{J1900, J1900, 27.6496, 50.66744, 0.00245, 148.1796, 161.3339, 1.0500},/* Hades Neely */
{J1900, J1900, 165.1232, 59.21436, 0.00120, 299.0440, 0.0000, 0.0000},/* Zeus Neely */
{J1900, J1900, 169.0193, 64.81960, 0.00305, 208.8801, 0.0000, 0.0000},/* Kronos Neely */
{J1900, J1900, 138.0533, 70.29949, 0.00000, 0.0000, 0.0000, 0.0000},/* Apollon Neely */
{J1900, J1900, 351.3350, 73.62765, 0.00000, 0.0000, 0.0000, 0.0000},/* Admetos Neely */
{J1900, J1900, 55.8983, 77.25568, 0.00000, 0.0000, 0.0000, 0.0000},/* Vulcanus Neely */
{J1900, J1900, 165.5163, 83.66907, 0.00000, 0.0000, 0.0000, 0.0000},/* Poseidon Neely */
#else
{J1900, J1900, 104.5959, 40.99837, 0, 0, 0, 0}, /* Cupido */
{J1900, J1900, 337.4517, 50.667443, 0, 0, 0, 0}, /* Hades */
{J1900, J1900, 104.0904, 59.214362, 0, 0, 0, 0}, /* Zeus */
{J1900, J1900, 17.7346, 64.816896, 0, 0, 0, 0}, /* Kronos */
{J1900, J1900, 138.0354, 70.361652, 0, 0, 0, 0}, /* Apollon */
{J1900, J1900, -8.678, 73.736476, 0, 0, 0, 0}, /* Admetos */
{J1900, J1900, 55.9826, 77.445895, 0, 0, 0, 0}, /* Vulkanus */
{J1900, J1900, 165.3595, 83.493733, 0, 0, 0, 0}, /* Poseidon */
#endif
/* Isis-Transpluto; elements from "Die Sterne" 3/1952, p. 70ff.
* Strubell does not give an equinox. 1945 is taken to best reproduce
* ASTRON ephemeris. (This is a strange choice, though.)
* The epoch is 1772.76. The year is understood to have 366 days.
* The fraction is counted from 1 Jan. 1772 */
{2368547.66, 2431456.5, 0.0, 77.775, 0.3, 0.7, 0, 0},
/* Nibiru, elements from Christian Woeltge, Hannover */
{1856113.380954, 1856113.380954, 0.0, 234.8921, 0.981092, 103.966, -44.567, 158.708},
/* Harrington, elements from Astronomical Journal 96(4), Oct. 1988 */
{2374696.5, J2000, 0.0, 101.2, 0.411, 208.5, 275.4, 32.4},
/* Leverrier's Neptune,
according to W.G. Hoyt, "Planets X and Pluto", Tucson 1980, p. 63 */
{2395662.5, 2395662.5, 34.05, 36.15, 0.10761, 284.75, 0, 0},
/* Adam's Neptune */
{2395662.5, 2395662.5, 24.28, 37.25, 0.12062, 299.11, 0, 0},
/* Lowell's Pluto */
{2425977.5, 2425977.5, 281, 43.0, 0.202, 204.9, 0, 0},
/* Pickering's Pluto */
{2425977.5, 2425977.5, 48.95, 55.1, 0.31, 280.1, 100, 15}, /**/
#if 0 /* Ceres JPL 1600, without perturbations from other minor planets,
* from following initial elements:
* 2450600.5 2000 0 1 164.7073602 73.0340746 80.5995101
* 10.5840296 0.07652422 0.0 2.770176095 */
{2305447.5, J2000, 0.5874558977449977e+02, 0.2766536058742327e+01,
0.7870946565779195e-01, 0.5809199028919189e+02,
0.8650119410725021e+02, 0.1066835622280712e+02},
/* Chiron, Bowell database 18-mar-1997 */
{2450500.5, J2000, 7.258191, 13.67387471, 0.38174778, 339.558345, 209.379239, 6.933360}, /**/
#endif
};
/* computes a planet from osculating elements *
* tjd julian day
* ipl body number
* ipli body number in planetary data structure
* iflag flags
*/
int swi_osc_el_plan(double tjd, double *xp, int ipl, int ipli, double *xearth, double *xsun, char *serr)
{
double pqr[9], x[6];
double eps, K, fac, rho, cose, sine;
double alpha, beta, zeta, sigma, M2, Msgn, M_180_or_0;
double tjd0, tequ, mano, sema, ecce, parg, node, incl, dmot;
double cosnode, sinnode, cosincl, sinincl, cosparg, sinparg;
double M, E;
struct plan_data *pedp = &swed.pldat[SEI_EARTH];
struct plan_data *pdp = &swed.pldat[ipli];
int32 fict_ifl = 0;
int i;
/* orbital elements, either from file or, if file not found,
* from above built-in set
*/
if (read_elements_file(ipl, tjd, &tjd0, &tequ,
&mano, &sema, &ecce, &parg, &node, &incl,
NULL, &fict_ifl, serr) == ERR)
return ERR;
dmot = 0.9856076686 * DEGTORAD / sema / sqrt(sema); /* daily motion */
if (fict_ifl & FICT_GEO)
dmot /= sqrt(SUN_EARTH_MRAT);
cosnode = cos(node);
sinnode = sin(node);
cosincl = cos(incl);
sinincl = sin(incl);
cosparg = cos(parg);
sinparg = sin(parg);
/* Gaussian vector */
pqr[0] = cosparg * cosnode - sinparg * cosincl * sinnode;
pqr[1] = -sinparg * cosnode - cosparg * cosincl * sinnode;
pqr[2] = sinincl * sinnode;
pqr[3] = cosparg * sinnode + sinparg * cosincl * cosnode;
pqr[4] = -sinparg * sinnode + cosparg * cosincl * cosnode;
pqr[5] = -sinincl * cosnode;
pqr[6] = sinparg * sinincl;
pqr[7] = cosparg * sinincl;
pqr[8] = cosincl;
/* Kepler problem */
E = M = swi_mod2PI(mano + (tjd - tjd0) * dmot); /* mean anomaly of date */
/* better E for very high eccentricity and small M */
if (ecce > 0.975) {
M2 = M * RADTODEG;
if (M2 > 150 && M2 < 210) {
M2 -= 180;
M_180_or_0 = 180;
} else
M_180_or_0 = 0;
if (M2 > 330)
M2 -= 360;
if (M2 < 0) {
M2 = -M2;
Msgn = -1;
} else
Msgn = 1;
if (M2 < 30) {
M2 *= DEGTORAD;
alpha = (1 - ecce) / (4 * ecce + 0.5);
beta = M2 / (8 * ecce + 1);
zeta = pow(beta + sqrt(beta * beta + alpha * alpha), 1/3);
sigma = zeta - alpha / 2;
sigma = sigma - 0.078 * sigma * sigma * sigma * sigma * sigma / (1 + ecce);
E = Msgn * (M2 + ecce * (3 * sigma - 4 * sigma * sigma * sigma))
+ M_180_or_0;
}
}
E = swi_kepler(E, M, ecce);
/* position and speed, referred to orbital plane */
if (fict_ifl & FICT_GEO)
K = KGAUSS_GEO / sqrt(sema);
else
K = KGAUSS / sqrt(sema);
cose = cos(E);
sine = sin(E);
fac = sqrt((1 - ecce) * (1 + ecce));
rho = 1 - ecce * cose;
x[0] = sema * (cose - ecce);
x[1] = sema * fac * sine;
x[3] = -K * sine / rho;
x[4] = K * fac * cose / rho;
/* transformation to ecliptic */
xp[0] = pqr[0] * x[0] + pqr[1] * x[1];
xp[1] = pqr[3] * x[0] + pqr[4] * x[1];
xp[2] = pqr[6] * x[0] + pqr[7] * x[1];
xp[3] = pqr[0] * x[3] + pqr[1] * x[4];
xp[4] = pqr[3] * x[3] + pqr[4] * x[4];
xp[5] = pqr[6] * x[3] + pqr[7] * x[4];
/* transformation to equator */
eps = swi_epsiln(tequ, 0);
swi_coortrf(xp, xp, -eps);
swi_coortrf(xp+3, xp+3, -eps);
/* precess to J2000 */
if (tequ != J2000) {
swi_precess(xp, tequ, 0, J_TO_J2000);
swi_precess(xp+3, tequ, 0, J_TO_J2000);
}
/* to solar system barycentre */
if (fict_ifl & FICT_GEO) {
for (i = 0; i <= 5; i++) {
xp[i] += xearth[i];
}
} else {
for (i = 0; i <= 5; i++) {
xp[i] += xsun[i];
}
}
if (pdp->x == xp) {
pdp->teval = tjd; /* for precession! */
pdp->iephe = pedp->iephe;
}
return OK;
}
#if 1
/* note: input parameter tjd is required for T terms in elements */
static int read_elements_file(int32 ipl, double tjd,
double *tjd0, double *tequ,
double *mano, double *sema, double *ecce,
double *parg, double *node, double *incl,
char *pname, int32 *fict_ifl, char *serr)
{
int i, iline, iplan, retc, ncpos;
FILE *fp = NULL;
char s[AS_MAXCH], *sp;
char *cpos[20], serri[AS_MAXCH];
AS_BOOL elem_found = FALSE;
double tt = 0;
/* -1, because file information is not saved, file is always closed */
if ((fp = swi_fopen(-1, SE_FICTFILE, swed.ephepath, serr)) == NULL) {
/* file does not exist, use built-in bodies */
if (ipl >= SE_NFICT_ELEM) {
if (serr != NULL)
sprintf(serr, "error no elements for fictitious body no %7.0f", (double) ipl);
return ERR;
}
if (tjd0 != NULL)
*tjd0 = plan_oscu_elem[ipl][0]; /* epoch */
if (tequ != NULL)
*tequ = plan_oscu_elem[ipl][1]; /* equinox */
if (mano != NULL)
*mano = plan_oscu_elem[ipl][2] * DEGTORAD; /* mean anomaly */
if (sema != NULL)
*sema = plan_oscu_elem[ipl][3]; /* semi-axis */
if (ecce != NULL)
*ecce = plan_oscu_elem[ipl][4]; /* eccentricity */
if (parg != NULL)
*parg = plan_oscu_elem[ipl][5] * DEGTORAD; /* arg. of peri. */
if (node != NULL)
*node = plan_oscu_elem[ipl][6] * DEGTORAD; /* asc. node */
if (incl != NULL)
*incl = plan_oscu_elem[ipl][7] * DEGTORAD; /* inclination */
if (pname != NULL)
strcpy(pname, plan_fict_nam[ipl]);
return OK;
}
/*
* find elements in file
*/
iline = 0;
iplan = -1;
while (fgets(s, AS_MAXCH, fp) != NULL) {
iline++;
sp = s;
while(*sp == ' ' || *sp == '\t')
sp++;
swi_strcpy(s, sp);
if (*s == '#')
continue;
if (*s == '\r')
continue;
if (*s == '\n')
continue;
if (*s == '\0')
continue;
if ((sp = strchr(s, '#')) != NULL)
*sp = '\0';
ncpos = swi_cutstr(s, ",", cpos, 20);
sprintf(serri, "error in file %s, line %7.0f:", SE_FICTFILE, (double) iline);
if (ncpos < 9) {
if (serr != NULL) {
sprintf(serr, "%s nine elements required", serri);
}
goto return_err;
}
iplan++;
if (iplan != ipl)
continue;
elem_found = TRUE;
/* epoch of elements */
if (tjd0 != NULL) {
sp = cpos[0];
for (i = 0; i < 5; i++)
sp[i] = tolower(sp[i]);
if (strncmp(sp, "j2000", 5) == OK)
*tjd0 = J2000;
else if (strncmp(sp, "b1950", 5) == OK)
*tjd0 = B1950;
else if (strncmp(sp, "j1900", 5) == OK)
*tjd0 = J1900;
else if (*sp == 'j' || *sp == 'b') {
if (serr != NULL) {
sprintf(serr, "%s invalid epoch", serri);
}
goto return_err;
} else
*tjd0 = atof(sp);
tt = tjd - *tjd0;
}
/* equinox */
if (tequ != NULL) {
sp = cpos[1];
while(*sp == ' ' || *sp == '\t')
sp++;
for (i = 0; i < 5; i++)
sp[i] = tolower(sp[i]);
if (strncmp(sp, "j2000", 5) == OK)
*tequ = J2000;
else if (strncmp(sp, "b1950", 5) == OK)
*tequ = B1950;
else if (strncmp(sp, "j1900", 5) == OK)
*tequ = J1900;
else if (strncmp(sp, "jdate", 5) == OK)
*tequ = tjd;
else if (*sp == 'j' || *sp == 'b') {
if (serr != NULL) {
sprintf(serr, "%s invalid equinox", serri);
}
goto return_err;
} else
*tequ = atof(sp);
}
/* mean anomaly t0 */
if (mano != NULL) {
retc = check_t_terms(tt, cpos[2], mano);
*mano = swe_degnorm(*mano);
if (retc == ERR) {
if (serr != NULL) {
sprintf(serr, "%s mean anomaly value invalid", serri);
}
goto return_err;
}
/* if mean anomaly has t terms (which happens with fictitious
* planet Vulcan), we set
* epoch = tjd, so that no motion will be added anymore
* equinox = tjd */
if (retc == 1) {
*tjd0 = tjd;
}
*mano *= DEGTORAD;
}
/* semi-axis */
if (sema != NULL) {
retc = check_t_terms(tt, cpos[3], sema);
if (*sema <= 0 || retc == ERR) {
if (serr != NULL) {
sprintf(serr, "%s semi-axis value invalid", serri);
}
goto return_err;
}
}
/* eccentricity */
if (ecce != NULL) {
retc = check_t_terms(tt, cpos[4], ecce);
if (*ecce >= 1 || *ecce < 0 || retc == ERR) {
if (serr != NULL) {
sprintf(serr, "%s eccentricity invalid (no parabolic or hyperbolic orbits allowed)", serri);
}
goto return_err;
}
}
/* perihelion argument */
if (parg != NULL) {
retc = check_t_terms(tt, cpos[5], parg);
*parg = swe_degnorm(*parg);
if (retc == ERR) {
if (serr != NULL) {
sprintf(serr, "%s perihelion argument value invalid", serri);
}
goto return_err;
}
*parg *= DEGTORAD;
}
/* node */
if (node != NULL) {
retc = check_t_terms(tt, cpos[6], node);
*node = swe_degnorm(*node);
if (retc == ERR) {
if (serr != NULL) {
sprintf(serr, "%s node value invalid", serri);
}
goto return_err;
}
*node *= DEGTORAD;
}
/* inclination */
if (incl != NULL) {
retc = check_t_terms(tt, cpos[7], incl);
*incl = swe_degnorm(*incl);
if (retc == ERR) {
if (serr != NULL) {
sprintf(serr, "%s inclination value invalid", serri);
}
goto return_err;
}
*incl *= DEGTORAD;
}
/* planet name */
if (pname != NULL) {
sp = cpos[8];
while(*sp == ' ' || *sp == '\t')
sp++;
swi_right_trim(sp);
strcpy(pname, sp);
}
/* geocentric */
if (fict_ifl != NULL && ncpos > 9) {
for (sp = cpos[9]; *sp != '\0'; sp++)
*sp = tolower(*sp);
if (strstr(cpos[9], "geo") != NULL)
*fict_ifl |= FICT_GEO;
}
break;
}
if (!elem_found) {
if (serr != NULL) {
sprintf(serr, "%s elements for planet %7.0f not found", serri, (double) ipl);
}
goto return_err;
}
fclose(fp);
return OK;
return_err:
fclose(fp);
return ERR;
}
#endif
static int check_t_terms(double t, char *sinp, double *doutp)
{
int i, isgn = 1, z;
int retc = 0;
char *sp;
double tt[5], fac;
tt[0] = t / 36525;
tt[1] = tt[0];
tt[2] = tt[1] * tt[1];
tt[3] = tt[2] * tt[1];
tt[4] = tt[3] * tt[1];
if ((sp = strpbrk(sinp, "+-")) != NULL)
retc = 1; /* with additional terms */
sp = sinp;
*doutp = 0;
fac = 1;
z = 0;
while (1) {
while(*sp != '\0' && strchr(" \t", *sp) != NULL)
sp++;
if (strchr("+-", *sp) || *sp == '\0') {
if (z > 0)
*doutp += fac;
isgn = 1;
if (*sp == '-')
isgn = -1;
fac = 1 * isgn;
if (*sp == '\0')
return retc;
sp++;
} else {
while(*sp != '\0' && strchr("* \t", *sp) != NULL)
sp++;
if (*sp != '\0' && strchr("tT", *sp) != NULL) {
/* a T */
sp++;
if (*sp != '\0' && strchr("+-", *sp))
fac *= tt[0];
else if ((i = atoi(sp)) <= 4 && i >= 0)
fac *= tt[i];
} else {
/* a number */
if (atof(sp) != 0 || *sp == '0')
fac *= atof(sp);
}
while (*sp != '\0' && strchr("0123456789.", *sp))
sp++;
}
z++;
}
return retc; /* there have been additional terms */
}
Reference in New Issue View Git Blame Copy Permalink