12 replies to this topic

[Jenab6]

I wrote this program so that I could predict where comet Ison would appear in the sky when it becomes bright at the end of 2013 and in early 2014. It will work with planets or asteroids, too.

The first line contains a calendar date in "YYYY.MMDD, hh:mm:ss" format, stored in Str 1.

The 2nd line contains the Keplerian elements of the observed object's orbit, stored in List 1.

The Keplerian elements are, in order, the semimajor axis in astronomical units, the eccentricity, the inclination to the ecliptic (degrees), the longitude of the ascending node (degrees), the argument of the perihelion (degrees), and the time of perihelion passage in Julian date format.

The semimajor axis for any hyperbolic orbit should be entered as an intrinsically negative number.

"2014.0106, 12:00:00"→Str 1↵
{-4693.4,1.0000026666,61.7852,295.7508,345.5009,2456625.3147}→List 1↵
Exp(StrLeft(Str 1,4))→Y↵
Exp(StrMid(Str 1,6,2))→M↵
Exp(StrMid(Str 1,8,2))→D↵
Exp(StrMid(Str 1,12,2))+Exp(StrMid(Str 1,15,2))÷60+Exp(StrMid(Str 1,18,2))÷3600→T↵
Int ((M−14)÷12)→A↵
Int ((1461(Y+4800+A))÷4)→B↵
Int ((367(M−2−12A))÷12)→C↵
Int ((Y+4900+A)÷100)→E↵
Int ((3E)÷4)→F↵
B+C−F+D−32075.5→J↵
J+T÷24→T¸
If List 1[1]=0 Or List 1[2]≤0 Or List 1[2]=1:Then Stop:IfEnd↵
If List 1[1]>0 And List 1[2]>1:Then Stop:IfEnd↵
If List 1[1]<0 And List 1[2]<1:Then Stop:IfEnd↵
(π÷180)List 1[3]→List 1[3]↵
(π÷180)List 1[4]→List 1[4]↵
(π÷180)List 1[5]→List 1[5]↵
If List 1[1]>0:Then Goto 2:IfEnd↵
1.32712440018ᴇ20→G↵
1.49597870691ᴇ11→O↵
86400√(G÷(-OList 1[1])^3)→M↵
M(T−List 1[6])→M↵
0→U↵
1→V↵
While Abs (V−U)>1ᴇ-12↵
U→V↵
List 1[2]sinh V−V−M→F↵
List 1[2]cosh V−1→G↵
List 1[2]sinh V→H↵
List 1[2]cosh V→I↵
-F÷G→A↵
-F÷(G+AH÷2)→B↵
-F÷(G+AH÷2+B²I÷6)→C↵
V+C→U↵
WhileEnd↵
cos⁻¹ ((List 1[2]−cosh U)÷(List 1[2]cosh U−1))→F↵
If U<0:Then 2π−F→F:IfEnd↵
List 1[1](1−List 1[2]cosh U)→R¸
Rcos F→List 2[1]↵
Rsin F→List 2[2]↵
List 2[1]cos List 1[5]−List 2[2]sin List 1[5]→List 3[1]↵
List 2[1]sin List 1[5]+List 2[2]cos List 1[5]→List 3[2]↵
List 3[1]→List 2[1]↵
List 3[2]cos List 1[3]→List 2[2]↵
List 3[2]sin List 1[3]→List 2[3]↵
List 2[1]cos List 1[4]−List 2[2]sin List 1[4]→List 3[1]↵
List 2[1]sin List 1[4]+List 2[2]cos List 1[4]→List 3[2]↵
List 2[3]→List 3[3]↵
Goto 3↵
Lbl 2↵
365.256898326List 1[1]^1.5→P↵
(T−List 1[6])÷P→M↵
M−Int M→M↵
If M<0:Then M+1→M:IfEnd↵
2πM→M↵
M→U↵
9→V↵
While Abs (V−U)>1ᴇ-12↵
U→V↵
V−List 1[2]sin V−M→F↵
1−List 1[2]cos V→G↵
List 1[2]sin V→H↵
List 1[2]cos V→I↵
-F÷G→A↵
-F÷(G+AH÷2)→B↵
-F÷(G+AH÷2+B²I÷6)→C↵
V+C→U↵
WhileEnd↵
List 1[1](cos U−List 1[2])→List 2[1]↵
List 1[1]√(1−List 1[2]²)sin U→List 2[2]↵
√(List 2[1]²+List 2[2]²)→R¸
List 2[1]cos List 1[5]−List 2[2]sin List 1[5]→List 3[1]↵
List 2[1]sin List 1[5]+List 2[2]cos List 1[5]→List 3[2]↵
List 3[1]→List 2[1]↵
List 3[2]cos List 1[3]→List 2[2]↵
List 3[2]sin List 1[3]→List 2[3]↵
List 2[1]cos List 1[4]−List 2[2]sin List 1[4]→List 3[1]↵
List 2[1]sin List 1[4]+List 2[2]cos List 1[4]→List 3[2]↵
List 2[3]→List 3[3]↵
Lbl 3↵
{1.000003,0.016701,0,0,103.154,2456294.541}→List 4↵
(π÷180)List 4[3]→List 4[3]↵
(π÷180)List 4[4]→List 4[4]↵
(π÷180)List 4[5]→List 4[5]↵
365.2563496155List 4[1]^1.5→P↵
(T−List 4[6])÷P→M↵
M−Int M→M↵
If M<0:Then M+1→M:IfEnd↵
2πM→M↵
M→U↵
9→V↵
While Abs (V−U)>1ᴇ-12↵
U→V↵
V−List 4[2]sin V−M→F↵
1−List 4[2]cos V→G↵
List 4[2]sin V→H↵
List 4[2]cos V→I↵
-F÷G→A↵
-F÷(G+AH÷2)→B↵
-F÷(G+AH÷2+B²I÷6)→C↵
V+C→U↵
WhileEnd↵
List 4[1](cos U−List 4[2])→List 2[1]↵
List 4[1]√(1−List 4[2]²)sin U→List 2[2]↵
√(List 2[1]²+List 2[2]²)→R¸
List 2[1]cos List 4[5]−List 2[2]sin List 4[5]→List 5[1]↵
List 2[1]sin List 4[5]+List 2[2]cos List 4[5]→List 5[2]↵
List 5[1]→List 2[1]↵
List 5[2]cos List 4[3]→List 2[2]↵
List 5[2]sin List 4[3]→List 2[3]↵
List 2[1]cos List 4[4]−List 2[2]sin List 4[4]→List 5[1]↵
List 2[1]sin List 4[4]+List 2[2]cos List 4[4]→List 5[2]↵
List 2[3]→List 5[3]↵
List 3[1]−List 5[1]→List 6[1]↵
List 3[2]−List 5[2]→List 6[2]↵
List 3[3]−List 5[3]→List 6[3]↵
(T−2451545)÷36525→W
84381.448−46.84024W−5.9ᴇ-4W²+1.813ᴇ-3W^3→θ↵
πθ÷648000→θ↵
List 6[1]→X↵
List 6[2]cos θ−List 6[3]sin θ→Y↵
List 6[2]sin θ+List 6[3]cos θ→Z↵
√(X²+Y²+Z²)→R↵
tan⁻¹ (Y÷X)→A↵
If X<0:Then A+π→A:IfEnd↵
If X>0 And Y<0:Then A+2π→A:IfEnd↵
(12÷π)A→A↵
(180÷π)sin⁻¹ (Z÷R)→D↵
R¸
θ¸
A¸
D¸
Locate 1,1,"t(JD)"↵
Locate 1,2,"r(AU)"↵
Locate 1,3,"R(AU)"↵
Locate 1,4,"ε(rad)"↵
Locate 1,5,"a(h)"↵
Locate 1,6,"δ(°)"↵

