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
|
/* getplane.cpp: functions for computing planetary ephems
Copyright (C) 2010, Project Pluto
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
02110-1301, USA. */
#include <math.h>
#include <string.h>
#include "watdefs.h"
#include "lunar.h"
#include "afuncs.h"
#define PI 3.1415926535897932384626433832795028841971693993751058209749445923
/* The compute_planet( ) function is supposed to provide a "standardized"
way to get positions for planets and the moon, in any of four systems.
You give it a pointer to the VSOP data, a planet number (0=sun,
1=mercury, ... 9 = Pluto, 10 = the Moon), and a time in centuries
from J2000. The returned array of fifteen doubles gives you the object
position in five different systems:
ovals[0, 1, 2] = lon, lat, r, heliocentric ecliptic of date;
ovals[3, 4, 5] = x, y, z, heliocentric ecliptic of date;
ovals[6, 7, 8] = x, y, z, heliocentric equatorial of date;
ovals[9, 10, 11] = x, y, z, heliocentric equatorial, J2000.0
ovals[12, 13, 14] = x, y, z, heliocentric ecliptical, J2000.0
I intended to rig this up so you could keep on going to 11=Io, 12=Europa,
etc. This would make all sorts of sense. But I've not done it yet. */
int DLL_FUNC compute_planet( const char FAR *vsop_data, const int planet_no,
const double t_c, double DLLPTR *ovals)
{
double lat, lon, r;
const double obliquit = mean_obliquity( t_c);
double matrix[9];
const double obliq_2000 = 23.4392911 * PI / 180.;
/* first, compute polar heliocentric in eclip of date */
if( planet_no != 10)
{
lon = calc_vsop_loc( vsop_data, planet_no, 0, t_c, 0.);
lat = calc_vsop_loc( vsop_data, planet_no, 1, t_c, 0.);
r = calc_vsop_loc( vsop_data, planet_no, 2, t_c, 0.);
}
else
{
double fund[N_FUND];
lunar_fundamentals( vsop_data, t_c, fund);
lat = lunar_lat( vsop_data, fund, 0L);
lunar_lon_and_dist( vsop_data, fund, &lon, &r, 0L);
lon *= PI / 180.;
lat *= PI / 180.;
r /= AU_IN_KM; /* from km to AU */
}
ovals[0] = lon;
ovals[1] = lat;
ovals[2] = r;
/* next, compute polar cartesian in eclip of date */
ovals[3] = cos( lon) * cos( lat) * r;
ovals[4] = sin( lon) * cos( lat) * r;
ovals[5] = sin( lat) * r;
/* next, compute polar cartesian in eclip of date, */
/* but in equatorial coords */
FMEMCPY( ovals + 6, ovals + 3, 3 * sizeof( double));
rotate_vector( ovals + 6, obliquit, 0);
/* next, precess to get J2000.0 equatorial values */
setup_precession( matrix, 2000. + t_c * 100., 2000.);
precess_vector( matrix, ovals + 6, ovals + 9);
/* Finally, rotate equatorial J2000.0 into ecliptical J2000 */
FMEMCPY( ovals + 12, ovals + 9, 3 * sizeof( double));
rotate_vector( ovals + 12, -obliq_2000, 0);
return( 0);
}
|
