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
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
|
SUBROUTINE slERTH (IY, ID, FD, PV)
*+
* - - - - - -
* E R T H
* - - - - - -
*
* Approximate heliocentric position and velocity of the Earth
*
* Given:
* IY I year
* ID I day in year (1 = Jan 1st)
* FD R fraction of day
*
* Returned:
* PV R(6) Earth position & velocity vector
*
* Notes:
*
* 1 The date and time is TDB (loosely ET) in a Julian calendar
* which has been aligned to the ordinary Gregorian
* calendar for the interval 1900 March 1 to 2100 February 28.
* The year and day can be obtained by calling slCAYD or
* slCLYD.
*
* 2 The Earth heliocentric 6-vector is mean equator and equinox
* of date. Position part, PV(1-3), is in AU; velocity part,
* PV(4-6), is in AU/sec.
*
* 3 Max/RMS errors 1950-2050:
* 13/5 E-5 AU = 19200/7600 km in position
* 47/26 E-10 AU/s = 0.0070/0.0039 km/s in speed
*
* 4 More accurate results are obtainable with the routines slEVP
* and slEPV.
*
* Last revision: 5 April 2005
*
* Copyright P.T.Wallace. All rights reserved.
*
* License:
* 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 (see SLA_CONDITIONS); if not, write to the
* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
* Boston, MA 02110-1301 USA
*
* Copyright (C) 1995 Association of Universities for Research in Astronomy Inc.
*-
IMPLICIT NONE
INTEGER IY,ID
REAL FD,PV(6)
INTEGER IY4
REAL TWOPI,SPEED,REMB,SEMB,YI,YF,T,ELM,GAMMA,EM,ELT,EPS0,
: E,ESQ,V,R,ELMM,COSELT,SINEPS,COSEPS,W1,W2,SELMM,CELMM
PARAMETER (TWOPI=6.28318530718)
* Mean orbital speed of Earth, AU/s
PARAMETER (SPEED=1.9913E-7)
* Mean Earth:EMB distance and speed, AU and AU/s
PARAMETER (REMB=3.12E-5,SEMB=8.31E-11)
* Whole years & fraction of year, and years since J1900.0
YI=FLOAT(IY-1900)
IY4=MOD(MOD(IY,4)+4,4)
YF=(FLOAT(4*(ID-1/(IY4+1))-IY4-2)+4.0*FD)/1461.0
T=YI+YF
* Geometric mean longitude of Sun
* (cf 4.881627938+6.283319509911*T MOD 2PI)
ELM=MOD(4.881628+TWOPI*YF+0.00013420*T,TWOPI)
* Mean longitude of perihelion
GAMMA=4.908230+3.0005E-4*T
* Mean anomaly
EM=ELM-GAMMA
* Mean obliquity
EPS0=0.40931975-2.27E-6*T
* Eccentricity
E=0.016751-4.2E-7*T
ESQ=E*E
* True anomaly
V=EM+2.0*E*SIN(EM)+1.25*ESQ*SIN(2.0*EM)
* True ecliptic longitude
ELT=V+GAMMA
* True distance
R=(1.0-ESQ)/(1.0+E*COS(V))
* Moon's mean longitude
ELMM=MOD(4.72+83.9971*T,TWOPI)
* Useful functions
COSELT=COS(ELT)
SINEPS=SIN(EPS0)
COSEPS=COS(EPS0)
W1=-R*SIN(ELT)
W2=-SPEED*(COSELT+E*COS(GAMMA))
SELMM=SIN(ELMM)
CELMM=COS(ELMM)
* Earth position and velocity
PV(1)=-R*COSELT-REMB*CELMM
PV(2)=(W1-REMB*SELMM)*COSEPS
PV(3)=W1*SINEPS
PV(4)=SPEED*(SIN(ELT)+E*SIN(GAMMA))+SEMB*SELMM
PV(5)=(W2-SEMB*CELMM)*COSEPS
PV(6)=W2*SINEPS
END
|