The output from the program are...

(1) The input date, converted to Julian Date format.
(2) The distance of the observed object from the sun in astronomical units.
(3) The distance of the observed object from Earth in astronomical units.
(4) Earth's obliquity (axis tilt) on the input date in radians.
(5) The right ascension of the object in Earth's sky in decimal hours.
(6) The declination of the object in Earth's sky in decimal degrees.

With the right ascension and declination known, you could look in a star atlas and find out which constellation the comet is in, then go outside and find the comet.

Edited by flyingfisch, 18 February 2013 - 03:46 PM.

[flyingfisch]

Hello Jenab6 and welcome to UCF! You should introduce yourself.

wow, thank you for sharing this. I tried to do something like this before but gave up long ago.

[Casimo]

This program looks very cool! Could you also make a downloadable file?

Edited by Casimo, 18 February 2013 - 03:52 PM.

[Jenab6]

Here is the downloadable file.
https://docs.google....dit?usp=sharing

There is an improved set of Keplerian orbital elements for comet Ison now. The second line in the program has been changed to

{-3618.577,1.0000034572,61.8122,295.7459,345.5039,2456625.2954}→List 1↵

Edited by Jenab6, 18 February 2013 - 07:15 PM.

[Jenab6]

Updated version. New download link.
https://docs.google....dit?usp=sharing

