mirror of https://github.com/XEphem/XEphem.git
71 lines
1.9 KiB
C
71 lines
1.9 KiB
C
#include <stdio.h>
|
|
#include <math.h>
|
|
|
|
#include "astro.h"
|
|
|
|
static void ecleq_aux (int sw, double mj, double x, double y,
|
|
double *p, double *q);
|
|
|
|
#define EQtoECL 1
|
|
#define ECLtoEQ (-1)
|
|
|
|
|
|
/* given the modified Julian date, mj, and an equitorial ra and dec, each in
|
|
* radians, find the corresponding geocentric ecliptic latitude, *lt, and
|
|
* longititude, *lg, also each in radians.
|
|
* correction for the effect on the angle of the obliquity due to nutation is
|
|
* not included.
|
|
*/
|
|
void
|
|
eq_ecl (double mj, double ra, double dec, double *lt, double *lg)
|
|
{
|
|
ecleq_aux (EQtoECL, mj, ra, dec, lg, lt);
|
|
}
|
|
|
|
/* given the modified Julian date, mj, and a geocentric ecliptic latitude,
|
|
* *lt, and longititude, *lg, each in radians, find the corresponding
|
|
* equitorial ra and dec, also each in radians.
|
|
* correction for the effect on the angle of the obliquity due to nutation is
|
|
* not included.
|
|
*/
|
|
void
|
|
ecl_eq (double mj, double lt, double lg, double *ra, double *dec)
|
|
{
|
|
ecleq_aux (ECLtoEQ, mj, lg, lt, ra, dec);
|
|
}
|
|
|
|
static void
|
|
ecleq_aux (
|
|
int sw, /* +1 for eq to ecliptic, -1 for vv. */
|
|
double mj,
|
|
double x, double y, /* sw==1: x==ra, y==dec. sw==-1: x==lg, y==lt. */
|
|
double *p, double *q) /* sw==1: p==lg, q==lt. sw==-1: p==ra, q==dec. */
|
|
{
|
|
static double lastmj = -10000; /* last mj calculated */
|
|
static double seps, ceps; /* sin and cos of mean obliquity */
|
|
double sx, cx, sy, cy, ty, sq;
|
|
|
|
if (mj != lastmj) {
|
|
double eps;
|
|
obliquity (mj, &eps); /* mean obliquity for date */
|
|
seps = sin(eps);
|
|
ceps = cos(eps);
|
|
lastmj = mj;
|
|
}
|
|
|
|
sy = sin(y);
|
|
cy = cos(y); /* always non-negative */
|
|
if (fabs(cy)<1e-20) cy = 1e-20; /* insure > 0 */
|
|
ty = sy/cy;
|
|
cx = cos(x);
|
|
sx = sin(x);
|
|
sq = (sy*ceps)-(cy*seps*sx*sw);
|
|
if (sq < -1) sq = -1;
|
|
if (sq > 1) sq = 1;
|
|
*q = asin(sq);
|
|
*p = atan(((sx*ceps)+(ty*seps*sw))/cx);
|
|
if (cx<0) *p += PI; /* account for atan quad ambiguity */
|
|
range (p, 2*PI);
|
|
}
|
|
|