 |
From: http://hubble.lamost.org/cpp.htmz1M
Computing planetary positions - a tutorial with worked examples>w:k!
By Paul Schlyter, Stockholm, Sweden^4YdS^
email: pausch@saaf.se or paul.schlyter@ausys.sez"iL$
WWW: http://welcome.to/pausch or http://hotel04.ausys.se/pausch\^fK
©緣生術數研究社 -- 術數研究 N
Break out of a frame'L#ae
©緣生術數研究社 -- 術數研究 `oxGjS
1. Fundamentals
2. Some useful functions
3. Rectangular and spherical coordinates
4. The time scale. A test date
5. The Sun's position
6. Sidereal time and hour angle. Altitude and azimuth
7. The Moon's position
8. The Moon's position with higher accuracy. Perturbations
9. The Moon's topocentric position
10. The orbital elements of the planets
11. The heliocentric positions of the planets
12. Higher accuracy - perturbations
13. Precession
14. Geocentric positions of the planets
15. The elongation and physical ephemerides of the planets
16. The positions of comets. Comet Encke and Levy
Today it's not really that difficult to compute a planet's position from its orbital elements. The only thing you'll need is a computer and a suitable program. If you want to write such a program yourself, this text contains the formulae you need. The aim here is to obtain the planetary positions at any date during the 20'th and 21'st century with an error of one or at the most two arc minutes, and to compute the position of an asteroid or a comet from its orbital elements./s
©緣生術數研究社 -- 術數研究 7%R
No programs are given, because different computers and calculators are programmed in different languages. It is easier to convert formulae to your favourite language than to translate a program from one programming language to another. Therefore formulae are presented instead. Each formula is accompanied with one numeric example - this will enable you to check your implementation of these formulae. Just remember that all numerical quantities contain rounding errors, therefore you may get slightly different results in your program compared to the numerical results presented here. When I checked these numerical values I used a HP-48SX pocket calculator, which uses 12 digits of accuracy.:N4Y
©緣生術數研究社 -- 術數研究 F
1. Fundamentalsa5G
A celestial body usually orbits the sun in an elliptical orbit. Perturbations from other planets causes small deviations from this elliptical orbit, but an unperturbed elliptical orbit can be used as a first approximation, and sometimes as the final approximation. If the distance from the Sun to the planet always is the same, then the planet follows a circular orbit. No planet does this, but the orbits of Venus and Neptune are very close to circles. Among the planets, Mercury and Pluto have orbits that deviate the most from a circle, i.e. are the most eccentric. Many asteroids have even more eccentric orbits, but the most eccentric orbits are to be found among the comets. Halley's comet, for instance, is closer to the Sun than Venus at perihelion, but farther away from the Sun than Neptune at aphelion. Some comets have even more eccentric orbits that are best approximated as a parabola. These orbits are not closed - a comet following a parabolic orbit passes the Sun only once, never to return. In reality these orbits are extremely elongated ellipses though, and these comets will eventually return, sometimes after many millennia.>~[
©緣生術數研究社 -- 術數研究 #FzWO
The perihelion and aphelion are the points in the orbit when the planet is closest to and most distant from the Sun. A parabolic orbit only has a perihemium of course.ucU
©緣生術數研究社 -- 術數研究 `=
The perigee and apogee are points in the Moon's orbit (or the orbit of an artificial Earth satellite) which are closest to and most distant from the Earth..V]
©緣生術數研究社 -- 術數研究 f=h
The celestial sphere is an imaginary sphere around the observer, at an arbitrary distance.?OOA
©緣生術數研究社 -- 術數研究 {:
The celestial equator is the Earth's equatorial plane projected on the celestial sphere.G
©緣生術數研究社 -- 術數研究 GENj
The ecliptic is the plane of the Earth's orbit. This is also the plane of the Sun's yearly apparent motion. The ecliptic is inclined by approximately 23.4_deg to the celestial equator. The ecliptic intersects the celestial equator at two points: The Vernal Point (or "the first point of Aries"), and the Autumnal Point. The Vernal Point is the point of origin for two different commonly used celestial coordinates: equatorial coordinates and ecliptic coordinates.544S
©緣生術數研究社 -- 術數研究 .F&l
Right Ascension and Declination are equatorial coordinates using the celestial equator as a fundamental plane. At the Vernal Point both the Right Ascension and the Declination are zero. The Right Ascension is usually measured in hours and minutes, where one revolution is 24 hours (which means 1 hour equals 15 degrees). It's counted countersunwise along the celestial equator. The Declination goes from +90 to -90 degrees, and it's positive north of, and negative south of, the celestial equator.9@]k7
©緣生術數研究社 -- 術數研究 %s
Longitude and Latitude are ecliptic coordinates, which use the ecliptic as a fundamental plane. Both are measured in degrees, and these coorinates too are both zero at the Vernal Point. The Longitude is counted countersunwise along the ecliptic. The Latitude is positive north of the ecliptic. Of course longitude and latitude are also used as terrestial coordinates, to measure a position of a point on the surface of the Earth.~:
©緣生術數研究社 -- 術數研究 u+$;Be
Heliocentric, Geocentric, Topocentric. A position relative to the Sun is heliocentric. If the position is relative to the center of the Earth, then it's geocentric. A topocentric position is relative to an observer on the surface of the Earth. Within the aim of our accuracy of 1-2 arc minutes, the difference between geocentric and topocentric position is negligible for all celestial bodies except the Moon (and some occasional asteroid which happens to pass very close to the Earth).{Rsd4U
©緣生術數研究社 -- 術數研究 `1/')D
The orbital elements consist of 6 quantities which completely define a circular, elliptic, parabolic or hyperbolic orbit. Three of these quantities describe the shape and size of the orbit, and the position of the planet in the orbit:jm^
©緣生術數研究社 -- 術數研究 P
a Mean distance, or semi-major axis8O1
e Eccentricitys
T Time at perihelionX7H
©緣生術數研究社 -- 術數研究 LGpO:
A cirular orbit has zero eccentricity. An elliptical orbit has an eccentricity between zero and one. A parabolic orbit has an eccentricity of exactly one. Finally, a hyperbolic orbit has an eccentricity larger than one. A parabolic orbit has an infinite semi-major axis, a, therefore one instead gives the perihelion distance, q, for a parabolic orbit:[var
©緣生術數研究社 -- 術數研究 dzH
q Perihelion distance = a * (1 - e)&
©緣生術數研究社 -- 術數研究 f([a
It is customary to give q instead of a also for hyperbolic orbit, and for elliptical orbits with eccentricity close to one.AaxC#G
©緣生術數研究社 -- 術數研究 !E`h
The three remaining orbital elements define the orientation of the orbit in space:NBof
©緣生術數研究社 -- 術數研究 qgVS7u
i Inclination, i.e. the "tilt" of the orbit relative to theKHkTG
ecliptic. The inclination varies from 0 to 180 degrees. If5y\6
the inclination is larger than 90 degrees, the planet is inc
a retrogade orbit, i.e. it moves "backwards". The mostTJ{T_H
well-known celestial body with retrogade motion is Comet Halley.^xD..N
©緣生術數研究社 -- 術數研究 I
N (usually written as "Capital Omega") Longitude of Ascending#<i%
Node. This is the angle, along the ecliptic, from the VernalQ
Point to the Ascending Node, which is the intersection betweenl&8a{<
the orbit and the ecliptic, where the planet moves from southTe
of to north of the ecliptic, i.e. from negative to positive8
latitudes.vr
©緣生術數研究社 -- 術數研究 ew
w (usually written as "small Omega") The angle from the Ascending3
node to the Perihelion, along the orbit.}T
©緣生術數研究社 -- 術數研究 >1p6
These are the primary orbital elements. From these many secondary orbital elements can be computed:tc~1r|
©緣生術數研究社 -- 術數研究 *R* <
q Perihelion distance = a * (1 - e)Zq
©緣生術數研究社 -- 術數研究 .xla
Q Aphelion distance = a * (1 + e)}*]XP
©緣生術數研究社 -- 術數研究 =)V`
P Orbital period = 365.256898326 * a**1.5/sqrt(1+m) days,x
where m = the mass of the planet in solar masses (0 for;]
comets and asteroids). sqrt() is the square root function.6k_
©緣生術數研究社 -- 術數研究 J$a
n Daily motion = 360_deg / P degrees/dayt
©緣生術數研究社 -- 術數研究 6@p[a*
t Some epoch as a day count, e.g. Julian Day Number. The TimeOnA+s.
at Perihelion, T, should then be expressed as the same day count.Ncw
©緣生術數研究社 -- 術數研究 0S$
t - T Time since Perihelion, usually in daysP:5!e?
©緣生術數研究社 -- 術數研究 .3[4{|
M Mean Anomaly = n * (t - T) = (t - T) * 360_deg / P#JAo
Mean Anomaly is 0 at perihelion and 180 degrees at aphelion5c'[4
©緣生術數研究社 -- 術數研究 ;%G5k
L Mean Longitude = M + w + NT'
©緣生術數研究社 -- 術數研究 <`
E Eccentric anomaly, defined by Kepler's equation: M = E - e * sin(E)ul6
An auxiliary angle to compute the position in an elliptic orbite
©緣生術數研究社 -- 術數研究 Yo2.J
v True anomaly: the angle from perihelion to the planet, as seen"W~$
from the Sun9P
©緣生術數研究社 -- 術數研究 +V|9q
r Heliocentric distance: the planet's distance from the Sun. l*V
©緣生術數研究社 -- 術數研究 Fs{
This relation is valid for an elliptic orbit:A(
©緣生術數研究社 -- 術數研究 8FM.K;
r * cos(v) = a * (cos(E) - e)oy;
r * sin(v) = a * sqrt(1 - e*e) * sin(E)~
©緣生術數研究社 -- 術數研究 '%
x,y,z Rectangular coordinates. Used e.g. when a heliocentricQvra)o
position (seen from the Sun) should be converted to au
corresponding geocentric position (seen from the Earth).6;
©緣生術數研究社 -- 術數研究 4
2. Some useful functionsT
When doing these orbital computations it's useful to have access to several utility functions. Some of them are in the standard library of the programming language, but others must be added by the programmer. Pocket calculators are often better equipped: they usually have sin/cos/tan and their inverses in both radians and degrees. Often one also finds functions to directly convert between rectangular and polar coordinates.0(<`2}
©緣生術數研究社 -- 術數研究 _
On microcomputers the situation is worse. Let's start with the programming language Basic: there we can count on having sin/cos/tan and atn (=arctan) in radians, but nothing more. arcsin and arccos are missing, and the trig functions don't work in degrees. However one can add one's own function library, like e.g. this:,a>
©緣生術數研究社 -- 術數研究 _:I>1H
95 rem Constants<~th9z
100 pi = 3.14159265359Skj4b
110 radeg = 180/pihnLZm
©緣生術數研究社 -- 術數研究 z=IE7H
115 rem arcsin and arccosP
120 def fnasin(x) = atn(x/sqr(1-x*x))@7
130 def fnacos(x) = pi/2 - fnasin(x)7lC_ri
©緣生術數研究社 -- 術數研究 9se
135 rem Trig. functions in degrees@~M
140 def fnsind(x) = sin(x/radeg) P:
150 def fncosd(x) = cos(x/radeg)BY
160 def fntand(x) = tan(x/radeg)A
170 def fnasind(x) = radeg*atn(x/sqr(1-x*x))*to
180 def fnacosd(x) = 180 - fnasind(x).qT
190 def fnatand(x) = radeg*atn(x)|/
©緣生術數研究社 -- 術數研究 E*)e;
195 rem arctan in all four quadrants_uH$@V
200 def fnatan2(y,x) = atn(y/x) - pi*(x<0)y<Yn:R
210 def fnatan2d(y,x) = radeg*atn(y/x) - 180*(x<0)*o0;/
©緣生術數研究社 -- 術數研究 -T
215 rem Normalize an angle between 0 and 360 degrees}3|I9h
217 rem Use Double Precision, if possible>G"~]
220 def fnrev#(x#) = x# - int(x#/360#)*360#S%izU
©緣生術數研究社 -- 術數研究 cfq#
225 rem Cube Root (needed for parabolic orbits)Ie:X$M
230 def fncbrt(x) = exp( log(x)/3 ){2Rh
©緣生術數研究社 -- 術數研究 .lza
The code above follows the convensions of traditional Microsoft Basic (MBASIC/BASICA/GWBASIC). If you use some other Basic dialect, you may want to modify this code.Y/rL
©緣生術數研究社 -- 術數研究 (<gX
The code above gives you sin/cos/tan and their inverses, in radians and degrees. The functions fnatan2() and fnatan2d() may need some explanation: they take two arguments, x and y, and they compute arctan(y/x) but puts the angle in the correct quadrant from -180 to +180 degrees. This is really part of a conversion from rectangular to polar coordinates, where the angle is computed. The distance is then computed by:[JB{Fw
©緣生術數研究社 -- 術數研究 ` E
sqrt( x*x + y*y )NqR`v
©緣生術數研究社 -- 術數研究 -I#
The function fnrev# normalizes an angle to between 0 and 360 degrees, by adding or subtracting even mutiples of 360 degrees until the final value falls between 0 and 360. The # sign means that the function and its argument use double precision. It is essential that this function uses more than single precision. More about this later.!\6
©緣生術數研究社 -- 術數研究 m\b
One warning: some of these functions may divide by zero. If one tries to compute fnasin(1.0), it's computed as: atn(1/sqrt(1-1.0*1.0)) = atn(1/sqr(0)) = atn(1/0). This is not such a big problem, since in practice one rarely tries to compute the arc sine of exactly 1.0. Also, some dialects of Microsoft Basic then just print the warning message "Overflow", compute 1/0 as the highest possible floating-point number, and then continue the program. When computing the arctan of this very large number, one will get pi/2 (or 90 degrees), which is the correct result. However, if you have access to a more modern, structured, Basic, with the ability to define multi-line function, then by all means use this to write a better version of arcsin, which treats arcsin(1.0) as a special case.2
©緣生術數研究社 -- 術數研究 Xa/8a]
There's a similar problem with the fncbrt() function - it only works for positive values. With a multi-line function definition it can be rewritten to work for negative values and zero as well, if one follows these simple rules: the Cube Root of zero is of course zero. The Cube Root of a negative number is computed by making the number positive, taking the Cube Root of that positive number, and then negating the result.g
©緣生術數研究社 -- 術數研究 OP]Bh?
Two other popular programming language are C and C++. The standard library of these languages are better equipped with trigonometric functions. You'll find sin/cos/tan and their inverses, and even an atan2() function among the standard functions. All you need to do is to define some macros to get the trig functions in degrees (include all parentheses in these macro definitions), and to define a rev() function which reduces an angle to between 0 and 360 degrees, and a cbrt() function which computes the cube root.MP
©緣生術數研究社 -- 術數研究 e;;;je
#define PI 3.141592653589793238461? Q
#define RADEG (180.0/PI)J
#define DEGRAD (PI/180.0)h#io
#define sind(x) sin((x)*DEGRAD)l**YjK
#define cosd(x) cos((x)*DEGRAD)Vd`(k6
#define tand(x) tan((x)*DEGRAD)!Nr
#define asind(x) (RADEG*asin(x))qCx~
#define acosd(x) (RADEG*acos(x))cPSRd
#define atand(x) (RADEG*atan(x))f2
#define atan2d(y,x) (RADEG*atan2((y),(x)))TR;
©緣生術數研究社 -- 術數研究 4
double rev( double x )<
{jzgf
return x - floor(x/360.0)*360.0;!U.!
}j%lYA
©緣生術數研究社 -- 術數研究 0^51-
double cbrt( double x )l_Z6?!
{iP;1
if ( x > 0.0 )s
return exp( log(x) / 3.0 );i}qKu
else if ( x < 0.0 )e3?^
return -cbrt(-x);J
else /* x == 0.0 */bmB`.[
return 0.0;Y&&f
}'y
©緣生術數研究社 -- 術數研究 TOy
In C++ the macros could preferably be defined as inline functions instead - this enables better type checking and also makes overloading of these function names possible.J?s^
©緣生術數研究社 -- 術數研究 x
The good ol' programming langauge FORTRAN is also well equipped with standard library trig functions: we find sin/cos/tan + inverses, and also an atan2 but only for radians. Therefore we need several function definitions to get the trig functions in degrees too. Below I give code only for SIND, ATAND, ATAN2D, plus REV and CBRT. The remaining functions COSD, TAND, ASIND and ACOSD are written in a similar way:.@
©緣生術數研究社 -- 術數研究 qyicD8
FUNCTION SIND(X)MK1c?
PARAMETER(RADEG=57.2957795130823)PD
SIND = SIN( X * (1.0/RADEG) )\N?
END4J
©緣生術數研究社 -- 術數研究 A^d2
FUNCTION ATAND(X)F
PARAMETER(RADEG=57.2957795130823)D?{7
ATAND = RADEG * ATAN(X)Mb
ENDHwQF)b
©緣生術數研究社 -- 術數研究 z>
FUNCTION ATAN2D(Y,X)W
PARAMETER(RADEG=57.2957795130823)F
ATAN2D = RADEG * ATAN2(Y,X):
END[U-=e
©緣生術數研究社 -- 術數研究 ynH
FUNCTION REV(X)Y+xKwT
REV = X - AINT(X/360.0)*360.04BQ.
IF (REV.LT.0.0) REV = REV + 360.0+orU
ENDo0CT{
©緣生術數研究社 -- 術數研究 xW`9^
FUNCTION CBRT(X)&A"=
IF (X.GE.0.0) THENARa
CBRT = X ** (1.0/3.0)zF=8
ELSE7+/fAu
CBRT = -((-X)**(1.0/3.0))EkRGPf
ENDIFdr!T9h
END##p,
©緣生術數研究社 -- 術數研究 f
The programming language Pascal is not as well equipped with trig functions. We have sin, cos, tan and arctan but nothing more. Therefore we need to write our own arcsin, arccos and arctan2, plus all the trig functions in degrees, and also the functions rev and cbrt. The trig functions in degrees are trivial when the others are defined, therefore I only define arcsin, arccos, arctan2, rev and cbrt:.89AyJ
©緣生術數研究社 -- 術數研究 EC8<K
const pi = 3.14159265358979323846;-<GZ8
half_pi = 1.57079632679489661923;_<yR
©緣生術數研究社 -- 術數研究 @jV5$q
function arcsin( x : real ) : real;W([
begin@uc|
if x = 1.0 theni3P&h8
arcsin := half_pi54z>UG
else if x = -1.0 thenk2)^O
arcsin := -half_piXpl{
else{ x
arcsin := arctan( x / sqrt( 1.0 - x*x ) )Jk=
end;Vx
©緣生術數研究社 -- 術數研究 GZ<)M
function arccos( x : real ) : real;)
beginU^R0
arccos := half_pi - arcsin(x);+{"
end;O oV
©緣生術數研究社 -- 術數研究 ./hSu
function arctan2( y, x : real ) : real;K&
beginlNoH.
if x = 0.0 then#0
begin&3_Y
if y = 0.0 thenfGk$d
(* Error! Give error message and stop program *)%w\7,8
else if y > 0.0 then+Q
arctan2 := half_pi]}
elsehk'D
arctan2 := -half_piH
endy-Ym
elsenrf&St
begin1
if x > 0.0 then"v7
arctan2 := arctan( y / x )G4
else if x < 0.0 then|AyKz
beginw.
if y >= 0.0 theng>gC!
arctan2 := arctan( y / x ) + pi'N]
elseRY`{q
arctan2 := arctan( y / x ) - pic-
end;7z
end;z3gI
end;M4|_)*
©緣生術數研究社 -- 術數研究 Z#IX,
function rev( x : real ) : real;&x#4p
var rv : real;TA
begind`'7~
rv := x - trunc(x/360.0)*360.0;f
if rv < 0.0 thenxMk2
rv := rv + 360.0;3JcE._
rev := rv;F1]E6G
end;VZ7R'9
©緣生術數研究社 -- 術數研究 S/r_[
function cbrt( x : real ) : real;^"{
begin}@<H?
if x > 0.0 then@^mx+T
cbrt := exp ( ln(x) / 3.0 );,
else if x < 0.0 then*NU8X
cbrt := -cbrt(-x)hb\
elseL
cbrt := 0.0W{lj
end;yYVq\
©緣生術數研究社 -- 術數研究 h"
It's well worth the effort to ensure that all these functions are available. Then you don't need to worry about these details which really don't have much to do with the problem of computing a planetary position.:~=+$
©緣生術數研究社 -- 術數研究 ICVJfR
3. Rectangular and spherical coordinatesJY<T'r
The position of a planet can be given in one of several ways. Two different ways that we'll use are rectangular and spherical coordinates.4!
©緣生術數研究社 -- 術數研究 =b
Suppose a planet is situated at some RA, Decl and r, where RA is the Right Ascension, Decl the declination, and r the distance in some length unit. If r is unknown or irrelevant, set r = 1. Let's convert this to rectangular coordinates, x,y,z:P'BG
©緣生術數研究社 -- 術數研究 PxF/
x = r * cos(RA) * cos(Decl)}q
y = r * sin(RA) * cos(Decl)|%L8
z = r * sin(Decl)V
©緣生術數研究社 -- 術數研究 y
(before we compute the sine/cosine of RA, we must first convert RA from hours/minutes/seconds to hours + decimals. Then the hours are converted to degrees by multiplying by 15);H&X
©緣生術數研究社 -- 術數研究 ;
If we know the rectangular coordinates, we can convert to spherical coordinates by the formulae below:|%_B<I
©緣生術數研究社 -- 術數研究 Z:a!?
r = sqrt( x*x + y*y + z*z )@<YY%|
RA = atan2( y, x )^T
Decl = asin( z / r ) = atan2( z, sqrt( x*x + y*y ) )-wQ:
©緣生術數研究社 -- 術數研究 mQ vY
At the north and south celestial poles, both x and y are zero. Since atan2(0,0) is undefined, the RA is undefined too at the celestial poles. The simplest way to handle this is to assign RA some arbitrary value, e.g. zero. Close to the celestial poles the formula asin(z/r) to compute the declination becomes sensitive to round-off errors - here the formula atan2(z,sqrt(x*x+y*y)) is preferable.(FJ#
©緣生術數研究社 -- 術數研究 JN&]
Not only equatorial coordinates can be converted between spherical and rectangular. These conversions can also be applied to ecliptic and horizontal coordinates. Just exchange RA,Decl with long,lat (ecliptic coordinates) or azimuth,altitude (horizontal coordinates).c^wF
©緣生術數研究社 -- 術數研究 /qsm
A coordinate system can be rotated. If a rectangular coordinate system is rotated around, say, the X axis, one can easily compute the new x,y,z coordinates. As an example, let's consider rotating an ecliptic x,y,z system to an equatorial x,y,z system. This rotation is done around the X axis (which points to the Vernal Point, the common point of origin in ecliptic and equatorial coordinates), through an angle of oblecl (the obliquity of the ecliptic, which is approximately 23.4 degrees):JSi
©緣生術數研究社 -- 術數研究 J
xequat = xeclipp
yequat = yeclip * cos(oblecl) - zeclip * sin(oblecl)We^
zequat = yeclip * sin(oblecl) + zeclip * cos(oblecl)eBvv:
©緣生術數研究社 -- 術數研究 x5Gn(
Now the x,y,z system is equatorial. It's easily rotated back to ecliptic coordinates by simply switching sign on oblecl:@~mW
©緣生術數研究社 -- 術數研究 o
xeclip = xequatxfO R
yeclip = yequat * cos(-oblecl) - zequat * sin(-oblecl)9|JgZ,
zeclip = yequat * sin(-oblecl) + zequat * cos(-oblecl)1P*P
©緣生術數研究社 -- 術數研究 eTg1tB
When computing sin and cos of -oblecl, one can use the identities:|br
©緣生術數研究社 -- 術數研究 i
cos(-x) = cos(x), sin(-x) = -sin(x)"
©緣生術數研究社 -- 術數研究
Now let's put this together to convert directly from spherical ecliptic coordinates (long, lat) to spherical equatorial coordinates (RA, Decl). Since the distance r is irrelevant in this case, let's set r=1 for simplicity.ZEF
©緣生術數研究社 -- 術數研究 >^*i r
Example: At the Summer Solstice the Sun's ecliptic longitude is 90 degrees. The Sun's ecliptic latitude is always very nearly zero. Suppose the obliquity of the ecliptic is 23.4 degrees:~W\PAh
©緣生術數研究社 -- 術數研究 X
xeclip = cos(90_deg) * cos(0_deg) = 0.0000Z;AVr
yeclip = sin(90_deg) * cos(0_deg) = 1.0000ZzWn"
zeclip = sin(0_deg) = 0.0000i>Pfo
©緣生術數研究社 -- 術數研究 E
Rotate through oblecl = 23.4_deg:k
©緣生術數研究社 -- 術數研究 J+d
xequat = 0.0000~
yequat = 1.0000 * cos(23.4_deg) - 0.0000 * sin(23.4_deg)8
zequat = 1.0000 * sin(23.4_deg) + 0.0000 * cos(23.4_deg)}
©緣生術數研究社 -- 術數研究 4dN+x
Our equatorial rectangular coordinates become:X\
©緣生術數研究社 -- 術數研究 R^*
x = 0jT/
y = cos(23.4_deg) = 0.9178q#eq
z = sin(23.4_deg) = 0.3971<
©緣生術數研究社 -- 術數研究 i
The "distance", r, becomes: sqrt( 0.8423 + 0.1577 ) = 1.0000 i.e. unchanged&~Hef!
©緣生術數研究社 -- 術數研究 Oe|C
RA = atan2( 0.9178, 0 ) = 90_deg)|9
Decl = asin( 0.3971 / 1.0000 ) = 23.40_deg0=_
©緣生術數研究社 -- 術數研究 x3"7KG
Alternatively:kQ*U#
©緣生術數研究社 -- 術數研究 *.i
Decl = atan2( 0.3971, sqrt( 0.8423 + 0.0000 ) ) = 23.40_deg==N]
©緣生術數研究社 -- 術數研究 D
Here we immediately see how simple it is to compute RA, thanks to the atan2() function: no need to consider in which quadrant it falls, the atan2() function handles this.349Ky8
©緣生術數研究社 -- 術數研究 Q3
4. The time scale. A test date._/
The time scale used here is a "day number" from 2000 Jan 0.0 TDT, which is the same as 1999 Dec 31.0 TDT, i.e. precisely at midnight TDT at the start of the last day of this century. With the modest accuracy we strive for here, one can usually disregard the difference between TDT (formerly canned ET) and UT.zA{~
©緣生術數研究社 -- 術數研究 HYw9H$
We call our day number d. It can be computed from a JD (Julian Day Number) or MJD (Modified Julian Day Number) like this:1LP
©緣生術數研究社 -- 術數研究 ?_D
d = JD - 2451543.5 = MJD - 51543.0)@)
©緣生術數研究社 -- 術數研究 {.7e
We can also compute d directly from the calendar date like this:dgCR#>
©緣生術數研究社 -- 術數研究 ur[]i
d = 367*Y - (7*(Y + ((M+9)/12)))/4 + (275*M)/9 + D - 730530J
©緣生術數研究社 -- 術數研究 ;9\zn
Y is the year (all 4 digits!), M the month (1-12) and D the date. In this formula all divisions should be INTEGER divisions. Use "div" instead of "/" in Pascal, and "\" instead of "/" in Microsoft Basic. In C/C++ and FORTRAN it's sufficient to ensure that both operands to "/" are integers.ne~
©緣生術數研究社 -- 術數研究 4,
This formula yields d as an integer, which is valid at the start (at midnight), in UT or TDT, of that calendar date. If you want d for some other time, add UT/24.0 (here the division is a floating-point division!) to the d obtained above.u]O|
©緣生術數研究社 -- 術數研究 4
Example: compute d for 19 april 1990, at 0:00 UT :oYEC0K
©緣生術數研究社 -- 術數研究 Y5]
We can look up, or compute the JD for this moment, and we'll get: JD = 2448000.5 which yields d = -3543.02X
©緣生術數研究社 -- 術數研究 A?H*Ff
Or we can compute d directly from the calendar date:Z4D]iF
©緣生術數研究社 -- 術數研究 E!o6B
d = 367*1990 - (7*(1990 + ((4+9)/12)))/4 + (275*4)/9 + 19 - 730530yq
©緣生術數研究社 -- 術數研究 H#F=w
d = 730330 - (7*(1990 + (13/12)))/4 + 1100/9 + 19 - 730530>7+')
©緣生術數研究社 -- 術數研究 ("O<}{
d = 730330 - (7*(1990 + 1))/4 + 122 + 19 - 730530%;fA
©緣生術數研究社 -- 術數研究 *dX@[3
d = 730330 - (7*1991)/4 + 122 + 19 - 730530=O]Yk!
©緣生術數研究社 -- 術數研究 ,Uscx|
d = 730330 - 13937/4 + 122 + 19 - 730530~D_O
©緣生術數研究社 -- 術數研究 >k
d = 730330 - 3484 + 122 + 19 - 730530 = -3543az0
©緣生術數研究社 -- 術數研究 Zrt
This moment, 1990 april 19, 0:00 UT/TDT, will be our test date for most numerical examples below. d is negative since our test date, 19 april 1990, is earlier than the "point of origin" of our day number, 31 dec 1999.`
©緣生術數研究社 -- 術數研究 ~G\~
5. The Sun's position.-y
Today most people know that the Earth orbits the Sun and not the other way around. But below we'll pretend as if it was the other way around. These orbital elements are thus valid for the Sun's (apparent) orbit around the Earth. All angular values are expressed in degrees:ZK@~
©緣生術數研究社 -- 術數研究 5I<~j
w = 282.9404_deg + 4.70935E-5_deg * d (longitude of perihelion)A5Kem$
a = 1.000000 (mean distance, a.u.)A-lS
e = 0.016709 - 1.151E-9 * d (eccentricity)BUhl(
M = 356.0470_deg + 0.9856002585_deg * d (mean anomaly)Lj_
©緣生術數研究社 -- 術數研究 f*l8%z
We also need the obliquity of the ecliptic, oblecl:{NQ
©緣生術數研究社 -- 術數研究 o
oblecl = 23.4393_deg - 3.563E-7_deg * d3TGXQ
©緣生術數研究社 -- 術數研究 T
and the Sun's mean longitude, L:2Z=I
©緣生術數研究社 -- 術數研究 dr/vho
L = w + M<-7~u
©緣生術數研究社 -- 術數研究 !v`*w
By definition the Sun is (apparently) moving in the plane of the ecliptic. The inclination, i, is therefore zero, and the longitude of the ascending node, N, becomes undefined. For simplicity we'll assign the value zero to N, which means that w, the angle between acending node and perihelion, becomes equal to the longitude of the perihelion.)Nu%
©緣生術數研究社 -- 術數研究 8x []
Now let's compute the Sun's position for our test date 19 april 1990. Earlier we've computed d = -3543.0 which yields:@z
©緣生術數研究社 -- 術數研究 46dHn
w = 282.7735_deg?K ~
a = 1.000000%2ysQz
e = 0.016713\uD
M = -3135.9347_deg"!v}
©緣生術數研究社 -- 術數研究 R-s3&
We immediately notice that the mean anomaly, M, will get a large negative value. We use our function rev() to reduce this value to between 0 and 360 degrees. To do this, rev() will need to add 9*360 = 3240 degrees to this angle:\"?
©緣生術數研究社 -- 術數研究 #76U2
M = 104.0653_degD4,
©緣生術數研究社 -- 術數研究 2sH5@=
We also compute:uQ
©緣生術數研究社 -- 術數研究 t
L = w + M = 386.8388_deg = 26.8388_deg9AWv
©緣生術數研究社 -- 術數研究 i5
oblecl = 23.4406_deg}
©緣生術數研究社 -- 術數研究 \~zH2
Let's go on computing an auxiliary angle, the eccentric anomaly. Since the eccentricity of the Sun's (i.e. the Earth's) orbit is so small, 0.017, a first approximation of E will be accurate enough. Below E and M are in degrees:i?(
©緣生術數研究社 -- 術數研究 n(
E = M + (180/pi) * e * sin(M) * (1 + e * cos(M))6g
©緣生術數研究社 -- 術數研究 #W
When we plug in M and e, we get:WhD
©緣生術數研究社 -- 術數研究 :O^y
E = 104.9904_degzW
©緣生術數研究社 -- 術數研究 +\/
Now we compute the Sun's rectangular coordinates in the plane of the ecliptic, where the X axis points towards the perihelion:iWR|B{
©緣生術數研究社 -- 術數研究 #<8%
x = r * cos(v) = cos(E) - e`Fl6
y = r * sin(v) = sin(E) * sqrt(1 - e*e)pU
©緣生術數研究社 -- 術數研究 pO+y
We plug in E and get:^yO1S
©緣生術數研究社 -- 術數研究 XOZ+R=
x = -0.275370&/S'=D
y = +0.965834p!Em
©緣生術數研究社 -- 術數研究 +3
Convert to distance and true anomaly:f9
©緣生術數研究社 -- 術數研究 2A[pG3
r = sqrt(x*x + y*y)9
v = arctan2( y, x )8FN"
©緣生術數研究社 -- 術數研究 efNzM
Numerically we get:+/\4N
©緣生術數研究社 -- 術數研究 z
r = 1.004323H
v = 105.9134_deg$
©緣生術數研究社 -- 術數研究 l.6
Now we can compute the longitude of the Sun:y?P
©緣生術數研究社 -- 術數研究 O
lon = v + w:'^Mq]
©緣生術數研究社 -- 術數研究 [[6c/
lon = 105.9134_deg + 282.7735_deg = 388.6869_deg = 28.6869_degjs
©緣生術數研究社 -- 術數研究 Mw@
We're done!k
©緣生術數研究社 -- 術數研究 -F
How close did we get to the correct values? Let's compare with the Astronomical Almanac:h'Tjyn
©緣生術數研究社 -- 術數研究 mEfjy5
Our results Astron. Almanac Differenceqh3&Z
©緣生術數研究社 -- 術數研究 Ek
lon 28.6869_deg 28.6813_deg +0.0056_deg = 20"ITk
r 1.004323 1.004311 +0.000012FSU?d@
©緣生術數研究社 -- 術數研究 (F[T
The error in the Sun's longitude was 20 arc seconds, which is well below our aim of an accuracy of one arc minute. The error in the distance was about 1/3 Earth radius. Not bad!TPYi
©緣生術數研究社 -- 術數研究 7`H
Finally we'll compute the Sun's ecliptic rectangular coordinates, rotate these to equatorial coordinates, and then compute the Sun's RA and Decl:G.
©緣生術數研究社 -- 術數研究 oL
x = r * cos(lon)l
y = r * sin(lon)pV
z = 0.0+)b1h
©緣生術數研究社 -- 術數研究 gC
We plug in our longitude:[y3
©緣生術數研究社 -- 術數研究 N.
x = 0.881048_
y = 0.482098cM^8
z = 0.0QTnJIk
©緣生術數研究社 -- 術數研究 'Ym[r~
We use oblecl = 23.4406 degrees, and rotate these coordinates:#Q
©緣生術數研究社 -- 術數研究 W~?8*'
xequat = 0.881048j{Z-
yequat = 0.482098 * cos(23.4406_deg) + 0.0 * sin(23.4406_deg)5a3|AO
zequat = 0.482098 * sin(23.4406_deg) + 0.0 * cos(23.4406_deg)[
©緣生術數研究社 -- 術數研究 y
which yields:"zF
©緣生術數研究社 -- 術數研究 Y`v.q
xequat = 0.88104881Ho*
yequat = 0.442312;Kt~M
zequat = 0.191778Ou
©緣生術數研究社 -- 術數研究 RaKmN
Convert to RA and Decl:XI
©緣生術數研究社 -- 術數研究 8L
r = 1.004323 (unchanged)=
RA = 26.6580_deg = 26.6580/15 h = 1.77720 h = 1h 46m 37.9sKv
Decl = +11.0084_deg = +11_deg 0' 30"{qSbc
©緣生術數研究社 -- 術數研究 <>
The Astronomical Almanac says:.qz`#_
©緣生術數研究社 -- 術數研究 /40H(
RA = 1h 46m 36.0s Decl = +11_deg 0' 22">
©緣生術數研究社 -- 術數研究 ecZ
6. Sidereal time and hour angle. Altitude and azimuth.?_
The Sidereal Time tells the Right Ascension of the part of the sky that's precisely south, i.e. in the meridian. Sidereal Time is a local time, which can be computed from:M@I^^
©緣生術數研究社 -- 術數研究 B"df)|
SIDTIME = GMST0 + UT + LON/15Op'w
©緣生術數研究社 -- 術數研究 u[
where SIDTIME, GMST0 and UT are given in hours + decimals. GMST0 is the Sidereal Time at the Greenwich meridian at 00:00 right now, and UT is the same as Greenwich time. LON is the terrestial longitude in degrees (western longitude is negative, eastern positive). To "convert" the longitude from degrees to hours we divide it by 15. If the Sidereal Time becomes negative, we add 24 hours, if it exceeds 24 hours we subtract 24 hours.o"c#
©緣生術數研究社 -- 術數研究 Yvs7
Now, how do we compute GMST0? Simple - we add (or subtract) 180 degrees to (from) L, the Sun's mean longitude, which we've already computed earlier. Then we normalise the result to between 0 and 360 degrees, by applying the rev() function. Finally we divide by 15 to convert degrees to hours:f-8j>*
©緣生術數研究社 -- 術數研究 lbT
GMST0 = ( L + 180_deg ) / 15 = L/15 + 12ha
©緣生術數研究社 -- 術數研究 AO
We've already computed L = 26.8388_deg, which yields:wGq
©緣生術數研究社 -- 術數研究 *^
GMST0 = 26.8388_deg/15 + 12h = 13.78925 hoursk7
©緣生術數研究社 -- 術數研究 e`{
Now let's compute the local Sidereal Time for the time meridian of Central Europe (at 15 deg east longitude = +15 degrees long) on 19 april 1990 at 00:00 UT:'3
©緣生術數研究社 -- 術數研究 RFl_\2
SIDTIME = GMST0 + UT + LON/15 = 13.78925h + 0 + 15_deg/15 = 14.78925 hours[QvN
©緣生術數研究社 -- 術數研究 ;;
SIDTIME = 14h 47m 21.3se\C>l
©緣生術數研究社 -- 術數研究 fn
To compute the altitude and azimuth we also need to know the Hour Angle, HA. The Hour Angle is zero when the clestial body is in the meridian i.e. in the south (or, from the southern heimsphere, in the north) - this is the moment when the celestial body is at its highest above the horizon.5l}fY8
©緣生術數研究社 -- 術數研究 Q>I
The Hour Angle increases with time (unless the object is moving faster than the Earth rotates; this is the case for most artificial satellites). It is computed from:nW
©緣生術數研究社 -- 術數研究 mXpXr
HA = SIDTIME - RAaR+
©緣生術數研究社 -- 術數研究 l
Here SIDTIME and RA must be expressed in the same unit, hours or degrees. We choose hours:DM4.
©緣生術數研究社 -- 術數研究 f[I&Y
HA = 14.78925h - 1.77720h = 13.01205h = 195.1808_deg:
©緣生術數研究社 -- 術數研究 vi=fe
If the Hour Angle is 180_deg the celestial body can be seen (or not be seen, if it's below the horizon) in the north (or in the south, from the southern hemisphere). We get HA = 195_deg for the Sun, which seems OK since it's around 01:00 local time.VJL$%i
©緣生術數研究社 -- 術數研究 t|0
Now we'll convert the Sun's HA = 195.1808_deg and Decl = +11.0084_deg to a rectangular (x,y,z) coordinate system where the X axis points to the celestial equator in the south, the Y axis to the horizon in the west, and the Z axis to the north celestial pole: The distance, r, is here irrelevant so we set r=1 for simplicity:4
©緣生術數研究社 -- 術數研究 _yi
x = cos(HA) * cos(Decl) = -0.947346Xb
y = sin(HA) * cos(Decl) = -0.257047^
z = sin(Decl) = +0.190953f[[quS
©緣生術數研究社 -- 術數研究 V&
Now we'll rotate this x,y,z system along an axis going east-west, i.e. the Y axis, in such a way that the Z axis will point to the zenith. At the North Pole the angle of rotation will be zero since there the north celestial pole already is in the zenith. At other latitudes the angle of rotation becomes 90_deg - latitude. This yields:dLg$a
©緣生術數研究社 -- 術數研究 8o
xhor = x * cos(90_deg - lat) - z * sin(90_deg - lat)K+QL-
yhor = y)
zhor = x * sin(90_deg - lat) + z * cos(90_deg - lat)~&N
©緣生術數研究社 -- 術數研究 0\vl
Since sin(90_deg-lat) = cos(lat) (and reverse) we'll get:])3F=J
©緣生術數研究社 -- 術數研究 K^
xhor = x * sin(lat) - z * cos(lat)`KJ6
yhor = y"E
zhor = x * cos(lat) + z * sin(lat))C`zE2
©緣生術數研究社 -- 術數研究 fWUr2
Finally we compute our azimuth and altitude:#iNnP
©緣生術數研究社 -- 術數研究 h1
azimuth = atan2( yhor, xhor ) + 180_deg4F'%9u
altitude = asin( zhor ) = atan2( zhor, sqrt(xhor*xhor+yhor*yhor) )c<7
©緣生術數研究社 -- 術數研究 L!c&,
Why did we add 180_deg to the azimuth? To adapt to the most common way to specify azimuth: from North (0_deg) through East (90_deg), South (180_deg), West (270_deg) and back to North. If we didn't add 180_deg the azimuth would be counted from South through West/etc instead. If you want to use that kind of azimuth, then don't add 180_deg above.;4Nh`?
©緣生術數研究社 -- 術數研究 o[t\+b
We select some place in central Scandinavia: the longitude is as before +15_deg (15_deg East), and the latitude is +60_deg (60_deg N):c}yz$%
©緣生術數研究社 -- 術數研究 -0Wz
xhor = -0.947346 * sin(60_deg) - (+0.190953) * cos(60_deg) = -0.915902L
yhor = -0.257047 = -0.257047g&
zhor = -0.947346 * cos(60_deg) + (+0.190953) * sin(60_deg) = -0.308303u%nkCe
©緣生術數研究社 -- 術數研究 [hy
Now we've computed the horizontal coordinates in rectangular form. To get azimuth and altitude we convert to spherical coordinates (r=1):_-
©緣生術數研究社 -- 術數研究 .ud\?!
azimuth = atan2(-0.257047,-0.915902) + 180_deg = 375.6767_deg = 15.6767_degD
altitude = asin( -0.308303 ) = -17.9570_deg%s
©緣生術數研究社 -- 術數研究 UgbZ9V
Let's round the final result to two decimals:3z
©緣生術數研究社 -- 術數研究 CJ_^Z
azimuth = 15.68_deg, altitude = -17.96_deg.AJjF
©緣生術數研究社 -- 術數研究 k
The Sun is thus 17.96_deg below the horizon at this moment and place. This is very close to astronomical twilight (18_deg below the horizon).`7uG
©緣生術數研究社 -- 術數研究 nv|,]{
7. The Moon's positionDg,aG
Let's continue by computing the position of the Moon. The computatons will become more complicated, since the Moon doesn't move in the plane of the ecliptic, but in a plane inclined somewhat more than 5 degrees to the ecliptic. Also, the Sun perturbs the Moon's motion significantly, an effect we must account for.4
©緣生術數研究社 -- 術數研究 3;
The orbital elements of the Moon are:<9|
©緣生術數研究社 -- 術數研究 H%
N = 125.1228_deg - 0.0529538083_deg * d (Long asc. node)IY
i = 5.1454_deg (Inclination)U+]u
w = 318.0634_deg + 0.1643573223_deg * d (Arg. of perigee)5&Mx
a = 60.2666 (Mean distance)Ov[8$
e = 0.054900 (Eccentricity)qA4
M = 115.3654_deg + 13.0649929509_deg * d (Mean anomaly)=Ycq
©緣生術數研究社 -- 術數研究 s`*
The Moon's ascending node is moving in a retrogade ("backwards") direction one revolution in about 18.6 years. The Moon's perigee (the point of the orbit closest to the Earth) moves in a "forwards" direction one revolution in about 8.8 years. The Moon itself moves one revolution in aboout 27.5 days. The mean distance, or semi-major axis, is expressed in Earth equatorial radii).1$@t
©緣生術數研究社 -- 術數研究 F#q-
Let's compute numerical values for the Moon's orbital elements on our test date 19 april 1990 (d = -3543):'rL588
©緣生術數研究社 -- 術數研究 :
N = 312.7381_degZUtnc
i = 5.1454_degiD#d
w = -264.2546_degH
a = 60.2666 (Earth equatorial radii)reV=x_
e = 0.054900a
M = -46173.9046_dege;A
©緣生術數研究社 -- 術數研究 ]t>Lo5
Now the need for sufficient numerical accuracy becomes obvious. If we would compute M with normal single precision, i.e. only 7 decimal digits of accuracy, then the error in M would here be about 0.01 degrees. Had we selected a date around 1901 or 2099 then the error in M would have been about 0.1 degrees, which is definitely worse than our aim of a maximum error of one or two arc minutes. Therefore, when computing the Moon's mean anomaly, M, it's important to use at least 9 or 10 digits of accuracy.x6
©緣生術數研究社 -- 術數研究 @bpsk
(This was a real problem around 1980, when microcomputers were a novelty. Around then, pocket calculators usually offered better precision than microcomputers. But these days are long gone: nowadays microcomputers routinely offer double precision (14-16 digits of accuracy) support in hardware; all you need to do is to select a compiler which really suports this.)/
©緣生術數研究社 -- 術數研究 K
All angular elements should be normalized to within 0-360 degrees, by calling the rev() function. We get:-41'(8
©緣生術數研究社 -- 術數研究 xNaW%
N = 312.7381_deg5
i = 5.1454_degaI
w = 95.7454_deg]s=
a = 60.2666 (Earth equatorial radii)Od+iR,
e = 0.054900[d
M = 266.0954_degs
©緣生術數研究社 -- 術數研究 b
To normalize M we had to add 129*360 = 46440 degrees.|,sO
©緣生術數研究社 -- 術數研究 @^p
Next, we compute E, the eccentric anomaly. We start with a first approximation (E0 and M in degrees):q
©緣生術數研究社 -- 術數研究 p3)&@
E0 = M + (180_deg/pi) * e * sin(M) * (1 + e * cos(M))wI=A
©緣生術數研究社 -- 術數研究 6Asd8B
The eccentricity of the Moon's orbit is larger than of the Earth's orbit. This means that our first approximation will have a bigger error - it'll be close to the limit of what we can tolerate within our accuracy aim. If you want to be careful, you should therefore use the iteration formula below: set E0 to our first approximation, compute E1, then set E0 to E1 and compute a new E1, until the difference between E0 and E1 becomes small enough, i.e. smaller than about 0.005 degrees. Then use the last E1 as the final value. In the formula below, E0, E1 and M are in degrees:/[p=
©緣生術數研究社 -- 術數研究 C91
E1 = E0 - (E0 - (180_deg/pi) * e * sin(E0) - M) / (1 - e * cos(E0))`B
©緣生術數研究社 -- 術數研究 (%
On our test date, the first approximation of E becomes: E=262.9689_deg The iterations then yield: E = 262.9735_deg, 262.9735_deg ......Q;E(
©緣生術數研究社 -- 術數研究 -
Now we've computed E - the next step is to compute the Moon's distance and true anomaly. First we compute rectangular (x,y) coordinates in the plane of the lunar orbit:HS
©緣生術數研究社 -- 術數研究 nGr5F)
x = r * cos(v) = a * (cos(E) - e)h
y = r * sin(v) = a * sqrt(1 - e*e) * sin(E)gFXtA
©緣生術數研究社 -- 術數研究 !
Our test date yields:d8:
©緣生術數研究社 -- 術數研究 g6
x = -10.68095"
y = -59.72377Gp<O"
©緣生術數研究社 -- 術數研究 (7IS}R
Then we convert this to distance and true anonaly:kMN^R
©緣生術數研究社 -- 術數研究 *ob]{'
r = sqrt( x*x + y*y ) = 60.67134 Earth radii>$[
v = atan2( y, x ) = 259.8605_deg8b/f
©緣生術數研究社 -- 術數研究 N
Now we know the Moon's position in the plane of the lunar orbit. To compute the Moon's position in ecliptic coordinates, we apply these formulae:p.E
©緣生術數研究社 -- 術數研究 r[`s
xeclip = r * ( cos(N) * cos(v+w) - sin(N) * sin(v+w) * cos(i) )+
yeclip = r * ( sin(N) * cos(v+w) + cos(N) * sin(v+w) * cos(i) )o
zeclip = r * sin(v+w) * sin(i)CXX|i
©緣生術數研究社 -- 術數研究 |0%n_
Our test date yields:\~a#.
©緣生術數研究社 -- 術數研究 ;H}
xeclip = +37.65311><:</
yeclip = -47.57180J$Yl)
zeclip = -0.41687W7B{2
©緣生術數研究社 -- 術數研究 H
Convert to ecliptic longitude, latitude, and distance:'
©緣生術數研究社 -- 術數研究 f
long = 308.3616_deg2CGev#
lat = -0.3937_deg1|V;
r = 60.6713V-r
©緣生術數研究社 -- 術數研究 M
According to the Astronomical Almanac, the Moon's position at this moment is 306.94_deg, and the latitude is -0.55_deg. This differs from our figures by 1.42_deg in longitude and 0.16_deg in latitude!!! This difference is much larger than our aim of an error of max 1-2 arc minutes. Why is this so?XKx
©緣生術數研究社 -- 術數研究 'CJ
8. The Moon's position with higher accuracy. Perturbations)jY6y
The big error in our computed lunar position is because we completely ignored the perturbations on the Moon. Below we'll compute the most important perturbation terms, and then add these as corrections to our previous figures. This will cut down the error to 1-2 arc minutes, or less.k1\
©緣生術數研究社 -- 術數研究 T
First we need several fundamental arguments:lp~d
©緣生術數研究社 -- 術數研究 )%biR
Sun's mean longitude: Ls (already computed)WW%e&*
Moon's mean longitude: Lm = N + w + M (for the Moon)]^pb1Q
Sun's mean anomaly: Ms (already computed)`sw7
Moon's mean anomaly: Mm (already computed)nLf
Moon's mean elongation: D = Lm - Lsw=#066
Moon's argument of latitude: F = Lm - N`ImS{
©緣生術數研究社 -- 術數研究 0?
Let's plug in the figures for our test date:F!
©緣生術數研究社 -- 術數研究 O;"kC=
Ms = 104.0653_degsNlC
Mm = 266.0954_degA?cFmZ
Ls = 26.8388_degJO]R
Lm = 312.7381_deg + 95.7454_deg + 266.0954_deg = 674.5789_deg
= 314.5789_degj{G
D = 314.5789_deg - 26.8388_deg = 287.7401_degu^Cn7.
F = 314.5789_deg - 312.7381_deg = 1.8408_degcX
©緣生術數研究社 -- 術數研究 \
Now it's time to compute and add up the 12 largest perturbation terms in longitude, the 5 largest in latitude, and the 2 largest in distance. These are all the perturbation terms with an amplitude larger than 0.01_deg in longitude resp latitude. In the lunar distance, only the perturbation terms larger than 0.1 Earth radii has been included:h|$
©緣生術數研究社 -- 術數研究 nCl *s
Perturbations in longitude (degrees):7
©緣生術數研究社 -- 術數研究 MJ0<V
-1.274_deg * sin(Mm - 2*D) (Evection)JV)Ns>
+0.658_deg * sin(2*D) (Variation)3;
-0.186_deg * sin(Ms) (Yearly equation)H3mGV{
-0.059_deg * sin(2*Mm - 2*D)q3m~r
-0.057_deg * sin(Mm - 2*D + Ms)6c.<
+0.053_deg * sin(Mm + 2*D)D*ju
+0.046_deg * sin(2*D - Ms)po!<[B
+0.041_deg * sin(Mm - Ms)T>
-0.035_deg * sin(D) (Parallactic equation)X
-0.031_deg * sin(Mm + Ms) t
-0.015_deg * sin(2*F - 2*D)g
+0.011_deg * sin(Mm - 4*D)vWjx@q
©緣生術數研究社 -- 術數研究 "*]uS
Perturbations in latitude (degrees):p|
©緣生術數研究社 -- 術數研究 `{#^M
-0.173_deg * sin(F - 2*D)KdR
-0.055_deg * sin(Mm - F - 2*D)n
-0.046_deg * sin(Mm + F - 2*D)yb*
+0.033_deg * sin(F + 2*D)Ht[
+0.017_deg * sin(2*Mm + F)PopI
©緣生術數研究社 -- 術數研究 OP0/!
Perturbations in lunar distance (Earth radii):5dN'
©緣生術數研究社 -- 術數研究 iE,gw=
-0.58 * cos(Mm - 2*D)<)P
-0.46 * cos(2*D)0q;b\:
©緣生術數研究社 -- 術數研究 0ftRIR
Some of the largest perturbation terms in longitude even have received individual names! The largest perturbation, the Evection, was discovered already by Ptolemy (he made it one of the epicycles in his theory for the Moon's motion). The two next largest perturbations, the Variation and the Yearly equation, were discovered by Tycho Brahe.E|p`IU
©緣生術數研究社 -- 術數研究 C92z*
If you don't need 1-2 arcmin accuracy, you don't need to compute all these perturbation terms. If you only include the two largest terms in longitude and the largest in latitude, the error in longitude rarely becomes larger than 0.25_deg, and in latitude rarely larger than 0.15_deg.M "
©緣生術數研究社 -- 術數研究 {9s F
Let's compute these perturbation terms for our test date:6
©緣生術數研究社 -- 術數研究 .T
longitude: -0.9847 - 0.3819 - 0.1804 + 0.0405 - 0.0244 + 0.0452 +'-RQ
0.0428 + 0.0126 - 0.0333 - 0.0055 - 0.0079 - 0.0029
= -1.4132_deg5tC+=G
©緣生術數研究社 -- 術數研究 rK;]^
latitude: -0.0958 - 0.0414 - 0.0365 - 0.0200 + 0.0018 = -0.1919_deg*`
©緣生術數研究社 -- 術數研究 KFO
distance: -0.3680 + 0.3745 = +0.0066 Earth radii9o
©緣生術數研究社 -- 術數研究 :e*
Add this to the ecliptic positions we earlier computed:7i <C3
©緣生術數研究社 -- 術數研究 {w5h5
long = 308.3616_deg - 1.4132_deg = 306.9484_degP+n;f
lat = -0.3937_deg - 0.1919_deg = -0.5856_degS[
dist = 60.6713 + 0.0066 = 60.6779 Earth radiiqs`|P{
©緣生術數研究社 -- 術數研究 6[~B~}
Let's compare with the Astronomical Almanac:]U.G,)
©緣生術數研究社 -- 術數研究 l
longitude 306.94_deg, latitude -0.55_deg, distance 60.793 Earth radiipA#!qs
©緣生術數研究社 -- 術數研究 \ujY:|
Now the agreement is much better, right?;'
©緣生術數研究社 -- 術數研究 t^KJ&o
Let's continue by converting these ecliptic coordinates to Right Ascension and Declination. We do as described earlier: convert the ecliptic longitude/latitude to rectangular (x,y,z) coordinates, rotate this x,y,z, system through an angle corresponding to the obliquity of the ecliptic, then convert back to spherical coordinates. The Moon's distance doesn't matter here, and one can therefore set r=1.0. We get:U6
©緣生術數研究社 -- 術數研究 U5,P3
RA = 309.5011_deg4v.3\a
Decl = -19.1032_degf`r
©緣生術數研究社 -- 術數研究 C
According to the Astronomical Almanac:5Nex
©緣生術數研究社 -- 術數研究 VXG$
RA = 309.4881_deg01
Decl = -19.0741_degOJ0
©緣生術數研究社 -- 術數研究 K+
9. The Moon's topocentric position.^lkSp9
The Moon's position, as computed earlier, is geocentric, i.e. as seen by an imaginary observer at the center of the Earth. Real observers dwell on the surface of the Earth, though, and they will see a different position - the topocentric position. This position can differ by more than one degree from the geocentric position. To compute the topocentric positions, we must add a correction to the geocentric position.6f6G^
©緣生術數研究社 -- 術數研究 f{
Let's start by computing the Moon's parallax, i.e. the apparent size of the (equatorial) radius of the Earth, as seen from the Moon:
©緣生術數研究社 -- 術數研究 n
mpar = asin( 1/r )%
©緣生術數研究社 -- 術數研究 kh"9(
where r is the Moon's distance in Earth radii. It's simplest to apply the correction in horizontal coordinates (azimuth and altitude): within our accuracy aim of 1-2 arc minutes, no correction need to be applied to the azimuth. One need only apply a correction to the altitude above the horizon:C
©緣生術數研究社 -- 術數研究 +
alt_topoc = alt_geoc - mpar * cos(alt_geoc)Ee:;'
©緣生術數研究社 -- 術數研究 G8l4K
Sometimes one needs to correct for topocentric position directly in equatorial coordinates though, e.g. if one wants to draw on a star map how the Moon passes in front of the Pleiades, as seen from some specific location. Then we need to know the Moon's geocentric Right Ascension and Declination (RA, Decl), the Local Sidereal Time (LST), and our latitude (lat).rr#;G
©緣生術數研究社 -- 術數研究 |v`kZ4
Our astronomical latitude (lat) must first be converted to a geocentric latitude (gclat) and distance from the center of the Earth (rho) in Earth equatorial radii. If we only want an approximate topocentric position, it's simplest to pretend that the Earth is a perfect sphere, and simply set:`n{Q
©緣生術數研究社 -- 術數研究 U
gclat = lat, rho = 1.0J)yU`0
©緣生術數研究社 -- 術數研究 ar_o
However, if we do wish to account for the flattening of the Earth, we instead compute:kr/bu6
©緣生術數研究社 -- 術數研究 ;h
gclat = lat - 0.1924_deg * sin(2*lat)x:(mP/
rho = 0.99883 + 0.00167 * cos(2*lat)"7_h
©緣生術數研究社 -- 術數研究 2.
Next we compute the Moon's geocentric Hour Angle (HA):-xW
©緣生術數研究社 -- 術數研究 ~EZ8
HA = LST - RA:ei*mx
©緣生術數研究社 -- 術數研究 .Gx
We also need an auxiliary angle, g:<
©緣生術數研究社 -- 術數研究 "
g = atan( tan(gclat) / cos(HA) )8h
©緣生術數研究社 -- 術數研究 aEZ1
Now we're ready to convert the geocentric Right Ascention and Declination (RA, Decl) to their topocentric values (topRA, topDecl):yc
©緣生術數研究社 -- 術數研究 ;T%xXS
topRA = RA - mpar * rho * cos(gclat) * sin(HA) / cos(Decl)QT
topDecl = Decl - mpar * rho * sin(gclat) * sin(g - Decl) / sin(g)O&
©緣生術數研究社 -- 術數研究 NZK
Let's do this correction for our test date and for the geographical position 15 deg E longitude (= +15_deg) and 60 deg N latitude (= +60_deg). Earlier we computed the Local Sidereal Time as LST = SIDTIME = 14.78925 hours. Multiply by 15 to get degrees: LST = 221.8388_degawNi
©緣生術數研究社 -- 術數研究 I:x~xd
The Moon's Hour Angle becomes:>
©緣生術數研究社 -- 術數研究 vcem
HA = LST - RA = -87.6623_deg = 272.3377_degU
©緣生術數研究社 -- 術數研究 M>L-as
Our latitude +60_deg yields the following geocentric latitude and distance from the Earth's center:R Zl
©緣生術數研究社 -- 術數研究 Bu
gclat = 59.83_deg rho = 0.9975YB,@
©緣生術數研究社 -- 術數研究 {
We've already computed the Moon's distance as 60.6779 Earth radii, which means the Moon's parallax is:Bt
©緣生術數研究社 -- 術數研究 ]s
mpar = 0.9443_deg;8@(
©緣生術數研究社 -- 術數研究 K6\%
The auxiliary angle g becomes::
©緣生術數研究社 -- 術數研究 "`>qv%
g = 88.642s(V3
©緣生術數研究社 -- 術數研究 ww*[
And finally the Moon's topocentric position becomes:MBp>H
©緣生術數研究社 -- 術數研究 dvHD[&
topRA = 309.5011_deg - (-0.5006_deg) = 310.0017_deg9]WV
topDecl = -19.1032_deg - (+0.7758_deg) = -19.8790_deg`U!`?
©緣生術數研究社 -- 術數研究 h@'
This correction to topocentric position can also be applied to the Sun and the planets. But since they're much farther away, the correction becomes much smaller. It's largest for Venus at inferior conjunction, when Venus' parallax is somewhat larger than 32 arc seconds. Within our aim of obtaining a final accuracy of 1-2 arc minutes, it might barely be justified to correct to topocentric position when Venus is close to inferior conjunction, and perhaps also when Mars is at a favourable opposition. But in all other cases this correction can safely be ignored within our accuracy aim. We only need to worry about the Moon in this case.a6zkm
©緣生術數研究社 -- 術數研究 ,3,O
If you want to compute topocentric coordinates for the planets anyway, you do it the same way as for the Moon, with one exception: the parallax of the planet (ppar) is computed from this formula:Wj
©緣生術數研究社 -- 術數研究 ^)~ET
ppar = (8.794/3600)_deg / rnO[0r
©緣生術數研究社 -- 術數研究 V`Z
where r is the distance of the planet from the Earth, in astronomical units.s;5q
©緣生術數研究社 -- 術數研究 TQ7N(
10. The orbital elements of the planets{;uk
To compute the positions of the major planets, we first must compute their orbital elements:ZSXWj
©緣生術數研究社 -- 術數研究 2qv0
Mercury:+0\A
©緣生術數研究社 -- 術數研究 ]!8T^{
N = 48.3313_deg + 3.24587E-5_deg * d (Long of asc. node)OR:LA
i = 7.0047_deg + 5.00E-8_deg * d (Inclination)D
w = 29.1241_deg + 1.01444E-5_deg * d (Argument of perihelion)"'<M
a = 0.387098 (Semi-major axis)$nV
e = 0.205635 + 5.59E-10 * d (Eccentricity)['
M = 168.6562_deg + 4.0923344368_deg * d (Mean anonaly)fs
©緣生術數研究社 -- 術數研究 -"WC
Venus:=a[ofO
©緣生術數研究社 -- 術數研究 lV
N = 76.6799_deg + 2.46590E-5_deg * dOz
i = 3.3946_deg + 2.75E-8_deg * df^
w = 54.8910_deg + 1.38374E-5_deg * dmRo_zc
a = 0.723330T/
e = 0.006773 - 1.302E-9 * d3UN9U,
M = 48.0052_deg + 1.6021302244_deg * d le$,
©緣生術數研究社 -- 術數研究 Uu
Mars:©緣生術數研究社 -- 術數研究 <`t9l~
©緣生術數研究社 -- 術數研究 :T
N = 49.5574_deg + 2.11081E-5_deg * d
i = 1.8497_deg - 1.78E-8_deg * dnIAOa9
w = 286.5016_deg + 2.92961E-5_deg * dkAy
a = 1.523688j
e = 0.093405 + 2.516E-9 * dwZ$Xb
M = 18.6021_deg + 0.5240207766_deg * dzT{D@
©緣生術數研究社 -- 術數研究 ZR-B
Jupiter:R">*
©緣生術數研究社 -- 術數研究 VHz
N = 100.4542_deg + 2.76854E-5_deg * d&s6
i = 1.3030_deg - 1.557E-7_deg * dM^
w = 273.8777_deg + 1.64505E-5_deg * dPb-n
a = 5.20256URM0
e = 0.048498 + 4.469E-9 * djrGbp&
M = 19.8950_deg + 0.0830853001_deg * d9A
©緣生術數研究社 -- 術數研究 )VAi{
Saturn:sH5-
©緣生術數研究社 -- 術數研究 `Q
N = 113.6634_deg + 2.38980E-5_deg * d_3h
i = 2.4886_deg - 1.081E-7_deg * dsn
w = 339.3939_deg + 2.97661E-5_deg * duz
a = 9.55475-bj
e = 0.055546 - 9.499E-9 * dTS
M = 316.9670_deg + 0.0334442282_deg * dnTS/
©緣生術數研究社 -- 術數研究 L#i3
Uranus:O9)[t
©緣生術數研究社 -- 術數研究 [
N = 74.0005_deg + 1.3978E-5_deg * dr&
i = 0.7733_deg + 1.9E-8_deg * d(g
w = 96.6612_deg + 3.0565E-5_deg * dJuIpTV
a = 19.18171 - 1.55E-8 * dnQS<9c
e = 0.047318 + 7.45E-9 * dY6!&
M = 142.5905_deg + 0.011725806_deg * dR%Y
©緣生術數研究社 -- 術數研究 %Y4k}n
Neptune:c-,&`M
©緣生術數研究社 -- 術數研究 Nmr
N = 131.7806_deg + 3.0173E-5_deg * d!s
i = 1.7700_deg - 2.55E-7_deg * dC
w = 272.8461_deg - 6.027E-6_deg * dnY
a = 30.05826 + 3.313E-8 * d*[i[8
e = 0.008606 + 2.15E-9 * dUA
M = 260.2471_deg + 0.005995147_deg * d+"*w!
©緣生術數研究社 -- 術數研究 3LPX
Let's compute all these elements for our test date, 19 april 1990 0h UT:i~7b~
©緣生術數研究社 -- 術數研究 |h@)
N i w a e MEX\8}
deg deg deg a.e. degGC.S=
©緣生術數研究社 -- 術數研究 '`J
Mercury 48.2163 7.0045 29.0882 0.387098 0.205633 69.5153#
Venus 76.5925 3.3945 54.8420 0.723330 0.006778 131.6578hBTm_
Mars 49.4826 1.8498 286.3978 1.523688 0.093396 321.9965)sz8qz
©緣生術數研究社 -- 術數研究 g
Jupiter 100.3561 1.3036 273.8194 5.20256 0.048482 85.5238N
Saturn 113.5787 2.4890 339.2884 9.55475 0.055580 198.4741}9
Uranus 73.9510 0.7732 96.5529 19.18176 0.047292 101.0460^RWV=/
Neptune 131.6737 1.7709 272.8675 30.05814 0.008598 239.0063 ,R
©緣生術數研究社 -- 術數研究 +?
11. The heliocentric positions of the planets%H5<qa
The heliocentric ecliptic coordinates of the planets are computed in the same way as we computed the geocentric ecliptic coordinates of the Moon: first we compute E, the eccentric anomaly. Several planetary orbits have quite high eccentricities, which means we must use the iteration formula to obtain an accurate value of E. When we know E, we compute, as earlier, the distance r ("radius vector") and the true anomaly, v. Then we compute ecliptic rectangular (x,y,z) coordinates as we did for the Moon. Since the Moon orbits the Earth, this position of the Moon was geocentric. The planets though orbit the Sun, which means we get heliocentric positions instead. Also the semi-major axis, a, and the distance, r, which was given in Earth radii for the Moon, are given in astronomical units for the planets, where one astronomical unit is 149.6 million kilometers.!`gZ9
©緣生術數研究社 -- 術數研究 !'iN_H
Let's do this for Mercury on our test date: the first approximation of E is 81.3464_deg, and following iterations yield 81.1572_deg, 81.1572_deg .... From this we find:z"X_D
©緣生術數研究社 -- 術數研究 {yA8
r = 0.374862l
v = 93.0727_deg`{Bsfr
©緣生術數研究社 -- 術數研究 K
Mercury's heliocentric ecliptic rectangular coordinates become:1E
©緣生術數研究社 -- 術數研究 ,Uvp
x = -0.367821Ekb
y = +0.061084L
z = +0.038699z$>-
©緣生術數研究社 -- 術數研究 JJ0Wa^
Convert to spherical coordinates:@7<
©緣生術數研究社 -- 術數研究 -m
lon = 170.5709_deg9~*
lat = +5.9255_degc3
r = 0.374862~I%A>^
©緣生術數研究社 -- 術數研究 %;l
The Astronomical Almanac gives these figures:)7
©緣生術數研究社 -- 術數研究 HiJ}wG
lon = 170.5701_degyq
lat = +5.9258_degkL4]_
r = 0.374856z/?D
©緣生術數研究社 -- 術數研究 _35x
The agreement is almost perfect! The discrepancy is only a few arc seconds. This is because it's quite easy to get an accurate position for Mercury: it's close to the Sun where the Sun's gravitational field is strongest, and therefore perturbations from the other planets will be smallest for Mercury.6KY;
©緣生術數研究社 -- 術數研究 80FO
If we compute the heliocentric longitude, latitude and distance for the other planets from their orbital elements, we get:8D
©緣生術數研究社 -- 術數研究 kf
Heliocentric`j
longitude latitude distance"qs!
lon lat r4hwD:
©緣生術數研究社 -- 術數研究 `ky
Mercury 170.5709_deg +5.9255_deg 0.374862bj
Venus 263.6570_deg -0.4180_deg 0.726607kZ[
Mars 290.6297_deg -1.6203_deg 1.417194jdG
Jupiter 105.2543_deg +0.1113_deg 5.19508[Y9{g
Saturn 289.4523_deg +0.1792_deg 10.06118%l8le_
Uranus 276.7999_deg -0.3003_deg 19.396282XBh
Neptune 282.7192_deg +0.8575_deg 30.19284v
©緣生術數研究社 -- 術數研究 ;_XO
For e.g. Saturn, the Astronomical Almanac says:@wspk
©緣生術數研究社 -- 術數研究 0
lon = 289.3864_degKb
lat = +0.1816_deg65)I ?
r = 10.018506~-mzm
©緣生術數研究社 -- 術數研究 J/i
The difference is here much larger! For Mercury our discrepancy was only a few arc seconds, but for Saturn it's up to four arc minutes! And still we've been lucky, since sometimes the discrepancy can be up to one full degree for Saturn. This is the planet that's perturbed most severely, mostly by Jupiter.5X'
©緣生術數研究社 -- 術數研究 '58V6#
12. Higher accuracy - perturbations}
To reach our aim of a final accuracy of 1-2 arc minutes, we must compute Jupiter's and Saturn's perturbations on each other, and their perturbations on Uranus. The perturbations on, and by, other planets can be ignored, with our aim for 1-2 arcmin accuracy.F6>;u
©緣生術數研究社 -- 術數研究 B
First we need three fundamental arguments:`%D3
©緣生術數研究社 -- 術數研究 @-o Qy
Jupiters mean anomaly: Mjn&j&
Saturn mean anomaly: Ms H
Uranus mean anomaly: Mu00
©緣生術數研究社 -- 術數研究 C.q
Then these terms should be added to Jupiter's heliocentric longitude:]7])-
©緣生術數研究社 -- 術數研究 [)U[f>
-0.332_deg * sin(2*Mj - 5*Ms - 67.6_deg):/
-0.056_deg * sin(2*Mj - 2*Ms + 21_deg)x*
+0.042_deg * sin(3*Mj - 5*Ms + 21_deg)6Dg
-0.036_deg * sin(Mj - 2*Ms)<+5f
+0.022_deg * cos(Mj - Ms)1X4p
+0.023_deg * sin(2*Mj - 3*Ms + 52_deg)j(ZvXE
-0.016_deg * sin(Mj - 5*Ms - 69_deg)K
©緣生術數研究社 -- 術數研究 QmJZ?
For Saturn we must correct both the longitude and latitude. Add this to Saturn's heliocentric longitude:ok
©緣生術數研究社 -- 術數研究 @
+0.812_deg * sin(2*Mj - 5*Ms - 67.6_deg)w/t
-0.229_deg * cos(2*Mj - 4*Ms - 2_deg)lSP
+0.119_deg * sin(Mj - 2*Ms - 3_deg)tCIBR
+0.046_deg * sin(2*Mj - 6*Ms - 69_deg)r#
+0.014_deg * sin(Mj - 3*Ms + 32_deg)2\D
©緣生術數研究社 -- 術數研究 (~:
and to Saturn's heliocentric latitude these terms should be added::"
©緣生術數研究社 -- 術數研究 3]}nIn
-0.020_deg * cos(2*Mj - 4*Ms - 2_deg)O7
+0.018_deg * sin(2*Mj - 6*Ms - 49_deg)Z
©緣生術數研究社 -- 術數研究 \Cw
Finally, add this to Uranus heliocentric longitude:xX0
©緣生術數研究社 -- 術數研究 Q5z50j
+0.040_deg * sin(Ms - 2*Mu + 6_deg)hU)
+0.035_deg * sin(Ms - 3*Mu + 33_deg)Q
-0.015_deg * sin(Mj - Mu + 20_deg)B~F4=.
©緣生術數研究社 -- 術數研究 "u>+
The perturbation terms above are all terms having an amplitude of 0.01 degrees or more. We ignore all perturbations in the distances of the planets, since a modest perturbation in distance won't affect the apparent position very much.*6'iOu
©緣生術數研究社 -- 術數研究 7OsY
The largest perturbation term, "the grand Jupiter-Saturn term" is the perturbation in longitude with the largest amplitude in both Jupiter and Saturn. Its period is 918 years, and its amplitude is a large part of a degree for both Jupiter and Saturn. There is also a "grand Saturn-Uranus term", which has a period of 560 years and an amplitude of 0.035 degrees for Uranus, but less than 0.01 degrees for Saturn. Other included terms have periods between 14 and 100 years. Finally we should mention the "grand Uranus-Neptune term", which has a period if 4200 years and an amplitude of almost one degree. It's not included in our perturbation terms here, instead its effects have been included in the orbital elements for Uranus and Neptune. This is why the mean distances of Uranus and Neptune are varying.Lm(
©緣生術數研究社 -- 術數研究 haw-b
If we compute the perturbations for our test date, we get:a0gH
©緣生術數研究社 -- 術數研究 a!12/?
Mj = 85.5238_deg Ms = 198.4741_deg Mu = 101.0460:5
©緣生術數研究社 -- 術數研究 \fQ
Perturbations in Jupiter's longitude:4:[f
©緣生術數研究社 -- 術數研究 W9
+ 0.0637_deg - 0.0236_deg + 0.0038_deg - 0.0270_deg - 0.0086_degd1L-3G
- 0.0049_deg - 0.0155_deg = -0.0120_deg6
©緣生術數研究社 -- 術數研究 N3T
Jupiter's heliocentric longitude, with perturbations: 105.2423_deg3Cx>
The Astronomical Almanac says: 105.2603_deg[qc3
©緣生術數研究社 -- 術數研究 8E
Perturbations in Saturn's longitude:H-'s
©緣生術數研究社 -- 術數研究 YAdU
-0.1560_deg + 0.0206_deg + 0.0850_deg - 0.0070_deg - 0.0124_deg
= - 0.0699_deg1e*
©緣生術數研究社 -- 術數研究 X
Perturbations in Saturn's latitude:Tc$Nep
©緣生術數研究社 -- 術數研究 Lkw\X"
+0.0018_deg + 0.0035_deg = +0.0053_degD
©緣生術數研究社 -- 術數研究 ,-
Saturns position, with perturbations: 289.3824_deg +0.1845_deg.+l!(q
The Astronomical Almanac says: 289.3864_deg +0.1816_deg(N2M^_
©緣生術數研究社 -- 術數研究 !D!AQn
Perturbations in Uranus' longitude:-Ka6<2
©緣生術數研究社 -- 術數研究 \_GQ/
+0.0017_deg - 0.0332_deg - 0.0012_deg = -0.0327_deg/
©緣生術數研究社 -- 術數研究 9wQ)Mg
Uranus heliocentric longitude, with perturbations: 276.7672_deg=q)7
The Astronomical Almanac says: 276.7706_degtp/*Ij
©緣生術數研究社 -- 術數研究 3<?AY
13. Precessionw
The planetary positions computed here are for "the epoch of the day", i.e. relative to the celestial equator and ecliptic at the moment. Sometimes you need to use some other epoch, e.g. some standard epoch like 1950.0 or 2000.0. Due to our modest accuracy requirement of 1-2 arc minutes, we need not distinguish J2000.0 from B2000.0, it's enough to simply use 2000.0.+VQXI
©緣生術數研究社 -- 術數研究 T-nzx
We will simplify the precession correction further by doing it in eliptic coordinates: the correction is simply done by addingj<$
©緣生術數研究社 -- 術數研究 D[&"
3.82394E-5_deg * ( 365.2422 * ( epoch - 2000.0 ) - d )]=d^
©緣生術數研究社 -- 術數研究 8
to the ecliptic longitude. We ignore precession in ecliptic latitude. "epoch" is the epoch we wish to precess to, and "d" is the "day number" we used when computing our planetary positions.f~
©緣生術數研究社 -- 術數研究 xgg]\\
Example: if we wish to precess computations done at our test date 19 April 1990, when d = -3543, we add the quantity below (degrees) to the ecliptic longitude:Vfd
©緣生術數研究社 -- 術數研究 (55R
3.82394E-5_deg * ( 365.2422 * ( 2000.0 - 2000.0 ) - (-3543) ) =KTB?v6
= 0.1355_deg0
©緣生術數研究社 -- 術數研究 nzCrE
So we simply add 0.1355_deg to our ecliptic longitude to get the position at 2000.0.$?F
©緣生術數研究社 -- 術數研究 xhp}
14. Geocentric positions of the planetsHQT7
To convert the planets' heliocentric positions to geocentric positions, we simply add the Sun's rectangular (x,y,z) coordinates to the rectangular (x,y,z) heliocentric coordinates of the planet:\
©緣生術數研究社 -- 術數研究 T
Let's do this for Mercury on our test date - we add the x, y and z coordinates separately:a
©緣生術數研究社 -- 術數研究 :$kq4i
xsun = +0.881048 ysun = +0.482098 zsun = 0.0HL
xplan = -0.367821 yplan = +0.061084 zplan = +0.038699jw9C#
-----------------------------------------------------------------l1]^m
xgeoc = +0.513227 ygeoc = +0.543182 zgeoc = +0.038699nB}R6
©緣生術數研究社 -- 術數研究 0Y2g`
Now we have rectangular geocentric coordinates of Mercury. If we wish, we can convert this to spherical coordinates - then we get geocentric ecliptic longitude and latitude. This is useful if we want to precess the position to some other epoch: we then simply add the approproate precessional value to the longitude. Then we can convert back to rectangular coordinates.}CyJE:
©緣生術數研究社 -- 術數研究 X(.!l
But for the moment we want the "epoch of the day": let's rotate the x,y,z, coordinates around the X axis, as described earlier. Then we'll get equatorial rectangular geocentric (whew!) coordinates:v.c
©緣生術數研究社 -- 術數研究 h,
xequat = +0.513227 yequat = +0.482961 zequat = 0.251582C[Dh}
©緣生術數研究社 -- 術數研究 <
We can convert these coordinates to spherical coordinates, and then we'll (finally!) get geocentric Right Ascension, Declination and distance for Mercury:+
©緣生術數研究社 -- 術數研究 $* {b
RA = 43.2598_deg Decl = +19.6460_deg R = 0.748296W~
©緣生術數研究社 -- 術數研究 )4vyQ
Note that the distance now is different. This is quite natural since the distance now is from the Earth and not, as earlier, from the Sun.b>9zIG
©緣生術數研究社 -- 術數研究 yL
Let's conclude by checking the values given by the Astronomical Almanac: H<2+p
©緣生術數研究社 -- 術數研究 :H
RA = 43.2535_deg Decl = +19.6458_deg R = 0.7482628"
©緣生術數研究社 -- 術數研究 YT|tb
15. The elongation and physical ephemerides of the planetsS4`8
When we finally have completed our computation of the heliocentric and geocentric coordinates of the planets, it could also be interesting to know what the planet will look like. How large will it appear? What are its phase and magnitude (brightness)? These computations are much simpler than the computations of the positions.oCG:R
©緣生術數研究社 -- 術數研究 E'2@\Z
Let's start by computing the apparent diameter of the planet:e
©緣生術數研究社 -- 術數研究 e``xVy
d = d0 / R`GC!
©緣生術數研究社 -- 術數研究 Ce{x
R is the planet's geocentric distance in astronomical units, and d is the planet's apparent diameter at a distance of 1 astronomical unit. d0 is of course different for each planet. The values below are given in seconds of arc. Some planets have different equatorial and polar diameters:n:j@
©緣生術數研究社 -- 術數研究 _fkVO
Mercury 6.74"@rh
Venus 16.92"DuEn/
Earth 17.59" equ 17.53" polI9
Mars 9.36" equ 9.28" polK6jy
Jupiter 196.94" equ 185.08" polQU3#8"
Saturn 165.6" equ 150.8" polu9L'^N
Uranus 65.8" equ 62.1" polRqIV
Neptune 62.2" equ 60.9" poll6
©緣生術數研究社 -- 術數研究 +/@
The Sun's apparent diameter at 1 astronomical unit is 1919.26". The Moon's apparent diameter is:{PVI=
©緣生術數研究社 -- 術數研究 bLm
d = 1873.7" * 60 / rC
©緣生術數研究社 -- 術數研究 WbX|^e
where r is the Moon's distance in Earth radii.,Yh{
©緣生術數研究社 -- 術數研究 ]}(
Two other quantities we'd like to know are the phase angle and the elongation.S2j=_
©緣生術數研究社 -- 術數研究 e9"d5
The phase angle tells us the phase: if it's zero the planet appears "full", if it's 90 degrees it appears "half", and if it's 180 degrees it appears "new". Only the Moon and the inferior planets (i.e. Mercury and Venus) can have phase angles exceeding about 50 degrees.ZS3t
©緣生術數研究社 -- 術數研究 {ey_+p
The elongation is the apparent angular distance of the planet from the Sun. If the elongation is smaller than ca 20 degrees, the planet is hard to observe, and if it's smaller than ca 10 degrees it's usually not possible to observe the planet.u|
©緣生術數研究社 -- 術數研究 mc
To compute phase angle and elongation we need to know the planet's heliocentric distance, r, its geocentric distance, R, and the distance to the Sun, s. Now we can compute the phase angle, FV, and the elongation, elong:8Do4(
©緣生術數研究社 -- 術數研究 i
elong = acos( ( s*s + R*R - r*r ) / (2*s*R) )?)gf>B
©緣生術數研究社 -- 術數研究 o
FV = acos( ( r*r + R*R - s*s ) / (2*r*R) )Pj+01a
©緣生術數研究社 -- 術數研究 :T
When we know the phase angle, we can easily compute the phase:;ge.k=
©緣生術數研究社 -- 術數研究 $-at
phase = ( 1 + cos(FV) ) / 2 = hav(180_deg - FV)Z5w
©緣生術數研究社 -- 術數研究 &CuSu<
hav(FV) is the "haversine" of the phase angle. The "haversine" (or "half versine") is an old and now obsolete trigonometric function; it's defined as:JAG
©緣生術數研究社 -- 術數研究 _BD)~
hav(x) = ( 1 - cos(x) ) / 2 = sin^2 (x/2)l7L8sD
©緣生術數研究社 -- 術數研究 h("
As usual we must use a different procedure for the Moon. Since the Moon is so close to the Earth, the procedure above would introduce too big errors. Instead we use the Moon's ecliptic longitude and latitude, mlon and mlat, and the Sun's ecliptic longitude, mlon, to compute first the elongation, then the phase angle, of the Moon:]ub"
©緣生術數研究社 -- 術數研究 EN
elong = acos( cos(slon - mlon) * cos(mlat) ),.y
g;,I
FV = 180_deg - elongZ
©緣生術數研究社 -- 術數研究 ;k
Finally we'll compute the magnitude (or brightness) of the planets. Here we need to use a formula that's different for each planet. The phase angle, FV, is in degrees:}
©緣生術數研究社 -- 術數研究 7zZ:-h
Mercury: -0.36 + 5*log10(r*R) + 0.027 * FV + 2.2E-13 * FV**6CP
Venus: -4.34 + 5*log10(r*R) + 0.013 * FV + 4.2E-7 * FV**32h1
Mars: -1.51 + 5*log10(r*R) + 0.016 * FV!7B
Jupiter: -9.25 + 5*log10(r*R) + 0.014 * FV#OS^
Saturn: -9.0 + 5*log10(r*R) + 0.044 * FV + ring_magn{&tr
Uranus: -7.15 + 5*log10(r*R) + 0.001 * FVG{x
Neptune: -6.90 + 5*log10(r*R) + 0.001 * FVs-he$
©緣生術數研究社 -- 術數研究 oAhI
** is the power operator, thus FV**6 is the phase angle (in degrees) raised to the sixth power. If FV is 150 degrees, then FV**6 becomes ca 1.14E+13, which is a quite large number.V8y
©緣生術數研究社 -- 術數研究 m>
Saturn needs special treatment due to its rings: when Saturn's rings are "open" then Saturn will appear much brighter than when we view Saturn's rings edgewise. We'll compute ring_mang like this:!LCkv
©緣生術數研究社 -- 術數研究 5F
ring_magn = -2.6 * sin(abs(B)) + 1.2 * (sin(B))**20
©緣生術數研究社 -- 術數研究 d_pE
Here B is the tilt of Saturn's rings which we also need to compute. Then we start with Saturn's geocentric ecliptic longitude and latitude (los, las) which we've already computed. We also need the tilt of the rings to the ecliptic, ir, and the "ascending node" of the plane of the rings, Nr:)z
©緣生術數研究社 -- 術數研究 W)~)Y
ir = 28.06_degV
Nr = 169.51_deg + 3.82E-5_deg * d67,2@6
©緣生術數研究社 -- 術數研究 +^y~
Here d is our "day number" which we've used so many times before. For our test date d = -3543. Now let's compute the tilt of the rings:86
©緣生術數研究社 -- 術數研究 C
B = asin( sin(las) * cos(ir) - cos(las) * sin(ir) * sin(los-Nr) );Ca
©緣生術數研究社 -- 術數研究 ,$f&
This concludes our computation of the magnitudes of the planets.dTt
©緣生術數研究社 -- 術數研究 sk
16. The positions of comets. Comet Encke and Levy.]6BZ
If you want to compute the position of a comet or an asteroid, you must have access to orbital elements that still are valid. One set of orbital elements isn't valid forever. For instance if you try to use the 1986 orbital elements of comet Halley to compute its appearance in either 1910 or 2061, you'll get very large errors in your computed positions - sometimes the errors will be 90 degrees or more.lcFdNV
©緣生術數研究社 -- 術數研究 `yy
Comets will usually have a new set of orbital elements computed for each perihelion. The comets are perturbed most severely when they're close to aphelion, far away from the gravity of the Sun but maybe much closer to Jupiter, Saturn, Uranus or Neptune. When the comet is passing through the inner solar system, the perturbations are usually so small that the same set of orbital elements can be used for the entire apparition.rp5
©緣生術數研究社 -- 術數研究 Zpt
Orbital elements for an asteroid should preferably not be more than about one year old. If your accuracy requirements are lower, you can of course use older elements. If you use orbital elements that are five years old for a main-belt asteroid, then your computed positions can be several degrees in error. If the orbital elements are less than one year old, the errors usually stay below approximately one arc minute, for a main-belt asteroid.0WMD}
©緣生術數研究社 -- 術數研究 7kxBZ
If you have access to valid orbital elements for a comet or an asteroid, proceed as below to compute its position at some date:Nv
©緣生術數研究社 -- 術數研究 V
1. If necessary, precess the angular elements N,w,i to the epoch of today. The simples way to do this is to add the precession angle to N, the longitude of the ascending node. This method is approximate, but it's good enough for our accuracy aim of 1-2 arc minutes.@1
©緣生術數研究社 -- 術數研究 v<;
2. Compute the day number for the time or perihelion, call it D. Then compute the number of days since perihelion, d - D (before perihelion this number is of course negative).>[Wk-
©緣生術數研究社 -- 術數研究 sj:EN
3. If the orbit is elliptical, compute the Mean Anomaly, M. Then compute r, the heliocentric distance, and v, the true anomaly..9=|iv
©緣生術數研究社 -- 術數研究 b^
4. If the orbit is a parabola, or close to a parabola (the eccentricity is 1.0 or nearly 1.0), then the algorithms for elliptical orbits will break down. Then use another algorithm, presented below, to compute r, the heliocentric distance, and v, the true anomaly, for near-parabolic orbits.d
©緣生術數研究社 -- 術數研究 d\|@yu
5. When you know r and v, proceed as with the planets: compute first the heliocentric, then the geocentric, position.=.e;
©緣生術數研究社 -- 術數研究 LZkZxY
6. If needed, precess the final position to the desired epoch, e.g. 2000.0PI
©緣生術數研究社 -- 術數研究 _
A quantity we'll encounter here is Gauss' gravitational constant, k. This constant links the Sun's mass with our time unit (the day) and the length unit (the astronomical unit). The EXACT value of Gauss' gravitational constant k is:RE>rud
©緣生術數研究社 -- 術數研究 S
k = 0.01720209895 (exactly!)It7h
©緣生術數研究社 -- 術數研究 y
If the orbit is elliptical, and if the perihelion distance, q, is given instead of the mean distance, a, we start by computing the mean distance a from the perihelion distance q and the eccentricity e:D[
©緣生術數研究社 -- 術數研究 z%2
a = q / (1 - e);.
©緣生術數研究社 -- 術數研究 Qt90
Now we compute the Mean Anomaly, M:zA
©緣生術數研究社 -- 術數研究 jS{
M = (180_deg/pi) * (d - D) * k / (a ** 1.5)qL:R>
©緣生術數研究社 -- 術數研究 +s;N+q
a ** 1.5 is most easily computed as: sqrt(a*a*a):6huZX
©緣生術數研究社 -- 術數研究 .s
Now we know the Mean Anomaly, M. We proceed as for a planetary orbit by computing E, the eccentric anomaly. Since comet and asteroid orbits often have high eccentricities, we must use the iteration formula given earlier, and be sure to iterate until we get convergence for the value of E.8>J5
©緣生術數研究社 -- 術數研究 X"
The orbital period for a comet or an asteroid in elliptic orbit is (P in days):HRX1
©緣生術數研究社 -- 術數研究 D^>J
P = 2 * pi * (a ** 1.5) / k<'
©緣生術數研究社 -- 術數研究 '\MRM
If the comet's orbit is a parabola, the algorithm for elliptic orbits will break down: the semi-major axis and the orbital period will be infinite, and the Mean Anomaly will be zero. Then we must proceed in a different way. For a parabolic orbit we start by computing the quantities a, b and w (where a is not at all related to a for an elliptic orbit):TS`tne
©緣生術數研究社 -- 術數研究 *5
a = 1.5 * (d - D) * k / sqrt(2 * q*q*q)+94kO
©緣生術數研究社 -- 術數研究 }c
b = sqrt( 1 + a*a )XZ3
©緣生術數研究社 -- 術數研究 <iB
w = cbrt(b + a) - cbrt(b - a)P:+C
©緣生術數研究社 -- 術數研究 c
cbrt is the Cubic Root function. Finally we compute the true anomaly, v, and the heliocentric distiance, r:yjm,
©緣生術數研究社 -- 術數研究 ,
v = 2 * atan(w)HYeY
r = q * ( 1 + w*w )@Ky9e\
©緣生術數研究社 -- 術數研究 EXra
From here we can proceed as usual.Y;4~
©緣生術數研究社 -- 術數研究 M*<Y
Finally we have the case that's most common for newly discovered comets: the orbit isn't an exact parabola, but very nearly so. It's eccentricity is slightly below, or slightly above, one. The algorithm presented here can be used for eccentricities between about 0.98 and 1.02. If the eccentricity is smaller than 0.98 the elliptic algorithm should be used instead. No known comet has an eccentricity exceeding 1.02.\O
©緣生術數研究社 -- 術數研究 REztS
As for the purely parabolic orbit, we start by computing the time since perihelion in days, d - D, and the perihelion distance, q. We also need to know the eccentricity, e. Then we can proceed as:VV"
©緣生術數研究社 -- 術數研究 /qqV
a = 0.75 * (d - D) * k * sqrt( (1 + e) / (q*q*q) )
b = sqrt( 1 + a*a )*D6y\N
W = cbrt(b + a) - cbrt(b - a)3e):[O
c = 1 + 1/(W*W)oxBacD
f = (1 - e) / (1 + e)db_?L
g = f / (c*c)+>j'.
©緣生術數研究社 -- 術數研究 0ujGn<
a1 = (2/3) + (2/5) * W*W?
a2 = (7/5) + (33/35) * W*W + (37/175) * W**4kkeg
a3 = W*W * ( 432/175) + (956/1125) * W*W + (84/1575) * W**4 )8r
©緣生術數研究社 -- 術數研究 E\Fg
w = W * ( 1 + g * c * ( a1 + a2*g + a3*g*g ) )~
©緣生術數研究社 -- 術數研究 >;QSn
v = 2 * atan(w)%{
r = q * ( 1 + w*w ) / ( 1 + w*w * f )#cC
©緣生術數研究社 -- 術數研究 Rzx4$
This algorithm yields the true anomaly, v, and the heliocentric distance, r, for a nearly-parabolic orbit.Vu]z&W
©緣生術數研究社 -- 術數研究 E
Now it's time for a practical example. Let's select two of the comets that were seen in the autumn of 1990: Comet Encke, a well-known periodic comet, and comet Levy, which was easily seen towards a dark sky in the autumn of 1990. When passing the inner solar system, the orbit of comet Levy was slightly hyperbolic.[xp6BZ
©緣生術數研究社 -- 術數研究 [@
According the the Handbook of the British Astronomical Association the orbital elements for comet Encke in 1990 are:vdYGG
©緣生術數研究社 -- 術數研究 "W2T0{
T = 1990 Oct 28.54502 TDT6|
e = 0.8502196)`NCZ
q = 0.3308858{Y8
w = 186.24444_deg@
N = 334.04096_deg 1950.0iK%hM,
i = 11.93911_deg|9
©緣生術數研究社 -- 術數研究 z;PT@|
The orbital elements for comet Levy are (BAA Circular 704):_NL
©緣生術數研究社 -- 術數研究 MdO+?y
T = 1990 Oct 24.6954 ET?36
e = 1.00027076S
q = 0.938586k&
w = 242.6797_deg!d7@
N = 138.6637_deg 1950.06'g0
i = 131.5856_degc'b/
©緣生術數研究社 -- 術數研究 wq^L
Let's also choose another test date, when both these comets were visible: 1990 Aug 22, 0t UT, which yields a "day number" d = -3418.0p[[:
©緣生術數研究社 -- 術數研究 rG`
Now we compute the day numbers at perihelion for these two comets. We get for comet Encke:a=b3,'
©緣生術數研究社 -- 術數研究 6,r
D = -3350.45498 d - D = -67.54502W@l
©緣生術數研究社 -- 術數研究 *xYI
and for comet Levy:97>wrV
©緣生術數研究社 -- 術數研究 K
D = -3354.3046 d - D = -63.6954t=
©緣生術數研究社 -- 術數研究 nHpJ
We'll continue by computing the Mean Anomaly for comet Encke:pN
©緣生術數研究社 -- 術數研究 J
M = -20.2751_deg = 339.7249_deg\
©緣生術數研究社 -- 術數研究 tT
The first approximation plus successive approximation for the Eccentric anomaly, E, becomes (degrees):mP5dm
©緣生術數研究社 -- 術數研究 $W
E = 309.3811 293.5105 295.8474 295.9061 295.9061_deg ....B[
©緣生術數研究社 -- 術數研究 {C
Here we clearly see the great need for iteration: the initial approximation differs from the final value by 14 degrees. Finally we compute the true anomaly, v, and heliocentric distance, r, for comet Encke:*q[
©緣生術數研究社 -- 術數研究 2>hW
v = 228.8837_deg9
r = 1.3885g-
©緣生術數研究社 -- 術數研究 \9EY
Now it's time for comet Levy: we'll compute the true anomaly, v, and the heliocentric distance, r, for Levy in two different ways. First we'll pretend that the orbit of Levy is an exact parabola. We get:ZFUmoQ
©緣生術數研究社 -- 術數研究 \
a = -1.2780823 b = 1.6228045 w = -0.7250189G-j
©緣生術數研究社 -- 術數研究 ?Y;
v = -71.8856_degj
r = 1.431947~C
©緣生術數研究社 -- 術數研究 SZv+P
Then we repeat the computation but accounts for the fact that Levy's orbit deviates slightly from a parabola. We get:I*nU
©緣生術數研究社 -- 術數研究 O}?4*d
a = -1.2781686 b = 1.6228724 W = -0.7250566j6mA=.
c = 2.9022000 f = -1.3498E-4 g = -1.60258E-5b<q
a1= 0.8769495 a2= 1.9540987 a3= 1.54957024IpY+C
w = -0.7250270gm
©緣生術數研究社 -- 術數研究 c2-
v = -71.8863_degPo
r = 1.432059Wa#
©緣生術數研究社 -- 術數研究 C@@K&
The difference is small in this case - only 0.0007 degrees or 2.5 arc seconds in true anomaly, and 0.000112 a.u. in heliocentric distance. Here it would have been sufficient to treat Levy's orbit as an exact parabola.IMu
©緣生術數研究社 -- 術數研究 X7)saW
Now we know the true anomaly, v, and the heliocentric distance, r, for both Encke and Levy. Next we proceed by precessing N, the longitude of the ascending node, from 1950.0 to the "epoch of the day". Let's compute the precession angle from 1950.0 to 1990 Aug 22:+
©緣生術數研究社 -- 術數研究 .3JY
prec = 3.82394E-5_deg * ( 365.2422 * ( 1950.0 - 2000.0 ) - (-3418) )pd0ib,
prec = -0.5676_degf
©緣生術數研究社 -- 術數研究 T]vwP_
To precess from 1990 Aug 22 to 1950.0, we should add this angle to N. But now we want to do the opposite: precess from 1950.0 to 1990 Aug 22, therefore we must instead subtract this angle:Ca:h
©緣生術數研究社 -- 術數研究 SB7
For comet Encke we get:{:GJ
©緣生術數研究社 -- 術數研究 DRM
N = 334.04096_deg - (-0.5676_deg) = 334.60856_degK7
©緣生術數研究社 -- 術數研究 Aa)h]Q
and for comet Levy we get:E6_^a'
©緣生術數研究社 -- 術數研究 O{y
N = 138.6637_deg - (-0.5676_deg) = 139.2313_degmD`(i
©緣生術數研究社 -- 術數研究 :Ms
Using this modified value for N we proceed just like for the planets. I won't repeat the details, but merely state some intermediate and final results:Q~S'
©緣生術數研究社 -- 術數研究 xGRM
Sun's position: x = -0.863890 y = +0.526123+f.=t_
©緣生術數研究社 -- 術數研究 bBez6<
Heliocentric: Encke Levy$C
©緣生術數研究社 -- 術數研究 x-.h
x +1.195087 +1.169908n+BG}
y +0.666455 -0.807922>
z +0.235663 +0.171375">^$
©緣生術數研究社 -- 術數研究 ^XWN(G
Geoc., eclipt.: Encke Levy*)y[o
©緣生術數研究社 -- 術數研究 0)c
x +0.331197 +0.306018\&|Y
y +1.192579 -0.281799hmp2x
z +0.235663 +0.171375pH%D"
©緣生術數研究社 -- 術數研究 RQ<
Geoc., equat.: Encke LevyIV=}
©緣生術數研究社 -- 術數研究 M*9.
x +0.331197 +0.306018wur
y +1.000414 -0.326716@o;qt5
z +0.690619 +0.045133X<\%d
©緣生術數研究社 -- 術數研究 |8
RA 71.6824_deg 313.1264_degt*O$
Decl +33.2390_deg +5.7572_deg4t{
R 1.259950 0.449919]+CL
©緣生術數研究社 -- 術數研究 K
These positions are for the "epoch of the day". If you want positions for some standard epoch, e.g. 2000.0, these positions must be precessed to that epoch.u,o^s
©緣生術數研究社 -- 術數研究 AfI{&{
Finally some notes about computing the magnitude of a comet. To accurately predict a comet's magnitude is usually hard and sometimes impossible. It's fairly common that a magnitude prediction is off by 1-2 magnitudes or even more. For comet Levy the magnitude formula looked like this:*^N
©緣生術數研究社 -- 術數研究 k^
m = 4.0 * 5*log10(R) + 10*log10(r)w
©緣生術數研究社 -- 術數研究 tC\aj
where R is the geocentric distance and r the heliocentric distance. The general case is:9~/"
©緣生術數研究社 -- 術數研究 Or)ChP
m = G * 5*log10(R) + H*log10(r)h@`Yi
©緣生術數研究社 -- 術數研究 6(O
where H usually is around 10. If H is unknown, it's usually assumed to be 10. Each comet has it's own G and H.?+wGc/
©緣生術數研究社 -- 術數研究 f,H
Some comets have a different magnitude formula. One good example is comet Encke, where the magnitude formula looks like this:@
©緣生術數研究社 -- 術數研究 n3'
m1 = 10.8 + 5*log10(R) + 3.55 * ( r**1.8 - 1 )Wo
©緣生術數研究社 -- 術數研究 N7
"m1" refers to the total magnitude of the comet. There is another cometary magnitude, "m2", which refers to the magnitude of the nucleus of the comet. The magnitude formula for Encke's m2 magnitude looks like this:e-w%&
©緣生術數研究社 -- 術數研究 JvAO
m2 = 14.5 + 5*log10(R) + 5*log10(r) + 0.03*FVdh(&
©緣生術數研究社 -- 術數研究 ]/:
Here FV is the phase angle. This kind of magnitude formula looks very much like the magnitude formula of asteroids, for a very good reason: when a comet is far away form the Sun, no gases are evaporated from the surface of the comet. Then the comet has no tail (of course) and no coma, only a nucleus. Which means the comet then behaves much like an asteroid.OsD.#
©緣生術數研究社 -- 術數研究 A
During the last few years it's become increasingly obvious that comets and asteroids often are similar kinds of solar-system objects. The asteroid (2060) Chiron has displayed cometary activity and is now also considered a comet. And in some cases comets that have "disappeared" have been re-discovered as asteroids! Apparently they "ran out of gas" and what remains of the former comet is only rock, i.e. an asteroid.;S c
| | |
|
|
|
頂端 
加到"我的最愛"
主題管理：
總固頂 
取消總固頂 
短訊息
查看
搜尋
通訊錄
引用
回覆
只看我
2004/12/21 01:56pm 
IP: 已設定保密