Changes are:

First approximation to the eccentric anomaly for elliptical orbits in the old version was the mean anomaly. Now the first approximation is found from the Danby Fourier series.

The new code retains the Danby method iteration loop to converge the eccentric anomaly, after the first approximation is found.

Mean elements for Earth on 1 January 2014 have been used to locate the observer's position. The earlier version had the elements of the Earth-moon barycenter for July 2012.

Edited by Jenab6, 20 February 2013 - 11:53 PM.

[Forty-Two]

Wow, this is really something. Where can I find the keplarian elements for a particular body?

[Jenab6]

The JPL Small Body Database is where you should get elements for the orbits of asteroids and comets.

I've written a program to return the mean Keplerian elements for the orbits of the planets. You can download MEANOE.G1M from my Google drive. This updated version of MEANOE includes the time of perihelion passage.

The user provides input for the program by editing the program and typing in the calendar date in the program's first line (replacing the example shown).

The program output is (first page), the calendar date and time (GMT), echoed from input, and the Julian date that corresponds to that calendar date and time. Then (second page) an array is shown with the planets in successive rows. The columns are (1st) semimajor axis in astronomical units, (2nd) eccentricity, (3rd) inclination to ecliptic, (4th) longitude of the ascending node, (5th) argument of the perihelion, and (6th) time of perihelion passage (Julian date).

Edited by Jenab6, 27 February 2013 - 10:49 PM.

[nsg]

Jenab,
what algorithm you are implementing in this program? I looked inside and saw sinh and cosh used. I did not expect them.
What do they do? Some kind of relativity correction?

[Jenab6]

nsg, no. There's no correction for general relativity in the COMET program. The sinh and cosh functions occur in the Danby's method iteration to converge the eccentric anomaly for hyperbolic orbits.

For ELLIPTICAL orbits, Kepler's equation relates the mean anomaly, m, and the eccentricity, e, to the eccentric anomaly, u, like this:

m = u − e sin u

Notice that the equation is transcendental in u. That means that if you know M and e, and you want to find u, you can't just use algebra in closed form. You have to use some iterative method that approaches the correct value for u by successive approximations. One such method is Newton's method. Another is Danby's method. Both of them are based on Taylor series, but whereas Newton's method uses only the first derivative of the Kepler equation, Danby's method uses the first, second, and third derivatives, which makes it converge faster.

f(u) = u − e sin u − m = 0
f'(u) = 1 − e cos u
f''(u) = e sin u
f'''(u) = e cos u

Now for HYPERBOLIC orbits, Kepler's equation is different.

m = e sinh u − u

f(u) = e sinh u − u − m = 0
f'(u) = e cosh u − 1
f''(u) = e sinh u
f'''(u) = e cosh u

The regular (spherical) trigonometry functions are used for the elliptical orbits, while the hyperbolic trig functions are used for the hyperbolic orbits because of the difference in the two kinds of curves.

[Jenab6]

Although I use Danby's method in the calculator program, once in a while it will hit an iteration that lands too close to a horizontal tangent on the function and instead of converging the procedure will "skitter," meaning that the next approximation goes way, way off. It doesn't happen often. It's a very rare problem with successive approximation methods like Newton's and Danby's. But when it does happen, there are others you can use that don't skitter.

k₀ = π/180 ≈ 0.01745329251994

n = -1

