1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
|
#include <math.h>
#include "observe.h"
/* Assorted functions useful in conjunction with the satellite code
library for determining where the _observer_, as well as the _target_,
happens to be. Combine the two positions, and you can get the
distance/RA/dec of the target as seen by the observer. */
#define PI 3.141592653589793238462643383279
#define EARTH_MAJOR_AXIS 6378140.
#define EARTH_MINOR_AXIS 6356755.
#define EARTH_AXIS_RATIO (EARTH_MINOR_AXIS / EARTH_MAJOR_AXIS)
/* function for Greenwich sidereal time, ripped from 'deep.cpp' */
static inline double ThetaG( double jd)
{
/* Reference: The 1992 Astronomical Almanac, page B6. */
const double omega_E = 1.00273790934;
/* Earth rotations per sidereal day (non-constant) */
const double UT = fmod( jd + .5, 1.);
const double seconds_per_day = 86400.;
const double jd_2000 = 2451545.0; /* 1.5 Jan 2000 = JD 2451545. */
double t_cen, GMST, rval;
t_cen = (jd - UT - jd_2000) / 36525.;
GMST = 24110.54841 + t_cen * (8640184.812866 + t_cen *
(0.093104 - t_cen * 6.2E-6));
GMST = fmod( GMST + seconds_per_day * omega_E * UT, seconds_per_day);
if( GMST < 0.)
GMST += seconds_per_day;
rval = 2. * PI * GMST / seconds_per_day;
return( rval);
} /*Function thetag*/
void DLL_FUNC observer_cartesian_coords( const double jd, const double lon,
const double rho_cos_phi, const double rho_sin_phi,
double *vect)
{
const double angle = lon + ThetaG( jd);
*vect++ = cos( angle) * rho_cos_phi * EARTH_MAJOR_AXIS / 1000.;
*vect++ = sin( angle) * rho_cos_phi * EARTH_MAJOR_AXIS / 1000.;
*vect++ = rho_sin_phi * EARTH_MAJOR_AXIS / 1000.;
}
void DLL_FUNC lat_alt_to_parallax( const double lat, const double ht_in_meters,
double *rho_cos_phi, double *rho_sin_phi)
{
const double u = atan( sin( lat) * EARTH_AXIS_RATIO / cos( lat));
*rho_sin_phi = EARTH_AXIS_RATIO * sin( u) +
(ht_in_meters / EARTH_MAJOR_AXIS) * sin( lat);
*rho_cos_phi = cos( u) + (ht_in_meters / EARTH_MAJOR_AXIS) * cos( lat);
}
void DLL_FUNC get_satellite_ra_dec_delta( const double *observer_loc,
const double *satellite_loc, double *ra,
double *dec, double *delta)
{
double vect[3], dist2 = 0.;
int i;
for( i = 0; i < 3; i++)
{
vect[i] = satellite_loc[i] - observer_loc[i];
dist2 += vect[i] * vect[i];
}
*delta = sqrt( dist2);
*ra = atan2( vect[1], vect[0]);
if( *ra < 0.)
*ra += PI + PI;
*dec = asin( vect[2] / *delta);
}
/* Formulae from Meeus' _Astronomical Algorithms_ for approximate precession.
More than accurate enough for our purposes. */
static void precess( const double t_centuries, double *ra, double *dec)
{
const double m = (3.07496 + .00186 * t_centuries / 2.) * (PI / 180.) / 240.;
const double n = (1.33621 - .00057 * t_centuries / 2.) * (PI / 180.) / 240.;
const double ra_rate = m + n * sin( *ra) * tan( *dec);
const double dec_rate = n * cos( *ra);
*ra -= t_centuries * ra_rate * 100.;
*dec -= t_centuries * dec_rate * 100.;
}
void DLL_FUNC epoch_of_date_to_j2000( const double jd, double *ra, double *dec)
{
const double t_centuries = (jd - 2451545.) / 36525.;
precess( t_centuries, ra, dec);
}
void DLL_FUNC j2000_to_epoch_of_date( const double jd, double *ra, double *dec)
{
const double t_centuries = (jd - 2451545.) / 36525.;
precess( -t_centuries, ra, dec);
}
|
