TheBackShed.com

Forum Index : Microcontroller and PC projects : MMX/RPI - Chebyshev ephemeris

Author

Message

cdeagle
Senior Member

Joined: 22/06/2014
Location: United States
Posts: 266

Posted: 08:36am 21 Apr 2017

Copy link to clipboard

Print this post

If you look at any of the ephemeris programs I've posted on this forum, you will notice that most of the algorithms consist of lots and lots of double precision numbers used in trig and series calculations. These calculations can involve lots of CPU time since an ephemeris algorithm will be called many times to compute events such as rise and set, eclipses and so forth.

This post describes four MicroMite BASIC program for the MMX and Raspberry Pi that can be used to create and evaluate ephemerides that consist of Chebyshev polynomials. These types of polynomials can be used to accurately represent the three components of the celestial body position vector. Also, they can be evaluated much quicker than an analytic algorithm.

Here are some guidelines for creating a Chebyshev ephemeris.

(1) don't create a Chebyshev ephemeris for a time duration longer than the orbital period of the celestial body of interest. For example, for the Moon the maximum fit duration should be 28 days, for the Earth 360 days, and so forth.

(2) the "order" of the Chebyshev fit is defined in the source code statement; nfit% = 8.

(3) the duration of the fit, in days, is defined in the line of source code which reads deltat = 30.

(4) the idea is to change the order and the fit duration until you are happy with the Chebyshev approximation. For the Moon you might use a fit of 30 while for Mars you might use 8 or perhaps 10. The motion of the Moon is quite complicated which requires a high order fit for an accurate approximation.

The four BASIC programs included in the zipped archive are

(1) make_chebyshev_moon.bas - creates lunar ephemeris data statements

(2) chebyshev_moon.bas - evaluates Chebyshev coefficients for the moon

(3) make_chebyshev_planet.bas - creates ephemeris data statements for one of the inner planets. Define the planet using the source code statement ibody% = #, where # is 1 for Mercury, 2 for Venus, 3 for Earth and 4 for Mars.

(4) chebyshev_planet.bas - evaluates Chebyshev coefficients for a planet

Within the two "make" programs, you must define the initial calendar date of the ephemeris like so

''''''''''''''''''''''''''''''
' define initial calendar date
' of the Chebyshev ephemeris
''''''''''''''''''''''''''''''

cmonth = 1

cday = 1

cyear = 2017

' initial tdb julian date

julian(cmonth, cday, cyear, jdtdb1)

PLEASE NOTE THAT THE FUNDAMENTAL TIME ARGUMENT FOR THE CHEBYSHEV EPHEMERIS, x, IS BASED ON THE BARYCENTRIC DYNAMICAL TIME (TDB) JULIAN DAY AND COMPUTED USING

x = (jdtdb - 2451545.0) / 36525.0

Here are typical DATA statements for the Moon, analytic and Chebyshev position vectors and their differences. To use this program for a user-defined Chebyshev ephemeris, simply copy the calculated DATA statements into the DATA statement area of the Chebyshev evaluation source code.

Chebyshev ephemeris of the Moon
-------------------------------

DATA 30 , 2.457754500000000e+06 , 2.457784500000000e+06

DATA -1.671212766283444e+05 , 2.339452743007602e+05 , -6.643820380953825e+04
DATA -1.013516615782387e+05 , -5.338868026854525e+04 , -3.252851865625804e+04
DATA -1.410488649048607e+05 , 3.198528042990981e+05 , -6.150080272239724e+04
DATA 2.504108815580476e+05 , 1.288933583288243e+05 , 7.969374208748670e+04
DATA 5.588685465655353e+04 , -1.387664335549331e+05 , 2.469063811957224e+04
DATA -5.357791059449854e+04 , -1.938342258785099e+04 , -1.734897073672480e+04
DATA -1.563117911126481e+03 , 1.962116581684537e+04 , -1.324979274227767e+03
DATA 6.033806923981408e+03 , -2.163099288612435e+03 , 2.132847784192905e+03
DATA -1.996995877841867e+03 , -2.076839904134528e+03 , -5.954315567404735e+02
DATA -6.101299362497032e+02 , 1.259631080869243e+03 , -2.578993245468610e+02
DATA 6.157834144601497e+02 , 1.789866309373343e+02 , 2.024828825228624e+02
DATA 1.397766147511673e+01 , -3.046738556139396e+02 , 1.667284166198618e+01
DATA -1.282736190179750e+02 , 2.496789357918427e+01 , -4.497879709269770e+01
DATA 2.756566127147638e+01 , 5.707144398157365e+01 , 7.092149271244830e+00
DATA 2.095669284160436e+01 , -2.071875123884303e+01 , 8.138655535849524e+00
DATA -1.229585583649403e+01 , -7.863442975289804e+00 , -3.633748605814500e+00
DATA -1.920784881261493e+00 , 6.959414879665635e+00 , -8.818695691720236e-01
DATA 3.437597807267711e+00 , 2.037702345105278e-01 , 9.140860024969198e-01
DATA -3.219809705323930e-01 , -1.642568770972561e+00 , -1.628938505196974e-01
DATA -7.579454218014894e-01 , 3.401845352284820e-01 , -1.184860230841291e-01
DATA 2.273044256970377e-01 , 3.331935564938542e-01 , 1.425602672792475e-01
DATA 1.755790245839424e-01 , -1.539275528470640e-01 , -2.516797414207230e-03
DATA -8.034554061915484e-02 , -8.772612029250950e-02 , -5.614954702351468e-02
DATA -5.802522909230432e-02 , 5.263075910625048e-02 , 2.241081730913251e-03
DATA 2.752260642497981e-02 , 3.340975976560355e-02 , 1.774626340756508e-02
DATA 2.174060563969936e-02 , -1.890506993438448e-02 , 1.081468033296387e-03
DATA -1.064998679231071e-02 , -1.200350847812520e-02 , -5.418672341376453e-03
DATA -6.752858228977739e-03 , 7.203210781294582e-03 , -7.713717958165123e-04
DATA 4.207065041672679e-03 , 3.157313824856593e-03 , 1.719649341764800e-03
DATA 1.498989865555608e-03 , -2.408632381223968e-03 , 1.931360317441079e-04
DATA -1.806502909828616e-03 , -4.825714913163064e-04 , -6.470521568131046e-04