REPEAT
. n = n+1
. U₀ = 10nk₀
. U₂ = 10(n+1)k₀
. M₀ = U₀ − e sin U₀
. M₂ = U₂ − e sin U₂
. D₀ = m−M₀
. D₂ = m−M₂
UNTIL {D₀D₂≤0} or {n>35}

if n>35 then STOP {procedure failed}

REPEAT
. U₁ = (U₀+U₂)/2
. M₀ = U₀ − e sin U₀
. M₁ = U₁ − e sin U₁
. D₀ = m−M₀
. D₁ = m−M₁
. if D₀D₁<0 then U₂=U₁ else U₀=U₁
UNTIL U₂−U₀ < 1e-12

M₀ = U₀ − e sin U₀
M₂ = U₂ − e sin U₂

u = U₀ + (U₂−U₀)(m−M₀) / (M₂−M₀)


This procedure uses a bisection and reverse interpolation method to find the eccentric anomaly when the mean anomaly and the eccentricity are given. It doesn't skitter, but the code is bigger and the convergence speed is slow. It's for elliptical orbits only. I could make one for hyperbolic orbits, I guess.

Edited by Jenab6, 02 March 2013 - 05:52 AM.

[nsg]

Thanks for the explanation.
I am curious because i do not see a lot of use for hyperbolic functions outside some special applications, like power transmission calculations.
It was always a mystery to me why they include them on most of the calculators, who is a target audience.

[Jenab6]

We celestial mechanics need hyperbolic functions because some kinds of trajectories are hyperbolas.

Okay, there is an update to the comet program. (It still works for planets and asteroids, too. The object being tracked can be in either a hyperbolic orbit or an elliptical orbit.)

1. I explicitly dimensioned the List variables, so that the program would run on the fx-9750gii. Apparently, the fx-9860gii can dimension its own lists, but the fx-9750gii needs it done in the program, so I added the "dim list" code lines.

2. The program now calculates the mean elements of the Earth at date. Sorry, you can no longer put yourself on Mars to see where Vesta is, unless you add a line after the calculation of the Earth's mean elements to overwrite them in List 4. But the benefit of this is keeping your orbit for the Earth as current as the mean elements can be.

3. The program now outputs the heliocentric position (x,y,z) and distance, r, for both the comet and the Earth in ecliptic coordinates.

4. The obliquity of the ecliptic now appears in decimal degrees, not in radians.

5. The right ascension and declination are now in partitioned sexigesimal format, instead of decimal.

As before the user modifies the date in the first line of the program as his most usual change of input data.

Then, if the user wants to enter the Keplerian orbital elements of a new object, he skips down over the DIM LIST lines to reach the line for input to List 1. There can be more than one line entering the six orbital elements into List 1, and the user must be sure that the elements for the orbit he wants to run with is the last of these input lines. (Each of the lines that enter data to List 1 overwrites the one preceding it.)

Then the user runs the program.

On the first page of output appears...
* The string containing the calendar date and Greenwich Mean Time, echoed from the input on the first line of the program.
* The Julian Date corresponding to that calendar date.
* The distance of the comet from Earth on that date (astronomical units).
* The geocentric right ascension of the comet in hhmmss format (just pretend the ° is an h, the ' is a m, and the " is a s).
* The geocentric declination of the comet in ddmmss format.

On the second page of output appears...
* The obliquity of the ecliptic on the input date (in degrees).
* The angular separation (total) between the sun and the comet (in degrees).

On the third page of output appears...
* The x component of the comet's position in heliocentric ecliptic coordinates.
* The y component of the comet's position in heliocentric ecliptic coordinates.
* The z component of the comet's position in heliocentric ecliptic coordinates.
* The distance, r, between the sun and the comet on the input date (astronomical units).

On the fourth page of output appears...
* The X component of Earth's position in heliocentric ecliptic coordinates.
* The Y component of Earth's position in heliocentric ecliptic coordinates.
* The Z component of Earth's position in heliocentric ecliptic coordinates.
* The distance, R, between the sun and Earth on the input date (astronomical units).

Edited by Jenab6, 08 March 2013 - 10:42 PM.

[Rhodie]

Me=mind:blown!
#floored

And this is only a calculator...?