Correction: Thank's to Yu Cheng Lai for pointing out an error in the position for Pluto. The Pluto mean elements and the downloadable spreadsheet have now been corrected - 27th Feb 1999
This page provides details of a range of user functions for astronomy for Microsoft Excel version 5 or later. These functions may also be of interest to land surveyors and construction technicians, as well as teachers of astronomy and amateur astronomers. The accuracy of the algorithms used is not high enough for any navigational use, and there may well be errors in the code, although I have checked the results with other packages.
The algorithms shown here should be easily
translated into most other 'macro languages', including
AUTOLisp
, and SketchBasic
which presents the
possibility of automatic planisphere or diagramatic calendar
construction.
The modern business spreadsheet provides a powerful macro language, as well as the calculational facilities on the spreadsheet grid itself. All spreadsheets use radian measure for angles, and so the formulas used to calculate the positions of planets and other objects become rather long. If you want to make a table of positions, or an ephemeris for some object, you can end up with a large and complex spreadsheet! Worse, when you want to write a spreadsheet to do something different, you have to copy bits of the original, and change references.
You can write 'macros' to take a date, and produce a worksheet showing (say) the position, aspect phase and paralactic angle of the Moon. Macros do not usually 'hotlink' to the data, so for a new date you have to run the macro again.
The solution I adopt here is to use the facility in Excel for defining User Functions. User functions have no side effects, you cannot write a function which changes anything other than the cell which contains the function. This means that I have to return coordinates one by one by using an 'index' variable. As an example, to find the (geocentric) Right Ascension in hours of the Moon on 1998 September 27th at 15:50 UT would need two functions;
=moon(day2000(1998,9,27,15,50,0), 3)/15in general, the function
=moon(instant, index)returns the spherical coordinates of the moon,
index = 1 returns the distance in Earth radii index = 2 returns the declination in degrees index = 3 returns the right ascension in degreesThe spherical coordinates of any body are returned as;
index equatorial ecliptic (helio or geo centric) 1 distance (a.u.) distance (a.u.) 2 declination latitude 3 right ascension longitude
Note that the planet positions are referred to the J2000.0 equinox and ecliptic, you will find differences if you compare these positions with the apparent positions (referred to equinox and ecliptic of date). These differences amount to degrees for dates 50 years away from J2000.0. The Sun and Moon positions are referred to mean ecliptic and equinox of date.
A full tutorial on using these functions will appear when I have added support for rigourous precession, some ready made 'macros' for common tasks, and a range of precision choices. I also need to improve the error trapping!
Below you will find the Visual Basic for Applications source code for the set of user defined functions. Before the definition of each function, there is a brief comment describing the function.
To add these user defined functions to your spreadsheet, you should;
=day2000(1998,9,27,0,0,0)
into a cell, then press enter.-461.5
, i.e. the number of days before J2000.0
You may want to download the example spreadsheet (see below), which has all the functions plus some basic 'help' in the Function Wizard, and some example spreadsheets showing how the functions are used.
'************************************************************************ ' A series of astronomical functions which may be ' useful. These are all 'user defined functions'. ' This means that you can paste them into spreadsheets ' just like the normal functions - see Insert|function, ' you can even use the function wizard. ' The disadvantage ' of Excel's 'user defined functions' is that they ' can only return a single value, and the function cannot alter ' the properties of the worksheet. Arguments you pass to ' the VBA functions you define are passed 'by value'. ' However, VBA defaults to 'passing arguments by reference' ' when a function is called from another VBA function! This ' can lead to a function giving a different answer when ' called in the VBA module compared with when called in the ' spreadsheet. Use the ByVal keyword to tag arguments you ' change later in functions. See smoon() for an example. ' define some numerical constants - these are not ' accessible in the spreadsheet. Public Const pi As Double = 3.14159265358979 Public Const tpi As Double = 6.28318530717958 Public Const degs As Double = 57.2957795130823 Public Const rads As Double = 1.74532925199433E-02 ' The trig formulas working in degrees. This just ' makes the spreadsheet formulas a bit easier to ' read. DegAtan2() has had the arguments swapped ' from the Excel order, so the order matches most ' textbooks Function DegSin(x As Double) As Double DegSin = Sin(rads * x) End Function Function DegCos(x As Double) As Double DegCos = Cos(rads * x) End Function Function DegTan(x As Double) As Double DegTan = Tan(rads * x) End Function Function DegArcsin(x As Double) As Double DegArcsin = degs * Application.Asin(x) End Function Function DegArccos(x As Double) As Double DegArccos = degs * Application.Acos(x) End Function Function DegArctan(x As Double) As Double DegArctan = degs * Atn(x) End Function Function DegAtan2(y As Double, x As Double) As Double ' this returns the angle in the range 0 to 360 ' instead of -180 to 180 - and swaps the arguments ' This format matches Meeus and Duffett-Smith DegAtan2 = degs * Application.Atan2(x, y) If DegAtan2 < 0 Then DegAtan2 = DegAtan2 + 360 End Function Private Function range2pi(x) ' ' returns an angle x in the range 0 to two pi rads ' This function is not available in the spreadsheet ' range2pi = x - tpi * Int(x / tpi) End Function Private Function range360(x) ' ' returns an angle x in the range 0 to 360 ' used to prevent the huge values of degrees ' that you get from mean longitude formulas ' ' this function is private to this module, ' you won't find it in the Function Wizard, ' and you can't use it on a spreadsheet. ' If you want it on the spreadsheet, just remove ' the 'private' keyword above. ' range360 = x - 360 * Int(x / 360) End Function Function degdecimal(d, m, s) ' converts from dms format to ddd format degdecimal = d + m / 60 + s / 3600 End Function ' ' calander functions. jday and jcentury work on the Julian day numbers. ' day2000 and century2000 work on the days to J2000 to reduce the ' number of significant figures needed ' Function jday(year As Integer, month As Integer, day As Integer, hour As Integer, _ min As Integer, sec As Double, Optional greg) As Double ' returns julian day number given date in gregorian calender (greg=1) ' or julian calendar (greg = 0). If greg ommited, then Gregorian is assumed. Dim a As Double Dim b As Integer a = 10000# * year + 100# * month + day If (a < -47120101) Then MsgBox "Warning: date too early for algorithm" If (IsMissing(greg)) Then greg = 1 If (month <= 2) Then month = month + 12 year = year - 1 End If If (greg = 0) Then b = -2 + Fix((year + 4716) / 4) - 1179 Else b = Fix(year / 400) - Fix(year / 100) + Fix(year / 4) End If a = 365# * year + 1720996.5 jday = a + b + Fix(30.6001 * (month + 1)) + day + (hour + min / 60 + sec / 3600) / 24 End Function Function jcentury(jd As Double) As Double ' finds how many julian centuries since J2000 given ' the julian day number. Not used below, I just add ' a line into the subroutines which then take days ' before J2000 as the time argument jcentury = (jd - 2451545) / 36525 End Function Function day2000(year As Integer, month As Integer, day As Integer, hour As Integer, _ min As Integer, sec As Double, Optional greg) As Double ' returns days before J2000.0 given date in gregorian calender (greg=1) ' or julian calendar (greg = 0). If you don't provide a value for greg, ' then assumed Gregorian calender Dim a As Double Dim b As Integer If (IsMissing(greg)) Then greg = 1 a = 10000# * year + 100# * month + day If (month <= 2) Then month = month + 12 year = year - 1 End If If (greg = 0) Then b = -2 + Fix((year + 4716) / 4) - 1179 Else b = Fix(year / 400) - Fix(year / 100) + Fix(year / 4) End If a = 365# * year - 730548.5 day2000 = a + b + Fix(30.6001 * (month + 1)) + day + (hour + min / 60 + sec / 3600) / 24 End Function Function century2000(day2000 As Double) As Double ' finds how many julian centuries since J2000 given ' the days before J2000 century2000 = day2000 / 36525 End Function ' ' Conversion to and from rectangular and polar coordinates. ' X,Y,Z form a left handed set of axes, and r is the radius vector ' of the point from the origin. Theta is the elevation angle of ' r with the XY plane, and phi is the angle anti-clockwise from the ' X axis and the projection of r in the X,Y plane. ' ' in astronomical coordinate systems, ' ' item equatorial ecliptic (helio or geo centric) ' z celestial pole ecliptic pole ' x,y equatorial plane ecliptic ' theta declination latitude ' phi right ascension longitude ' Function rectangular(r As Double, theta As Double, phi As Double, _ index As Integer) As Double ' takes spherical coordinates in degrees and returns the rectangular ' coordinate shown by index, 1 = x, 2 = y, 3 = z ' ' x = r.cos(theta).cos(phi) ' y = r.cos(theta).sin(phi) ' z = r.sin(theta) ' Dim r_cos_theta As Double r_cos_theta = r * DegCos(theta) Select Case index Case 1 rectangular = r_cos_theta * DegCos(phi) 'returns x coord Case 2 rectangular = r_cos_theta * DegSin(phi) 'returns y coord Case 3 rectangular = r * DegSin(theta) 'returns z coord End Select End Function Function rlength(x As Double, y As Double, z As Double) As Double ' returns radius vector given the rectangular coords rlength = Sqr(x * x + y * y + z * z) End Function Function spherical(x As Double, y As Double, z As Double, index As Integer) As Double ' ' Takes the rectangular coordinates and returns the shperical ' coordinate selected by index - 1 = r, 2 = theta, 3 = phi ' ' r = sqrt(x*x + y*y + z*z) ' tan(phi) = y/x - use atan2 to get in correct quadrant ' tan(theta) = z/sqrt(x*x + y*y) - likewise ' Dim rho As Double rho = x * x + y * y Select Case index Case 1 spherical = Sqr(rho + z * z) 'returns r Case 2 rho = Sqr(rho) spherical = DegArctan(z / rho) 'returns theta Case 3 rho = Sqr(rho) spherical = DegAtan2(y, x) 'returns phi End Select End Function ' ' returns the obliquity of the ecliptic in degrees given the number ' of julian centuries from J2000 ' ' Most textbooks will give the IAU formula for the obliquity of ' the ecliptic below; ' ' obliquity = 23.43929111 - 46.8150"t - 0.00059"t^2 + 0.001813*t^3 ' ' as explained in Meeus or Numerical Recipes, it is more efficient and ' accurate to use the nested brackets shown in the function. If you ' multiply the brackets out, they come to the same. ' Function obliquity(d As Double) As Double Dim t As Double t = d / 36525 'julian centuries since J2000.0 obliquity = 23.43929111 - (46.815 + (0.00059 - 0.001813 * t) * t) * t / 3600# End Function ' ' functions for converting between equatorial and ecliptic ' geocentric coordinates, both polar and rectangular coords ' ' ' Converts geocentric ecliptic coordinates into geocentric equatorial ' coordinates. Expects rectangular coordinates. ' Function requatorial(x As Double, y As Double, z As Double, d As Double, _ index As Integer) As Double Dim obl As Double obl = obliquity(d) Select Case index Case 1 requatorial = x Case 2 requatorial = y * DegCos(obl) - z * DegSin(obl) Case 3 requatorial = y * DegSin(obl) + z * DegCos(obl) End Select End Function ' ' converts geocentric equatorial coordinates into geocentric ecliptic ' coordinates. Expects rectangular coordinates. ' Function recliptic(x As Double, y As Double, z As Double, d As Double, _ index As Integer) As Double Dim obl As Double obl = obliquity(d) Select Case index Case 1 recliptic = x Case 2 recliptic = y * DegCos(obl) + z * DegSin(obl) Case 3 recliptic = -y * DegSin(obl) + z * DegCos(obl) End Select End Function ' ' Converts geocentric ecliptic coordinates into geocentric equatorial ' coordinates. Expects spherical coordinates. ' Function sequatorial(r As Double, theta As Double, phi As Double, d As Double, _ index As Integer) As Double Dim x As Double, y As Double, z As Double x = rectangular(r, theta, phi, 1) y = rectangular(r, theta, phi, 2) z = rectangular(r, theta, phi, 3) sequatorial = spherical(requatorial(x, y, z, d, 1), requatorial(x, y, z, d, 2), _ requatorial(x, y, z, d, 3), index) End Function ' Converts geocentric equatorial coordinates into geocentric ecliptic ' coordinates. Expects spherical coordinates. ' Function secliptic(r As Double, theta As Double, phi As Double, d As Double, _ index As Integer) As Double Dim x As Double, y As Double, z As Double x = rectangular(r, theta, phi, 1) y = rectangular(r, theta, phi, 2) z = rectangular(r, theta, phi, 3) secliptic = spherical(recliptic(x, y, z, d, 1), recliptic(x, y, z, d, 2), _ recliptic(x, y, z, d, 3), index) End Function ' ' precession (approximate formula) from Meeus Algorithms p124 ' d1 is the epoch to precess from, d2 is the epoch to precess ' to, and index selects ra or dec. The function takes optional ' arguments dra and ddec to represent the proper motion of a ' star in seconds of arc per year. ' ' ra and dec must BOTH be in decimal degrees. This formula is ' different to the one elsewhere on the Web site! ' Function precess(d1 As Double, d2 As Double, dec As Double, _ ra As Double, index As Integer, _ Optional ddec, Optional dra) As Double Dim m As Double, n As Double, t As Double If (IsMissing(dra)) Then dra = 0 If (IsMissing(ddec)) Then ddec = 0 t = d1 / 36525 'years since J2000 m = 0.01281233333333 + 0.00000775 * t n = 0.005567527777778 - 2.361111111111E-06 * t t = (d2 - d1) / 365.25 'difference in julian _years_, not centuries! Select Case index Case 1 'dec precess = dec + (n * DegCos(ra) + ddec / 3600) * t Case 2 'ra precess = ra + (m + n * DegSin(ra) * DegTan(dec) + dra / 3600) * t End Select End Function ' ' The function below returns the geocentric ecliptic coordinates of the sun ' to an accuracy corresponding to about 0.01 degree over ' 100 years either side of J2000. Coordinates returned ' in spherical form. From page C24 of the 1996 Astronomical ' Alamanac. Comparing accuracy with Planeph, a DOS ephemeris ' by Chapront, we get; ' ' Sun error ' RA sec DEC arcsec ' Max within 3 year 0.6 8.9 ' Min within 3 year -2.1 -8.2 ' Max within 10 year 0.6 10.9 ' Min within 10 year -2.6 -12.5 ' Max within 50 year 1.0 16.8 ' Min within 50 year -2.9 -16.1 ' ' Error = C24 low precision method - Planeph ' ' Note: Planeph was set to give output referred to mean ' ecliptic and equinox of date. ' ' The accuracy of this routine is good enough for sunrise ' and shadow direction calculations, and for referring ' low precision planetary and comet positions to the Earth, ' but is no good for accurate coordinate conversion or ' for eclipse or occultation use. ' ' Coordinates are referred to the ecliptic and mean equinox of date ' Function ssun(d As Double, index As Integer) As Double Dim g As Double Dim l As Double g = range360(357.528 + 0.9856003 * d) l = range360(280.461 + 0.9856474 * d) Select Case index Case 1 ' radius vector of Sun ssun = 1.00014 - 0.01671 * DegCos(g) - 0.00014 * DegCos(2 * g) Case 2 ssun = 0 'ecliptic latitude of Sun is zero to very good accuracy Case 3 'longitude of Sun ssun = range360(l + 1.915 * DegSin(g) + 0.02 * DegSin(2 * g)) End Select End Function ' returns the geocentric ecliptic coordinates of the sun ' to an accuracy corresponding to about 0.01 degree over ' 100 years either side of J2000. Assumes ssun() exists. ' ' rectangular form is easier for converting positions from helio- ' centric to geocentric, but beware low accuracy (roughly 0.01 degree) ' of values ' Function rsun(d As Double, index As Integer) As Double Dim x As Double Dim y As Double Dim z As Double rsun = rectangular(ssun(d, 1), ssun(d, 2), ssun(d, 3), index) End Function ' ' sun() returns the geocentric ra and dec and radius vector ' of the moon - calls smoon three times, and sequatorial ' three times - sequatorial calls rectangular three times ' each! ' Function sun(d As Double, index As Integer) As Double sun = sequatorial(ssun(d, 1), ssun(d, 2), ssun(d, 3), d, index) End Function ' ' The function below implements Paul Schlyter's simplification ' of van Flandern and Pulkkinen's method for finding the geocentric ' ecliptic positions of the Moon to an accuracy of about 1 to 4 arcmin. ' ' I can probably reduce the number of variables, and there must ' be a quicker way of declaring variables! ' ' The VBA trig functions have been used throughout for speed, ' note how the atan function returns values in domain -pi to pi ' Function smoon(ByVal d As Double, index As Integer) As Double Dim Nm As Double, im As Double, wm As Double, am As Double, ecm As Double, _ Mm As Double, em As Double, Ms As Double, ws As Double, xv As Double, _ yv As Double, vm As Double, rm As Double, x As Double, y As Double, _ z As Double, lon As Double, lat As Double, ls As Double, lm As Double, _ dm As Double, F As Double, dlong As Double, dlat As Double ' Paul's routine uses a slightly different definition of ' the day number - I adjust for it below. Remember that VBA ' defaults to 'pass by reference' so this change in d ' will be visible to other functions unless you set d to 'ByVal' ' to force it to be passed by value! d = d + 1.5 ' moon elements Nm = range360(125.1228 - 0.0529538083 * d) * rads im = 5.1454 * rads wm = range360(318.0634 + 0.1643573223 * d) * rads am = 60.2666 '(Earth radii) ecm = 0.0549 Mm = range360(115.3654 + 13.0649929509 * d) * rads ' position of Moon em = Mm + ecm * Sin(Mm) * (1# + ecm * Cos(Mm)) xv = am * (Cos(em) - ecm) yv = am * (Sqr(1# - ecm * ecm) * Sin(em)) vm = Application.Atan2(xv, yv) ' If vm < 0 Then vm = tpi + vm rm = Sqr(xv * xv + yv * yv) x = rm * (Cos(Nm) * Cos(vm + wm) - Sin(Nm) * Sin(vm + wm) * Cos(im)) y = rm * (Sin(Nm) * Cos(vm + wm) + Cos(Nm) * Sin(vm + wm) * Cos(im)) z = rm * (Sin(vm + wm) * Sin(im)) ' moons geocentric long and lat lon = Application.Atan2(x, y) If lon < 0 Then lon = tpi + lon lat = Atn(z / Sqr(x * x + y * y)) ' mean longitude of sun ws = range360(282.9404 + 0.0000470935 * d) * rads Ms = range360(356.047 + 0.9856002585 * d) * rads ' perturbations ' first calculate arguments below, 'Ms, Mm Mean Anomaly of the Sun and the Moon 'Nm Longitude of the Moon's node 'ws, wm Argument of perihelion for the Sun and the Moon ls = Ms + ws 'Mean Longitude of the Sun (Ns=0) lm = Mm + wm + Nm 'Mean longitude of the Moon dm = lm - ls 'Mean elongation of the Moon F = lm - Nm 'Argument of latitude for the Moon ' then add the following terms to the longitude ' note amplitudes are in degrees, convert at end Select Case index Case 1 ' distance terms earth radii rm = rm - 0.58 * Cos(Mm - 2 * dm) rm = rm - 0.46 * Cos(2 * dm) smoon = rm Case 2 ' latitude terms dlat = -0.173 * Sin(F - 2 * dm) dlat = dlat - 0.055 * Sin(Mm - F - 2 * dm) dlat = dlat - 0.046 * Sin(Mm + F - 2 * dm) dlat = dlat + 0.033 * Sin(F + 2 * dm) dlat = dlat + 0.017 * Sin(2 * Mm + F) smoon = lat * degs + dlat Case 3 ' longitude terms dlon = -1.274 * Sin(Mm - 2 * dm) '(the Evection) dlon = dlon + 0.658 * Sin(2 * dm) '(the Variation) dlon = dlon - 0.186 * Sin(Ms) '(the Yearly Equation) dlon = dlon - 0.059 * Sin(2 * Mm - 2 * dm) dlon = dlon - 0.057 * Sin(Mm - 2 * dm + Ms) dlon = dlon + 0.053 * Sin(Mm + 2 * dm) dlon = dlon + 0.046 * Sin(2 * dm - Ms) dlon = dlon + 0.041 * Sin(Mm - Ms) dlon = dlon - 0.035 * Sin(dm) '(the Parallactic Equation) dlon = dlon - 0.031 * Sin(Mm + Ms) dlon = dlon - 0.015 * Sin(2 * F - 2 * dm) dlon = dlon + 0.011 * Sin(Mm - 4 * dm) smoon = lon * degs + dlon End Select End Function ' ' rmoon uses smoon to return the geocentric ecliptic rectangular coordinates ' of the moon to lowish accuracy. ' ' Function rmoon(d As Double, index As Integer) As Double rmoon = rectangular(smoon(d, 1), smoon(d, 2), smoon(d, 3), index) End Function ' ' moon() returns the geocentric ra and dec and radius vector ' of the moon - calls smoon three times, and sequatorial ' three times - sequatorial calls rectangular three times ' each! ' Function moon(d As Double, index As Integer) As Double Dim nigel As Double If index = 4 Then nigel = smoon(d, 2) moon = d Else moon = sequatorial(smoon(d, 1), smoon(d, 2), smoon(d, 3), d, index) End If End Function ' ' Solutions of the kepler equation using a Newton's method approach. ' See Meeus or Duffett-Smith (calculator) ' Function kepler(m As Double, ecc As Double, Optional eps) ' solves the equation e - ecc*sin(e) = m for e given an m ' returns the value of the 'true anomaly' in rads ' m - the 'mean anomaly' in rads ' ecc - the eccentricity of the orbit ' eps - the precision parameter - solution will be ' within 10^-eps of the true value. ' don't set eps above 14, as convergence ' can't be guaranteed. If not specified, then ' taken as 10^-8 or 10 nano radians! ' Dim delta As Double, e As Double, v As Double e = m 'first guess delta = 0.05 'set delta equal to a dummy value If (IsMissing(eps)) Then eps = 8 'if no eps then assume 10^-8 Do While Abs(delta) >= 10 ^ -eps 'converged? delta = e - ecc * Sin(e) - m 'new error e = e - delta / (1 - ecc * Cos(e)) 'corrected guess Loop v = 2 * Atn(((1 + ecc) / (1 - ecc)) ^ 0.5 * Tan(0.5 * e)) If v < 0 Then v = v + tpi kepler = v End Function ' ' The functions below return the heliocentric ecliptic coordinates ' of the planets to an accuracy of a few minutes of arc. The coordinates ' are referred to the equinox of J2000.0 ' ' The functions use a simple Kepler ellipse, but with ' mean elements which change slightly with the time since ' J2000. See 'Explanatory supplement to the Astronomical ' Almanac' 1992, page 316, table 5.8.1. Worst case errors ' over the period 1800 - 2050 AD in arcsec are below; ' ' Ra Dec ' Mercury 20" 5" ' Venus 20" 5" ' Earth 20" 5" ' Mars 25" 30" ' Jupiter 300" 100" ' Saturn 600" 200" ' Uranus 60" 25" ' Neptune 40" 20" ' Pluto 40" 10" ' ' ' The rplanet() function returns the ecliptic heliocentric coordinates ' of each of the major planets. You select the planet you want using ' pnumber, and the coordinate you want with index as usual. ' Function rplanet(d As Double, pnumber As Integer, index As Integer) As Double Dim x As Double, y As Double, z As Double, v As Double, m As Double, _ i As Double, o As Double, p As Double, a As Double, e As Double, _ l As Double, r As Double ' ' get elements of the planet ' element i, o, p, a, e, l, d, pnumber ' ' position of planet in its orbit ' m = range2pi(l - p) v = kepler(m, e, 8) r = a * (1 - e * e) / (1 + e * Cos(v)) ' ' heliocentric rectangular coordinates of planet ' Select Case index Case 1 'x coordinate rplanet = r * (Cos(o) * Cos(v + p - o) - Sin(o) * Sin(v + p - o) * Cos(i)) Case 2 'y coordinate rplanet = r * (Sin(o) * Cos(v + p - o) + Cos(o) * Sin(v + p - o) * Cos(i)) Case 3 'z coordinate rplanet = r * (Sin(v + p - o) * Sin(i)) End Select End Function ' ' ' The planet() function returns the equatorial geocentric coordinates ' of each of the major planets. You select the planet you want using ' pnumber, and the coordinate you want with index as usual. Code is ' duplicated from rplanet() to reduce the number of calls to kepler() ' ' Function planet(d As Double, pnumber As Integer, index As Integer) As Double Dim x As Double, y As Double, z As Double, v As Double, m As Double, _ i As Double, o As Double, p As Double, a As Double, e As Double, _ l As Double, r As Double, xe As Double, ye As Double, ze As Double, _ s1 As Double, si As Double, so As Double, c1 As Double, ci As Double, _ co As Double ' ' elements of planet - select from the values ' element i, o, p, a, e, l, d, pnumber ' ' position of planet in its orbit ' m = range2pi(l - p) v = kepler(m, e, 8) r = a * (1 - e * e) / (1 + e * Cos(v)) ' ' heliocentric rectangular coordinates of planet ' s1 = Sin(v + p - o) si = Sin(i) so = Sin(o) c1 = Cos(v + p - o) ci = Cos(i) co = Cos(o) x = r * (co * c1 - so * s1 * ci) y = r * (so * c1 + co * s1 * ci) z = r * (s1 * si) ' ' elements of earth (reusing variables) ' element i, o, p, a, e, l, d, 3 ' ' position of earth in its orbit ' m = range2pi(l - p) v = kepler(m, e, 8) r = a * (1 - e * e) / (1 + e * Cos(v)) ' ' heliocentric rectangular coordinates of earth ' s1 = Sin(v + p - o) si = Sin(i) so = Sin(o) c1 = Cos(v + p - o) ci = Cos(i) co = Cos(o) xe = r * (co * c1 - so * s1 * ci) ye = r * (so * c1 + co * s1 * ci) ze = r * (s1 * si) ' ' convert to geocentric rectangular coordinates ' x = x - xe y = y - ye ' z = z ' ' rotate around x axis from ecliptic to equatorial coords ' ecl = 23.429292 * rads 'value for J2000.0 frame xe = x ye = y * Cos(ecl) - z * Sin(ecl) ze = y * Sin(ecl) + z * Cos(ecl) ' ' find the RA and DEC from the rectangular equatorial coords ' Select Case index Case 3 ' RA in degrees planet = Application.Atan2(xe, ye) * degs If planet < 0 Then planet = 360 + planet Case 2 ' DEC in degrees planet = Atn(ze / Sqr(xe * xe + ye * ye)) * degs Case 1 ' Radius vector in au planet = Sqr(xe * xe + ye * ye + ze * ze) End Select End Function ' ' The subroutine below replaces the values of i,o,p,a,e,L ' with the values for the planet selected by pnum. You could ' always add planet like objects, but watch the value of ' the inclination i. The method used in planet is only ' good for orbits 'near' the ecliptic. ' Sub element(i As Double, o As Double, p As Double, _ a As Double, e As Double, l As Double, _ ByVal d As Double, ByVal pnum As Integer) Select Case pnum Case 1 'mercury i = (7.00487 - 0.000000178797 * d) * rads o = (48.33167 - 0.0000033942 * d) * rads p = (77.45645 + 0.00000436208 * d) * rads a = 0.38709893 + 1.80698E-11 * d e = 0.20563069 + 0.000000000691855 * d l = range2pi(rads * (252.25084 + 4.092338796 * d)) Case 2 'venus i = (3.39471 - 0.0000000217507 * d) * rads o = (76.68069 - 0.0000075815 * d) * rads p = (131.53298 - 0.000000827439 * d) * rads a = 0.72333199 + 2.51882E-11 * d e = 0.00677323 - 0.00000000135195 * d l = range2pi(rads * (181.97973 + 1.602130474 * d)) Case 3 'earth i = (0.00005 - 0.000000356985 * d) * rads o = (-11.26064 - 0.00013863 * d) * rads p = (102.94719 + 0.00000911309 * d) * rads a = 1.00000011 - 1.36893E-12 * d e = 0.01671022 - 0.00000000104148 * d l = range2pi(rads * (100.46435 + 0.985609101 * d)) Case 4 'mars i = (1.85061 - 0.000000193703 * d) * rads o = (49.57854 - 0.0000077587 * d) * rads p = (336.04084 + 0.00001187 * d) * rads a = 1.52366231 - 0.000000001977 * d e = 0.09341233 - 0.00000000325859 * d l = range2pi(rads * (355.45332 + 0.524033035 * d)) Case 5 'jupiter i = (1.3053 - 0.0000000315613 * d) * rads o = (100.55615 + 0.00000925675 * d) * rads p = (14.75385 + 0.00000638779 * d) * rads a = 5.20336301 + 0.0000000166289 * d e = 0.04839266 - 0.00000000352635 * d l = range2pi(rads * (34.40438 + 0.083086762 * d)) Case 6 'saturn i = (2.48446 + 0.0000000464674 * d) * rads o = (113.71504 - 0.0000121 * d) * rads p = (92.43194 - 0.0000148216 * d) * rads a = 9.53707032 - 0.0000000825544 * d e = 0.0541506 - 0.0000000100649 * d l = range2pi(rads * (49.94432 + 0.033470629 * d)) Case 7 'uranus i = (0.76986 - 0.0000000158947 * d) * rads o = (74.22988 + 0.0000127873 * d) * rads p = (170.96424 + 0.0000099822 * d) * rads a = 19.19126393 + 0.0000000416222 * d e = 0.04716771 - 0.00000000524298 * d l = range2pi(rads * (313.23218 + 0.011731294 * d)) Case 8 'neptune i = (1.76917 - 0.0000000276827 * d) * rads o = (131.72169 - 0.0000011503 * d) * rads p = (44.97135 - 0.00000642201 * d) * rads a = 30.06896348 - 0.0000000342768 * d e = 0.00858587 + 0.000000000688296 * d l = range2pi(rads * (304.88003 + 0.0059810572 * d)) Case 9 'pluto i = (17.14175 + 0.0000000841889 * d) * rads o = (110.30347 - 0.0000002839 * d) * rads p = (224.06676 - 0.00000100578 * d) * rads a = 39.48168677 - 0.0000000210574 * d e = 0.24880766 + 0.00000000177002 * d l = range2pi(rads * (238.92881 + 3.97557152635181E-03 * d)) End Select End Sub '****************************************************************
You can download an Excel 95 spreadsheet with a VBA module containing all the functions
above, and some example applications, as
astrovba.zip (32 Kb)
.
Below is a list of all the Excel spreadsheets available from this Web site;
oscxls.zip
, mercury worksheet, cell G14
will convince you that user functions can make spreadsheets neater to
write and easier to follow!
Meeus, Jean
Astronomical Algorithms
Willmann-Bell
1st English edition, 1991
ISBN 0-943396-35-2
Duffett-Smith, Peter
Practical Astronomy with your calculator
Cambridge University Press
3rd edition 1988
ISBN 0-521-35699-7
Montenbruck, Pfleger
Astronomy on the personal computer
Springer
3rd edition 1988
ISBN 3-540-63521-1
[Root]
Last Modified 27th February 1999
Keith Burnett