Chebyshev position vector and magnitude
---------------------------------------

r_x = 353767.14756075 kilometers

r_y = 115299.02153829 kilometers

r_z = 24816.66403370 kilometers

rmag = 372908.73665646 kilometers

analytic position vector and magnitude
--------------------------------------

r_x = 353767.14704308 kilometers

r_y = 115299.02207918 kilometers

r_z = 24816.66418116 kilometers

rmag = 372908.73634242 kilometers

position vector and magnitude differences (analytic - Chebyshev)
----------------------------------------------------------------

delta r_x = -0.5177 meters

delta r_y = 0.5409 meters

delta r_z = 0.1475 meters

delta_rmag = -0.3140 meters

Here's the source code for the program that demonstrates how to evaluate a typical Chebyshev ephemeris.

' program chebyshev_moon.bas April 18, 2017

' this MMBASIC eXtreme program demonstrates how
' to evaluate a Chebyshev ephemeris of the Moon

' MicroMite eXtreme version

'''''''''''''''''''''''''''

option default float

option base 1

dim rmoon(3)

' global constants

const pi2 = 2.0 * pi, pidiv2 = 0.5 * pi, dtr = pi / 180.0, rtd = 180.0 / pi

''''''''''''''''''''''''''''''''
' read Chebyshev fit information
''''''''''''''''''''''''''''''''

read nfit%, jdtdb1, jdtdb2

dim rx_poly(nfit% + 1), ry_poly(nfit% + 1), rz_poly(nfit% + 1)

for i% = 1 to nfit% + 1

read rx_poly(i%), ry_poly(i%), rz_poly(i%)

next i%

DATA 30 , 2.457754500000000e+06 , 2.457784500000000e+06

DATA -1.671212766283444e+05 , 2.339452743007602e+05 , -6.643820380953825e+04
DATA -1.013516615782387e+05 , -5.338868026854525e+04 , -3.252851865625804e+04
DATA -1.410488649048607e+05 , 3.198528042990981e+05 , -6.150080272239724e+04
DATA 2.504108815580476e+05 , 1.288933583288243e+05 , 7.969374208748670e+04
DATA 5.588685465655353e+04 , -1.387664335549331e+05 , 2.469063811957224e+04
DATA -5.357791059449854e+04 , -1.938342258785099e+04 , -1.734897073672480e+04
DATA -1.563117911126481e+03 , 1.962116581684537e+04 , -1.324979274227767e+03
DATA 6.033806923981408e+03 , -2.163099288612435e+03 , 2.132847784192905e+03
DATA -1.996995877841867e+03 , -2.076839904134528e+03 , -5.954315567404735e+02
DATA -6.101299362497032e+02 , 1.259631080869243e+03 , -2.578993245468610e+02
DATA 6.157834144601497e+02 , 1.789866309373343e+02 , 2.024828825228624e+02
DATA 1.397766147511673e+01 , -3.046738556139396e+02 , 1.667284166198618e+01
DATA -1.282736190179750e+02 , 2.496789357918427e+01 , -4.497879709269770e+01
DATA 2.756566127147638e+01 , 5.707144398157365e+01 , 7.092149271244830e+00
DATA 2.095669284160436e+01 , -2.071875123884303e+01 , 8.138655535849524e+00
DATA -1.229585583649403e+01 , -7.863442975289804e+00 , -3.633748605814500e+00
DATA -1.920784881261493e+00 , 6.959414879665635e+00 , -8.818695691720236e-01
DATA 3.437597807267711e+00 , 2.037702345105278e-01 , 9.140860024969198e-01
DATA -3.219809705323930e-01 , -1.642568770972561e+00 , -1.628938505196974e-01
DATA -7.579454218014894e-01 , 3.401845352284820e-01 , -1.184860230841291e-01
DATA 2.273044256970377e-01 , 3.331935564938542e-01 , 1.425602672792475e-01
DATA 1.755790245839424e-01 , -1.539275528470640e-01 , -2.516797414207230e-03
DATA -8.034554061915484e-02 , -8.772612029250950e-02 , -5.614954702351468e-02
DATA -5.802522909230432e-02 , 5.263075910625048e-02 , 2.241081730913251e-03
DATA 2.752260642497981e-02 , 3.340975976560355e-02 , 1.774626340756508e-02
DATA 2.174060563969936e-02 , -1.890506993438448e-02 , 1.081468033296387e-03
DATA -1.064998679231071e-02 , -1.200350847812520e-02 , -5.418672341376453e-03
DATA -6.752858228977739e-03 , 7.203210781294582e-03 , -7.713717958165123e-04
DATA 4.207065041672679e-03 , 3.157313824856593e-03 , 1.719649341764800e-03
DATA 1.498989865555608e-03 , -2.408632381223968e-03 , 1.931360317441079e-04
DATA -1.806502909828616e-03 , -4.825714913163064e-04 , -6.470521568131046e-04

' beginning and end of valid time interval
' in Julian centuries wrt 1/1.5/2000 TDB

ta = (jdtdb1 - 2451545.0) / 36525.0

tb = (jdtdb2 - 2451545.0) / 36525.0

''''''''''''''''''''''''''''''''
' evaluate ephemeris on 1/6/2017
''''''''''''''''''''''''''''''''

jdtdb = jdtdb1 + 5.0

x = (jdtdb - 2451545.0) / 36525.0

teval(ta, tb, nfit%, ry_poly(), x, rmoon(1))

teval(ta, tb, nfit%, rx_poly(), x, rmoon(2))

teval(ta, tb, nfit%, rz_poly(), x, rmoon(3))

' print position vector and magnitude

print " "
print "Moon position vector and magnitude"
print "----------------------------------"
print " "

print "r_x = ", str$(rmoon(1), 0, 8), " kilometers"
print " "
print "r_y = ", str$(rmoon(2), 0, 8), " kilometers"
print " "
print "r_z = ", str$(rmoon(3), 0, 8), " kilometers"
print " "
print "rmag = ", str$(vecmag(rmoon()), 0, 8), " kilometers"
print " "

end

'''''''''''''''''''''''''''''''
'''''''''''''''''''''''''''''''

sub teval(fa, fb, n%, c(), x, tf)

' evaluate chebyshev series subroutine

''''''''''''''''''''''''''''''''''''''

f1 = 0.0

f2 = 0.0

xx = (2.0 * x - fa - fb) / (fb - fa)

xx2 = 2.0 * xx

for i% = n% to 1 step -1

oldf1 = f1

f1 = xx2 * f1 - f2 + c(i% + 1)

f2 = oldf1

next i%

tf = xx * f1 - f2 + 0.5 * c(1)

end sub

''''''''''''''''''
''''''''''''''''''

function vecmag(a())

' vector magnitude function

' input

' { a } = column vector ( 3 rows by 1 column )

' output

' vecmag = scalar magnitude of vector { a }

vecmag = sqr(a(1) * a(1) + a(2) * a(2) + a(3) * a(3))

end function

The zipped archive is located here 2017-04-21_180252_chebyshev.zip

Grogster

Admin Group

Joined: 31/12/2012
Location: New Zealand
Posts: 9748

Posted: 04:24pm 21 Apr 2017

Copy link to clipboard

Print this post

I have seen these posts, and am totally lost.
I am in awe of your mathematical abilities - the forum has apparently gained a maths professor.

I know who to ask next time I have a tricky mathematical problem to solve!

Smoke makes things work. When the smoke gets out, it stops!

cdeagle
Senior Member

Joined: 22/06/2014
Location: United States
Posts: 266

Posted: 08:16am 05 Jul 2017

Copy link to clipboard

Print this post

I've recently updated two of the MMBASIC-MMX programs for the Moon mentioned in the first post about this subject.

The create program now writes the ephemeris information to a CSV file.

The program that evaluates the ephemeris reads and uses the data file.

I'm happy to say that these two programs run on the MMX, Raspberry Pi and DOS versions of MMBASIC without changing any of the code. A rare happening during software development indeed.

Here's a link to the forum-hosted zipped file with the software.

2017-07-05_181144_chebyshev_moon.zip

Tinine
Guru

Joined: 30/03/2016
Location: United Kingdom
Posts: 1646

Posted: 10:33pm 13 Aug 2017

Copy link to clipboard

Print this post

Grogster said I have seen these posts, and am totally lost.
I am in awe of your mathematical abilities - the forum has apparently gained a maths professor.

I know who to ask next time I have a tricky mathematical problem to solve!

Print this page

To reply to this topic, you need to log in.