Planning a Java astronomy library for programmers
Kerry Shetline
on this problem goes over.
I'm implementing and extending the methods found in "Low-precision
Formulae for Planetary Positions" by T.C. Van Flandern and K. F.
Pulkkinen (U.S. Naval Observatory, 1979) in Java. I've seen other works
based on this paper, but none that implement the series to their full
precision -- within one arc-second for +/- 300 years. It's a lot of
tedious typing to get all of those series terms in that paper entered.
So far I've got the data entered for the sun and moon, as well as
Mercury, Venus, and Mars. I've been spending more time coding in Java
than data entry, but I should soon have all of the planets out to Pluto
entered.
My goal is to create a flexible, extensible set of classes and methods
that are easily incorporated into any Java program. Beyond computing the
planetary elements themselves as given in the Van Flandern/Pulkkinen
paper, I've developed methods for computing moon phase events, solstice
and equinox events, rise and set times, and a number of utility
functions. In the future, I'd like to add elongation and opposition
events, Jovian moons, and any other good suggestions that come my way,
that I feel capable of implementing. Once the code is further developed,
I intend to make it all publicly available.
I'd even like to implement the full precision of the VSOP87 series some
day, but I'd need to find better resources than I've found so far on how
to compute and use these series. In keeping with that level of accuracy,
I'd need to learn more about how to compute things like stellar
aberration and the difference between Ephemeris Time and Universal Time.
Any help and suggestions here would be greatly appreciated.
Other resources used to create this code:
http://193.12.249.96/pausch/comp/ppcomp.html
http://www.skypub.com/
The API I have created so far:
Class: VPSolarSystem
Package: org.shetline.astronomy
* public double computeElement(int whichPlanet, int whichElement, double
tdays)
This returns the value of one of a series (PLON, BETA, RP, V, U, or W)
for a given planet (MOON, MERCURY, VENUS, SUN, MARS, JUPITER, SATURN,
URANUS, NEPTUNE, or PLUTO) as described in the Van Flandern/Pulkkinen
paper. The time tdays is given in days from January 1.5 (noon) 2000 ET.
At the given precision, ET, UT and GMT time scales are considered
interchangable. Many of the values returned are accurate enough,
however, that it might be beneficial to make the distinction. Conversion
methods for these time scales, anyone?
* public double latitude(int whichPlanet, double tdays)
Returns latitude of planet in degrees, [-90, 90].
* public double longitude(int whichPlanet, double tdays)
Returns longitude of planet in degrees, [0, 360).
* public double rightAscension(int whichPlanet, double tdays)
Returns right ascension of planet as an hour angle, [0, 24).
* public double declination(int whichPlanet, double tdays)
Returns declination of planet in degrees, [-90, 90].
* public double meanLongitudeOfSun(double tdays)
Returns the mean longitude of the sun in degrees, [0, 360).
* public double greenwichMeanSiderealTime(double tdays)
Returns sidereal time as an hour angle, [0, 24).
* double localSiderealTime(double tdays)
Returns sidereal time as an hour angle, [0, 24).
* public double[] altitudeAndAzimuth(int whichPlanet, double tdays,
double latitude, double longitude)
Given latitude and longitude in degrees, returns: result[0], altitude in
degrees, [-90, 90]; result[1], azimuth in degrees, [0, 360). NOTE: The
position of the Moon *is* corrected for topocentric displacement.
* public double momentaryLunarPhase(double tdays)
Returns continuously-variable value in degrees, [0, 360). Key values: 0
- new, 90 - first quarter, 180 - full, 270 - last quarter.
* public double percentLunarIllumination(double tdays)
Returns percentage of lunar illumination, [0, 100].
* public int enumeratedLunarPhase(int year, int month, int day, int[]
hourAndMinute)
Returns IN_BETWEEN_PHASES (-1) if no special phase occurs on given date
(GMT). Otherwise, returns NEW_MOON (0), FIRST_QUARTER (1), FULL_MOON
(2), or LAST_QUARTER (3), with the hour and minute of the event saved in
hourAndMinute[0] and hourAndMinute[1], respectively.
I've found these times to be typically within one minute of almanac
times for recent years, within 1-3 minutes going back to 1900.
(Routines for getting and displaying moon phase info by whole years
omitted.)
* public int equinoxOrSoltice(int year, int month, int day, int[]
hourAndMinute)
Returns NOT_EQUINOX_OR_SOLSTICE (-1) if no equinox or solstic occurs on
given date (GMT). Otherwise, returns VERNAL_EQUINOX (0), SUMMER_SOLSTICE
(1), AUTUMNAL_EQUINOX (2), or WINTER_SOLSTICE (3), with the hour and
minute of the event saved in hourAndMinute[0] and hourAndMinute[1],
respectively.
I've found these times to be typically within 1-3 minutes of almanac
times that I've checked.
(Routines for getting and displaying yearly equinox/solstice info
omitted.)
* public Vector getRiseAndSetTimes(int whichPlanet, int year, int month,
int day,
int timeZoneAdjInMinutes, double observerLatitude, double
observerLongitude)
Returns a vector containing 0 or more rise or set events. Each event is
and object of the class AstroEvent, which gives the time of the event,
and an integer value which tells you the type of event (RISE_EVENT,
SET_EVENT, VISIBLE_ALL_DAY, or UNSEEN_ALL_DAY). This method
automatically adjusts for typical atmospheric refraction, and another
quarter of a degree for the angular radius of the sun or moon. It is
possible to have a rise without a set, a set without a rise, and at
higher latitudes two rises or two sets can occur within one day.
* public Vector getRiseAndSetTimes(int whichPlanet, int year, int month,
int day,
int timeZoneAdjInMinutes, double observerLatitude, double
observerLongitude,
double angularSizeAndRefractionAdj)
Same as above, but allows you to explicitly specify the angle below the
mathematical horizon for the center of the planet. You can use 6, 12, or
18 degrees, for instance, to compute civil, nautical, and astronomical
twilight.
Class: AstronomyUtil
Package: org.shetline.astronomy
* public static int julianDate(int year, int month, int day)
* public static double timeInDays(int year, int month, int day,
int hour, int minute, double second)
* public static int daysInMonth(int year, int month)
* public static double[] altitudeAndAzimuth(double rightAscension,
double declination,
double localSiderealTime, double latitude, double longitude)
* public static double refractedAltitude(double trueAltitude)
Other utility methods:
public static String formatAngle(double degrees)
public static String formatHourAngle(double hourAngle)
public static String formatShortDate(int year, int month, int day)
public static String formatShortDate(int year, int month, int day, char
separator)
public static String formatLongDate(int year, int month, int day)
public static String formatLongDate(int year, int month, int day, char
separator)
public static String formatTime(int hour, int minute)
public static String formatTime(int hour, int minute, int second)
-Kerry
Mike Witters
writing a 'lite' planetarium software app. Your classes will surely come in
handy.
Mike Witters
Software Engineer
The Build-It-Yourself Astronomer Website
http://www.siscom.net/~mikew/astro/diy.htm
Paul Schlyter
Kerry Shetline <shet...@cyberzone.net> wrote:
> I know this isn't a completely new idea, but I'd like to see how my take
> on this problem goes over.
>
> I'm implementing and extending the methods found in "Low-precision
> Formulae for Planetary Positions" by T.C. Van Flandern and K. F.
> Pulkkinen (U.S. Naval Observatory, 1979) in Java. I've seen other works
> based on this paper, but none that implement the series to their full
> precision -- within one arc-second for +/- 300 years.
These series, published back in 1979, do have a claimed accuracy to
only one arc minute (except Pluto, which is accurate to only some 15
arc min). In practice, the heliocentric series (except Pluto) do
have better accuracy - perhaps 10-15 arc seconds.
Tom van Flandern did write in that paper from 1979 that he planned to
develop longer series, which had an accuracy of 1 arc second or
better, over a larger time interval; however that project never got
completed. And now it's overdue, since we have VSOP87 and ELP-2000.
> It's a lot of
> tedious typing to get all of those series terms in that paper entered.
> So far I've got the data entered for the sun and moon, as well as
> Mercury, Venus, and Mars. I've been spending more time coding in Java
> than data entry, but I should soon have all of the planets out to Pluto
> entered.
Perhaps I could save you some time. Below are Van Flandern-Pulkkinen's
series, coded by me in C, back in 1987. I only coded the heliocentric
series, since the geocentric series were both longer and less accurate,
and one can compute the geocentric positions from the heliocentric
positions anyway. I ended up extensively using only the series for
Jupiter and Saturn; for the other celestial bodies I picked the series
by Jean Meeus in his "Astronomical formulae for Calculators" instead,
which I found to be slightly more accurate.
> I'd even like to implement the full precision of the VSOP87 series some
> day, but I'd need to find better resources than I've found so far on how
> to compute and use these series.
Then you may need to learn some FORTRAN, since there are implementations
in FORTRAN of these series, to their full precision, on the Net. All
the ASCII source code is some 35 MBytes in size....
> In keeping with that level of accuracy,
> I'd need to learn more about how to compute things like stellar
> aberration and the difference between Ephemeris Time and Universal Time.
> Any help and suggestions here would be greatly appreciated.
>
> Other resources used to create this code:
> http://193.12.249.96/pausch/comp/ppcomp.html
Nowadays instead use:
http://hotel04.ausys.se/pausch/comp/ppcomp.html
Below is my old code for these series. I could have zipped and
uu-coded it to make it smaller, but then I would have been posting a
binary, which isn't allowed here. Therefore, to obey the charter of
this newsgroup, I instead post it in ASCII -- yes, this will make
the post perhaps 3 times larger, but I will then obey the charter
of this newsgroup... :-)
Enjoy!
==================================cut===================================
/*
TRIG.H
Trigonometric and angle funtions, sign & int funct:s for Turbo C
Paul Schlyter, 1987-06-13
Macros DM and DMS added, 1989-05-08
Macros M and S changed to MIN and SEC
*/
#ifndef _TRIG_H
#define _TRIG_H
#define PI 3.1415926535897932384
#define TWOPI ( 2.0 * PI )
#define HALFPI ( PI / 2.0 )
#define RADEG ( 180.0 / PI )
#define DEGRAD ( PI / 180.0 )
#define INV360 ( 1.0 / 360.0 )
#define ONEHALF 0.5
#define DM(d,m) ( (d) + (m)/60.0 )
#define DMS(d,m,s) ( (d) + (m)/60.0 + (s)/3600.0 )
#define MIN(m) ( (m)/60.0 )
#define SEC(s) ( (s)/3600.0 )
#define S(s) ( (s)/3600.0 )
/* Sign function + a better Int function */
#define sign(sgn,magn) ( sgn==0 ? 0 : sgn>0 ? magn : -magn )
/* Trig macros for degree values */
#define sind(x) sin((x)*DEGRAD)
#define cosd(x) cos((x)*DEGRAD)
#define tand(x) tan((x)*DEGRAD)
#define atand(x) (RADEG*atan(x))
#define asind(x) (RADEG*asin(x))
#define acosd(x) (RADEG*acos(x))
#define atan2d(y,x) (RADEG*atan2(y,x))
/* VARV.C, normalize angle within one revoultion */
double varv180(double angle);
double varv360(double angle);
==================================cut===================================
/*
VARV.C
Reduce angles within one revolution
Paul Schlyter, 1987-06-13
*/
#include <math.h>
#include "trig.h"
#define INV360 ( 1.0 / 360.0 )
#define ONEHALF 0.5
double varv180 ( double x )
/*************************/
/* Reduce angle to within -180..+180 degrees */
{
return x - 360.0 * floor( x * INV360 + ONEHALF );
}
double varv360 ( double x )
/*************************/
/* Reduce angle to within 0..360 degrees */
{
return x - 360.0 * floor( x * INV360 );
}
==================================cut===================================
/*
SOLPOSVF.C
Suns position
Referens: van Flandern & Pulkkinen
Astrophys. J. Suppl., 41:391-411, nov 1979
Paul Schlyter, 1987-06-14
*/
#include <stddef.h>
#include <math.h>
#include "trig.h"
static double frac( double x )
/****************************/
{
return( x - floor(x) );
} /* frac */
void solpos_vf( double JD, double *slon, double *slat, double *sr )
/*****************************************************************/
{
double dnr, t, ls, gs, lm, dm, g2, g4, g5,
sgs, cgs, sg2, cg2, sg5, cg5, s5, c5, s2, c2, s2gs, c2gs;
dnr = JD - 2451545.0;
t = dnr/36525.0 + 1.0;
ls = 360.0 * frac( 0.779072 + 0.00273790931 * dnr );
gs = 360.0 * frac( 0.993126 + 0.00273777850 * dnr );
/* Moon */
lm = 360.0 * frac( 0.606434 + 0.03660110129 * dnr );
dm = lm - ls;
/* Venus, Mars, Jupiter: */
g2 = 360.0 * frac( 0.140023 + 0.00445036173 * dnr );
g4 = 360.0 * frac( 0.053856 + 0.00145561327 * dnr );
g5 = 360.0 * frac( 0.056531 + 0.00023080893 * dnr );
/* Suns position: */
sgs = sind(gs);
cgs = cosd(gs);
sg2 = sind(g2);
cg2 = cosd(g2);
sg5 = sind(g5);
cg5 = cosd(g5);
s5 = sgs*cg5-cgs*sg5;
c5 = cgs*cg5+sgs*sg5;
s2 = sgs*cg2-cgs*sg2;
c2 = cgs*cg2+sgs*sg2;
s2gs = 2*sgs*cgs;
c2gs = 2*cgs*cgs-1.0;
*slon = ls + (
(6910.-17.*t) * sgs
+ 72. * s2gs
- 7. * c5
+ 6. * sind(dm)
+ 6.4 * sind(4*gs-8*g4+3*g5+38.9)
- 5. * (2*c2*c2-1.0)
- 4. * s2
+ 3. * 2*s2*c2
- 3. * sg5
- 3. * 2*s5*c5
) / 3600.0;
*sr = 1.00014 - 0.01675 * cgs - 0.00014 * c2gs;
if ( slat != NULL )
*slat = 0.0;
} /* solpos_vf */
==================================cut===================================
/*
MONPOSVF.C
Moons position
Referens: van Flandern & Pulkkinen
Astrophys. J. Suppl., 41:391-411, nov 1979
Paul Schlyter, 1987-06-14
*/
#include <math.h>
#include "trig.h"
static double frac( double x )
/****************************/
{
return( x - floor(x) );
} /* frac */
void monpos_vf( double JD, double *lon, double *lat, double *par )
/****************************************************************/
{
double dnr, t, ls, gs, lm, gm, fm, nm, dm, venarg, r,
sgs, cgs, sgm, cgm, sfm, cfm, sdm, cdm, snm, cnm, s2gs, c2gs,
s2gm, c2gm, s3gm, c3gm, s4gm, /* c4gm, */
s2fm, c2fm, s3fm, c3fm,
s2dm, c2dm, s3dm, c3dm, s4dm, c4dm,
sgmpgs, cgmpgs, sgmmgs, cgmmgs, sgmp2d, cgmp2d, sgmm2g,
cgmm2g, s2mmgs, c2mmgs, s2fm2d, c2fm2d, sgmm2d, cgmm2d,
s2gm2d, c2gm2d, sfmp2d, cfmp2d, sfmm2d, cfmm2d,
s2gpgs, c2gpgs, s2gp2d, c2gp2d,
sfmp4d, cfmp4d, sfmm4d, cfmm4d;
dnr = JD - 2451545.0;
t = dnr / 36525.0 + 1.0;
/* Sun: */
ls = 360.0 * frac( 0.779072 + 0.00273790931 * dnr );
gs = 360.0 * frac( 0.993126 + 0.00273777850 * dnr );
/* Moon: */
lm = 360.0 * frac( 0.606434 + 0.03660110129 * dnr );
gm = 360.0 * frac( 0.374897 + 0.03629164709 * dnr );
fm = 360.0 * frac( 0.259091 + 0.03674819520 * dnr );
nm = lm - fm;
dm = lm - ls;
/* Planetary perturbation,
argument: gm + 16*ls - 18*l2, l2 = Venus mean long. */
venarg = 360.0 * frac( 0.741085 - 1.024001E-05 * dnr );
sgs = sind(gs);
cgs = cosd(gs);
sgm = sind(gm);
cgm = cosd(gm);
sfm = sind(fm);
cfm = cosd(fm);
sdm = sind(dm);
cdm = cosd(dm);
snm = sind(nm);
cnm = cosd(nm);
s2gs = 2*sgs*cgs;
c2gs = 2*cgs*cgs-1.0;
s2gm = 2*sgm*cgm;
c2gm = 2*cgm*cgm-1.0;
s3gm = s2gm*cgm+c2gm*sgm;
c3gm = c2gm*cgm-s2gm*sgm;
s4gm = 2*s2gm*c2gm;
/* c4gm = 2*c2gm*c2gm-1.0; */
s2fm = 2*sfm*cfm;
c2fm = 2*cfm*cfm-1.0;
s3fm = s2fm*cfm+c2fm*sfm;
c3fm = c2fm*cfm-s2fm*sfm;
s2dm = 2*sdm*cdm;
c2dm = 2*cdm*cdm-1.0;
s3dm = s2dm*cdm+c2dm*sdm;
c3dm = c2dm*cdm-s2dm*sdm;
s4dm = 2*s2dm*c2dm;
c4dm = 2*c2dm*c2dm-1.0;
sgmpgs = sgm*cgs+cgm*sgs;
cgmpgs = cgm*cgs-sgm*sgs;
sgmmgs = sgm*cgs-cgm*sgs;
cgmmgs = cgm*cgs+sgm*sgs;
sgmp2d = sgm*c2dm+cgm*s2dm;
cgmp2d = cgm*c2dm-sgm*s2dm;
sgmm2g = sgm*c2gs-cgm*s2gs;
cgmm2g = cgm*c2gs+sgm*s2gs;
s2mmgs = s2gm*cgs-c2gm*sgs;
c2mmgs = c2gm*cgs+s2gm*sgs;
s2fm2d = s2fm*c2dm-c2fm*s2dm;
c2fm2d = c2fm*c2dm+s2fm*s2dm;
sgmm2d = sgm*c2dm-cgm*s2dm;
cgmm2d = cgm*c2dm+sgm*s2dm;
s2gm2d = s2gm*c2dm-c2gm*s2dm;
c2gm2d = c2gm*c2dm+s2gm*s2dm;
sfmp2d = sfm*c2dm+cfm*s2dm;
cfmp2d = cfm*c2dm-sfm*s2dm;
sfmm2d = sfm*c2dm-cfm*s2dm;
cfmm2d = cfm*c2dm+sfm*s2dm;
s2gpgs = s2gm*cgs+c2gm*sgs;
c2gpgs = c2gm*cgs-s2gm*sgs;
s2gp2d = s2gm*c2dm+c2gm*s2dm;
c2gp2d = c2gm*c2dm-s2gm*s2dm;
sfmp4d = sfm*c4dm+cfm*s4dm;
cfmp4d = cfm*c4dm-sfm*s4dm;
sfmm4d = sfm*c4dm-cfm*s4dm;
cfmm4d = cfm*c4dm+sfm*s4dm;
*lon = lm + (
22640. * sgm
- 4586. * sgmm2d
+ 2370. * s2dm
+ 769. * s2gm
- 668. * sgs
- 412. * s2fm
- 212. * s2gm2d
- 206. * (sgmpgs*c2dm-cgmpgs*s2dm)
+ 192. * sgmp2d
+ 165. * (s2dm*cgs-c2dm*sgs)
+ 148. * sgmmgs
- 125. * sdm
- 110. * sgmpgs
- 55. * s2fm2d
- 45. * (sgm*c2fm+cgm*s2fm)
+ 40. * (sgm*c2fm-cgm*s2fm)
- 38. * (sgm*c4dm-cgm*s4dm)
+ 36. * s3gm
- 31. * (s2gm*c4dm-c2gm*s4dm)
+ 28. * (sgmmgs*c2dm-cgmmgs*s2dm)
- 24. * (s2dm*cgs+c2dm*sgs)
+ 19. * (sgm*cdm-cgm*sdm)
+ 18. * (sdm*cgs+cdm*sgs)
+ 15. * (sgmmgs*c2dm+cgmmgs*s2dm)
+ 14. * (s2gm*c2dm+c2gm*s2dm)
+ 14. * s4dm
- 13. * (s3gm*c2dm-c3gm*s2dm)
+ (14.2+0.55*t) * sind(venarg+140.72-22.43*t)
+ 10. * (s2gm*cgs-c2gm*sgs)
+ 9. * (sgmm2d*c2fm-cgmm2d*s2fm)
- 9. * (s2gm2d*cgs+c2gm2d*sgs)
- 8. * (sgm*cdm+cgm*sdm)
+ 8. * (s2dm*c2gs-c2dm*s2gs)
- 8. * (s2gm*cgs+c2gm*sgs)
- 7. * s2gs
- 7. * (sgmm2d*c2gs+cgmm2d*s2gs)
+ 7. * snm
- 6. * (sgmp2d*c2fm-cgmp2d*s2fm)
- 6. * (s2fm*c2dm+c2fm*s2dm)
- 4. * (sgmpgs*c4dm-cgmpgs*s4dm)
- 4. * (s2gm*c2fm+c2gm*s2fm)
+ 3. * (sgm*c3dm-cgm*s3dm)
- 3. * (sgmpgs*c2dm+cgmpgs*s2dm)
- 3. * (s2gpgs*c4dm-c2gpgs*s4dm)
+ 3. * sgmm2g
+ 3. * (sgmm2g*c2dm-cgmm2g*s2dm)
- 2. * (s2mmgs*c2dm-c2mmgs*s2dm)
- 2. * (s2fm2d*cgs+c2fm2d*sgs)
+ 2. * (sgm*c4dm+cgm*s4dm)
+ 2. * s4gm
+ 2. * (s4dm*cgs-c4dm*sgs)
+ 2. * (s2gm*cdm-c2gm*sdm)
) / 3600.0;
*lat = (
18461. * sfm
+ 1010. * (sgm*cfm+cgm*sfm)
+ 1000. * (sgm*cfm-cgm*sfm)
- 624. * (sfm*c2dm-cfm*s2dm)
- 199. * (sgmm2d*cfm-cgmm2d*sfm)
- 167. * (sgmm2d*cfm+cgmm2d*sfm)
+ 117. * (sfm*c2dm+cfm*s2dm)
+ 62. * (s2gm*cfm+c2gm*sfm)
+ 33. * (sgmp2d*cfm-cgmp2d*sfm)
+ 32. * (s2gm*cfm-c2gm*sfm)
- 30. * (sfmm2d*cgs+cfmm2d*sgs)
- 16. * (s2gm*cfmm2d+c2gm*sfmm2d)
+ 15. * (sgm*cfmp2d+cgm*sfmp2d)
+ 12. * (sfmm2d*cgs-cfmm2d*sgs)
- 9. * (sgmpgs*cfmp2d-cgmpgs*sfmp2d)
- 8. * (sfm*cnm+cfm*snm)
+ 8. * (sfmp2d*cgs-cfmp2d*sgs)
- 7. * (sgmpgs*cfmm2d+cgmpgs*sfmm2d)
+ 7. * (sgmmgs*cfm+cgmmgs*sfm)
- 7. * (sgm*cfmm4d+cgm*sfmm4d)
- 6. * (sfm*cgs+cfm*sgs)
- 6. * s3fm
+ 6. * (sgmmgs*cfm-cgmmgs*sfm)
- 5. * (sfm*cdm+cfm*sdm)
- 5. * (sgmpgs*cfm+cgmpgs*sfm)
- 5. * (sgmpgs*cfm-cgmpgs*sfm)
+ 5. * (sfm*cgs-cfm*sgs)
+ 5. * (sfm*cdm-cfm*sdm)
+ 4. * (s3gm*cfm+c3gm*sfm)
- 4. * sfmm4d
- 3. * (sgm*cfmp4d-cgm*sfmp4d)
+ 3. * (sgm*c3fm-cgm*s3fm)
- 2. * (s2gm*cfmp4d-c2gm*sfmp4d)
- 2. * (s3fm*c2dm-c3fm*s2dm)
+ 2. * (s2gp2d*cfm-c2gp2d*sfm)
+ 2. * (sgmmgs*cfmm2d-cgmmgs*sfmm2d)
+ 2. * (s2gm*cfmp2d-c2gm*sfmp2d)
+ 2. * (s3gm*cfm-c3gm*sfm)
) / 3600.0;
r = 60.36298
- 3.27746 * cgm
- 0.57994 * cgmm2d
- 0.46357 * c2dm
- 0.08904 * c2gm
+ 0.03865 * c2gm2d
- 0.03237 * (c2dm*cgs+s2dm*sgs)
- 0.02688 * (cgm*c2dm-sgm*s2dm)
- 0.02358 * (cgmm2d*cgs-sgmm2d*sgs)
- 0.02030 * cgmmgs
+ 0.01719 * cdm
+ 0.01671 * cgmpgs
+ 0.01247 * (cgm*c2fm+sgm*s2fm)
+ 0.00704 * cgs
+ 0.00529 * (c2dm*cgs-s2dm*sgs)
- 0.00524 * (cgm*c4dm+sgm*s4dm)
+ 0.00398 * (cgmmgs*c2dm+sgmmgs*s2dm)
- 0.00366 * c3gm
- 0.00295 * (c2gm*c4dm+s2gm*s4dm)
- 0.00263 * (cdm*cgs-sdm*sgs)
+ 0.00249 * (c3gm*c2dm+s3gm*s2dm)
- 0.00221 * (cgmmgs*c2dm-sgmmgs*s2dm)
+ 0.00185 * (c2fm*c2dm+s2fm*s2dm)
- 0.00161 * (c2dm*c2gs+s2dm*s2gs)
+ 0.00147 * (cgmm2d*c2fm-sgmm2d*s2fm)
- 0.00142 * c4dm
+ 0.00139 * (c2gm2d*cgs-s2gm2d*sgs)
- 0.00118 * (cgmpgs*c4dm+sgmpgs*s4dm)
- 0.00116 * (c2gm*c2dm-s2gm*s2dm)
- 0.00110 * (c2gm*cgs+s2gm*sgs) ;
*par = asind( 1.0 / r );
} /* monpos_vf */
==================================cut===================================
/*
VF_MER.C - Mercurys position
Referens: van Flandern & Pulkkinen
Astrophys. J. Suppl., 41:391-411, nov 1979
Paul Schlyter, 1987-06-14
*/
#include <math.h>
#include "trig.h"
static double frac( double x )
/****************************/
{
return( x - floor(x) );
} /* frac */
void merpos_vf( double JD, double *lon, double *lat, double *r )
/**************************************************************/
{
double dnr, t,
l1, g1, f1, sg1, cg1, sf1, cf1, s2f1, c2f1, s3f1, c3f1,
s2g1, c2g1, s3g1, c3g1, s4g1, c4g1, s5g1, c5g1,
g2;
dnr = JD - 2451545.0;
t = dnr / 36525.0 + 1.0;
/* Venus mean anomaly: */
g2 = 360.0 * frac( 0.140023 + 0.00445036173 * dnr );
/* Mercurys heliocentric position: */
l1 = 360.0 * frac( 0.700695 + 0.01136771400 * dnr );
g1 = 360.0 * frac( 0.485541 + 0.01136759566 * dnr );
f1 = 360.0 * frac( 0.566441 + 0.01136762384 * dnr );
sg1 = sind(g1);
cg1 = cosd(g1);
sf1 = sind(f1);
cf1 = cosd(f1);
s2g1 = 2*sg1*cg1;
c2g1 = 2*cg1*cg1-1.0;
s3g1 = s2g1*cg1+c2g1*sg1;
c3g1 = c2g1*cg1-s2g1*sg1;
s4g1 = 2*s2g1*c2g1;
c4g1 = 2*c2g1*c2g1-1.0;
s5g1 = s4g1*cg1+c4g1*sg1;
c5g1 = c4g1*cg1-s4g1*sg1;
s2f1 = 2*sf1*cf1;
c2f1 = 2*cf1*cf1-1.0;
s3f1 = s2f1*cf1+c2f1*sf1;
c3f1 = c2f1*cf1-s2f1*sf1;
*lon = l1 + (
(84378.+8.*t) * sg1
+10733. * s2g1
+ 1892. * s3g1
- 646. * s2f1
+ 381. * s4g1
- 306. * (sg1*c2f1-cg1*s2f1)
- 274. * (sg1*c2f1+cg1*s2f1)
- 92. * (s2g1*c2f1+c2g1*s2f1)
+ 83. * s5g1
- 28. * (s3g1*c2f1+c3g1*s2f1)
+ 25. * (s2g1*c2f1-c2g1*s2f1)
+ 19. * (2*s3g1*c3g1)
- 9. * (s4g1*c2f1+c4g1*s2f1)
+ 7. * cosd(2*g1-5*g2)
) / 3600.0;
*lat = ( 24134. * sf1
+ 5180. * (sg1*cf1-cg1*sf1)
+ 4910. * (sg1*cf1+cg1*sf1)
+ 1124. * (s2g1*cf1+c2g1*sf1)
+ 271. * (s3g1*cf1+c3g1*sf1)
+ 132. * (s2g1*cf1-c2g1*sf1)
+ 67. * (s4g1*cf1+c4g1*sf1)
+ 18. * (s3g1*cf1-c3g1*sf1)
+ 17. * (s5g1*cf1+c5g1*sf1)
- 10. * s3f1
- 9. * (sg1*c3f1-cg1*s3f1)
) / 3600.0;
*r = 0.39528
- 0.07834 * cg1
- 0.00795 * c2g1
- 0.00121 * c3g1
- 0.00022 * c4g1;
} /* merpos_vf */
==================================cut===================================
/*
VF_VEN.C - Venus position
Referens: van Flandern & Pulkkinen
Astrophys. J. Suppl., 41:391-411, nov 1979
Paul Schlyter, 1987-06-14
*/
#include <math.h>
#include "trig.h"
static double frac( double x )
/****************************/
{
return( x - floor(x) );
} /* frac */
void venpos_vf( double JD, double *lon, double *lat, double *r )
/**************************************************************/
{
double dnr, t, gs, sgs, cgs,
l2, g2, f2, sg2, cg2, sf2, cf2, s2g2, c2g2, s3g2, c3g2,
s2gs, c2gs, s3gs, c3gs;
dnr = JD - 2451545.0;
t = dnr / 36525.0 + 1.0;
/* Fundamental arguments Sun: */
gs = 360.0 * frac( 0.993126 + 0.00273777850 * dnr );
sgs = sind(gs);
cgs = cosd(gs);
/* Venus mean anomaly: */
g2 = 360.0 * frac( 0.140023 + 0.00445036173 * dnr );
/* Venus heliocentric position: */
l2 = 360.0 * frac( 0.505498 + 0.00445046867 * dnr );
g2 = 360.0 * frac( 0.140023 + 0.00445036173 * dnr );
f2 = 360.0 * frac( 0.292498 + 0.00445040017 * dnr );
sg2 = sind(g2);
cg2 = cosd(g2);
sf2 = sind(f2);
cf2 = cosd(f2);
s2g2 = 2*sg2*cg2;
c2g2 = 2*cg2*cg2-1.0;
s3g2 = s2g2*cg2+c2g2*sg2;
c3g2 = c2g2*cg2-s2g2*sg2;
s2gs = 2*sgs*cgs;
c2gs = 2*cgs*cgs-1.0;
s3gs = s2gs*cgs+c2gs*sgs;
c3gs = c2gs*cgs-s2gs*sgs;
*lon = l2 + (
(2814.-20.*t) * sg2
- 181. * (2*sf2*cf2)
+ 12. * s2g2
- 10. * (c2gs*c2g2+s2gs*s2g2)
+ 7. * (c3gs*c3g2+s3gs*s3g2)
) / 3600.0;
*lat = ( 12215. * sf2
+ 83. * (sg2*cf2+cg2*sf2)
+ 83. * (sg2*cf2-cg2*sf2)
) / 3600.0;
*r = 0.72335 - 0.00493 * cg2;
} /* venpos_vf */
==================================cut===================================
/*
VF_MAR.C - Mars position
Referens: van Flandern & Pulkkinen
Astrophys. J. Suppl., 41:391-411, nov 1979
Paul Schlyter, 1987-06-14
*/
#include <math.h>
#include "trig.h"
static double frac( double x )
/****************************/
{
return( x - floor(x) );
} /* frac */
void marpos_vf( double JD, double *lon, double *lat, double *r )
/**************************************************************/
{
double dnr, t, gs, sgs, cgs, g5,
g2, sg2, cg2,
s2gs, c2gs,
l4, g4, f4, sg4, cg4, sf4, cf4, sg5, cg5,
s2g4, c2g4, s3g4, c3g4, s4g4, c4g4, s2f4, c2f4, s2g5, c2g5;
dnr = JD - 2451545.0;
t = dnr / 36525.0 + 1.0;
/* Fundamental arguments Sun: */
gs = 360.0 * frac( 0.993126 + 0.00273777850 * dnr );
sgs = sind(gs);
cgs = cosd(gs);
/* Frequently used perturbation arguments: */
/* Venus: */
g2 = 360.0 * frac( 0.140023 + 0.00445036173 * dnr );
/* Mars: */
g4 = 360.0 * frac( 0.053856 + 0.00145561327 * dnr );
/* Jupiter: */
g5 = 360.0 * frac( 0.056531 + 0.00023080893 * dnr );
/* Mars heliocentric position: */
l4 = 360.0 * frac( 0.987353 + 0.00145575328 * dnr );
g4 = 360.0 * frac( 0.053856 + 0.00145561327 * dnr );
f4 = 360.0 * frac( 0.849694 + 0.00145569465 * dnr );
sg4 = sind(g4);
cg4 = cosd(g4);
sf4 = sind(f4);
cf4 = cosd(f4);
sg2 = sind(g2);
cg2 = cosd(g2);
sg5 = sind(g5);
cg5 = cosd(g5);
s2gs = 2*sgs*cgs;
c2gs = 2*cgs*cgs-1.0;
s2g4 = 2*sg4*cg4;
c2g4 = 2*cg4*cg4-1.0;
s3g4 = s2g4*cg4+c2g4*sg4;
c3g4 = c2g4*cg4-s2g4*sg4;
s4g4 = 2*s2g4*c2g4;
c4g4 = 2*c2g4*c2g4-1.0;
s2f4 = 2*sf4*cf4;
c2f4 = 2*cf4*cf4-1.0;
s2g5 = 2*sg5*cg5;
c2g5 = 2*cg5*cg5-1.0;
*lon = l4 + (
(38451.+37.*t) * sg4
+ (2238.+4.*t) * s2g4
+ 181. * s3g4
- 52. * s2f4
- 22. * (cg4*c2g5+sg4*s2g5)
- 19. * (sg4*cg5-cg4*sg5)
+ 17. * (cg4*cg5+sg4*sg5)
+ 17. * s4g4
- 16. * (c2g4*c2g5+s2g4*s2g5)
+ 13. * (cgs*c2g4+sgs*s2g4)
- 10. * (sg4*c2f4-cg4*s2f4)
- 10. * (sg4*c2f4+cg4*s2f4)
+ 7. * (cgs*cg4+sgs*sg4)
- 7. * (c2gs*c3g4+s2gs*s3g4)
- 5. * (sg2*c3g4-cg2*s3g4)
- 5. * (sgs*cg4-cgs*sg4)
- 5. * (sgs*c2g4-cgs*s2g4)
- 4. * (c2gs*c4g4+s2gs*s4g4)
+ 4. * cg5
+ 3. * (cg2*c3g4+sg2*s3g4)
+ 3. * (s2g4*c2g5-c2g4*s2g5)
) / 3600.0;
*lat = ( 6603. * sf4
+ 622. * (sg4*cf4-cg4*sf4)
+ 615. * (sg4*cf4+cg4*sf4)
+ 64. * (s2g4*cf4+c2g4*sf4)
) / 3600.0;
*r = 1.53031
- 0.14170 * cg4
- 0.00660 * c2g4
- 0.00047 * c3g4;
} /* marpos_vf */
==================================cut===================================
/*
VF_JUP.C - Jupiters position
Referens: van Flandern & Pulkkinen
Astrophys. J. Suppl., 41:391-411, nov 1979
Paul Schlyter, 1987-06-14
*/
#include <math.h>
#include "trig.h"
static double frac( double x )
/****************************/
{
return( x - floor(x) );
} /* frac */
void juppos_vf( double JD, double *lon, double *lat, double *r )
/**************************************************************/
{
double dnr, t, g5, g6, g7,
sg5, cg5, sg6, cg6, sg7, cg7,
s2g5, c2g5, s3g5, c3g5, s4g5, c4g5, s5g5, c5g5,
s2g6, c2g6, s3g6, c3g6, s4g6, c4g6, s5g6, c5g6,
s6g6, c6g6, s10g6, c10g6, /* s7g6, c7g6, s9g6, c9g6, s11g6, c11g6, */
s2g7, c2g7, s3g7, c3g7,
/*s2gm4g,*/ c2gm4g, s2gm5g, c2gm5g, s2gm6g, c2gm6g,
l5, sl5, cl5, s2l5, c2l5;
dnr = JD - 2451545.0;
t = dnr/36525.+1.0;
g5 = 360.0 * frac( 0.056531 + 0.00023080893 * dnr );
g6 = 360.0 * frac( 0.882987 + 0.00009294371 * dnr );
g7 = 360.0 * frac( 0.400589 + 0.00003269438 * dnr );
sg5 = sind(g5);
cg5 = cosd(g5);
sg6 = sind(g6);
cg6 = cosd(g6);
sg7 = sind(g7);
cg7 = cosd(g7);
s2g5 = 2*sg5*cg5;
c2g5 = 2*cg5*cg5-1.0;
s3g5 = s2g5*cg5+c2g5*sg5;
c3g5 = c2g5*cg5-s2g5*sg5;
s4g5 = 2*s2g5*c2g5;
c4g5 = 2*c2g5*c2g5-1.0;
s5g5 = s4g5*cg5+c4g5*sg5;
c5g5 = c4g5*cg5-s4g5*sg5;
s2g6 = 2*sg6*cg6;
c2g6 = 2*cg6*cg6-1.0;
s3g6 = s2g6*cg6+c2g6*sg6;
c3g6 = c2g6*cg6-s2g6*sg6;
s4g6 = 2*s2g6*c2g6;
c4g6 = 2*c2g6*c2g6-1.0;
s5g6 = s4g6*cg6+c4g6*sg6;
c5g6 = c4g6*cg6-s4g6*sg6;
s6g6 = 2*s3g6*c3g6;
c6g6 = 2*c3g6*c3g6-1.0;
/* s7g6 = s6g6*cg6+c6g6*sg6; */
/* c7g6 = c6g6*cg6-s6g6*sg6; */
/* s9g6 = s7g6*c2g6+c7g6*s2g6; */
/* c9g6 = c7g6*c2g6-s7g6*s2g6; */
s10g6 = 2*s5g6*c5g6;
c10g6 = 2*c5g6*c5g6-1.0;
/* s11g6 = s10g6*cg6+c10g6*sg6; */
/* c11g6 = c10g6*cg6-s10g6*sg6; */
s2g7 = 2*sg7*cg7;
c2g7 = 2*cg7*cg7-1.0;
s3g7 = s2g7*cg7+c2g7*sg7;
c3g7 = c2g7*cg7-s2g7*sg7;
/* s2gm4g = s2g5*c4g6-c2g5*s4g6; */
c2gm4g = c2g5*c4g6+s2g5*s4g6;
s2gm5g = s2g5*c5g6-c2g5*s5g6;
c2gm5g = c2g5*c5g6+s2g5*s5g6;
s2gm6g = s2g5*c6g6-c2g5*s6g6;
c2gm6g = c2g5*c6g6+s2g5*s6g6;
l5 = varv180(g5 + 11.9077);
sl5 = sind(l5);
cl5 = cosd(l5);
s2l5 = 2*sl5*cl5;
c2l5 = 2*cl5*cl5-1.0;
*lon = l5 + (
(19934.+68.*t) * sg5
+ 5023. * t
+ 2511.
+ (1093.-19.*t) * c2gm5g
+ (601.+3.*t) * s2g5
- (479.+43.*t) * s2gm5g
- 185. * (s2g5*c2g6-c2g5*s2g6)
+ (137.-2.*t) * (s3g5*c5g6-c3g5*s5g6)
- 131. * (sg5*c2g6-cg5*s2g6)
+ 79. * (cg5*cg6+sg5*sg6)
- 76. * (c2g5*c2g6+s2g5*s2g6)
- (37.+74.*t) * cg5
+ 66. * (c2g5*c3g6+s2g5*s3g6)
+ 63. * (c3g5*c5g6+s3g5*s5g6)
+ 53. * (cg5*c5g6+sg5*s5g6)
+ 49. * (s2g5*c3g6-c2g5*s3g6)
+ 25. * s2l5
+ 25. * s3g5
- (23.+2.*t) * (sg5*c5g6-cg5*s5g6)
+ 17. * c2gm4g
+ 17. * (c3g5*c3g6+s3g5*s3g6)
- 14. * (sg5*cg6-cg5*sg6)
- 13. * (s3g5*c4g6-c3g5*s4g6)
- 9. * c2l5
+ 9. * cg6
- 9. * sg6
- 9. * (s3g5*c2g6-c3g5*s2g6)
+ 9. * (s4g5*c5g6-c4g5*s5g6)
+ 9. * (s2gm6g*c3g7+c2gm6g*s3g7)
- 8. * (c4g5*c10g6+s4g5*s10g6)
+ 7. * (c3g5*c4g6+s3g5*s4g6)
- 7. * (cg5*c3g6+sg5*s3g6)
- 7. * (s4g5*c10g6-c4g5*s10g6)
- 7. * (sg5*c3g6-cg5*s3g6)
+ 6. * (c4g5*c5g6+s4g5*s5g6)
- 6. * (s3g5*c3g6-c3g5*s3g6)
+ 5. * c2g6
- 4. * (s4g5*c4g6-c4g5*s4g6)
- 4. * c3g6
+ 4. * (c2g5*cg6+s2g5*sg6)
- 4. * (c3g5*c2g6+s3g5*s2g6)
- 4.*t * c2g5
+ 3. * c5g6
+ 3. * (c5g5*c10g6+s5g5*s10g6)
+ 3. * s2g6
- 2. * (s2l5*cg5-c2l5*sg5)
+ 2. * (s2l5*cg5+c2l5*sg5)
) / 3600.0;
*lat = (
(-4692.+21.*t) * cg5
+ (259.+30.*t) * sg5
+ 227.
- 227. * c2g5
+ 16. * (s3g5*c5g6-c3g5*s5g6)
- 13. * (sg5*c5g6-cg5*s5g6)
- 12. * c3g5
+ 12. * s2g5
+ 7. * (c3g5*c5g6+s3g5*s5g6)
- 5. * (cg5*c5g6+sg5*s5g6)
) / 3600.0;
*r = 5.20883
- (0.25122+0.00084*t) * cg5
- 0.00604 * c2g5
+ 0.00260 * (c2g5*c2g6+s2g5*s2g6)
- 0.00170 * (c3g5*c5g6+s3g5*s5g6)
- 0.00106 * (s2g5*c2g6-c2g5*s2g6)
- (0.00046+0.00091*t) * sg5
+ 0.00069 * (s2g5*c3g6-c2g5*s3g6)
- 0.00067 * (sg5*c5g6-cg5*s5g6)
+ 0.00066 * (s3g5*c5g6-c3g5*s5g6)
+ 0.00063 * (sg5*cg6-cg5*sg6)
- 0.00051 * (c2g5*c3g6+s2g5*s3g6)
- 0.00029 * (cg5*c5g6+sg5*s5g6)
+ 0.00027 * (cg5*c2g6+sg5*s2g6)
- 0.00022 * c3g5
- 0.00021 * s2gm5g;
} /* juppos_vf */
==================================cut===================================
/*
VF_SAT.C - Saturnus position
Referens: van Flandern & Pulkkinen
Astrophys. J. Suppl., 41:391-411, nov 1979
Paul Schlyter, 1987-06-14
*/
#include <math.h>
#include "trig.h"
static double frac( double x )
/****************************/
{
return( x - floor(x) );
} /* frac */
void satpos_vf( double JD, double *lon, double *lat, double *r )
/**************************************************************/
{
double dnr, t, g5, g6, g7,
sg5, cg5, sg6, cg6, sg7, cg7,
s2g5, c2g5, s3g5, c3g5, s4g5, c4g5, /* s5g5, c5g5, */
s2g6, c2g6, s3g6, c3g6, s4g6, c4g6, s5g6, c5g6,
s6g6, c6g6, s7g6, c7g6, s9g6, c9g6, s10g6, c10g6, s11g6, c11g6,
s2g7, c2g7, s3g7, c3g7,
s2gm4g, c2gm4g, s2gm5g, c2gm5g, s2gm6g, c2gm6g,
l6, sl6, cl6, s2l6, c2l6, s2gm7g, c2gm7g,
lona, lonb;
dnr = JD - 2451545.0;
t = dnr/36525.+1.0;
g5 = 360.0 * frac( 0.056531 + 0.00023080893 * dnr );
g6 = 360.0 * frac( 0.882987 + 0.00009294371 * dnr );
g7 = 360.0 * frac( 0.400589 + 0.00003269438 * dnr );
sg5 = sind(g5);
cg5 = cosd(g5);
sg6 = sind(g6);
cg6 = cosd(g6);
sg7 = sind(g7);
cg7 = cosd(g7);
s2g5 = 2*sg5*cg5;
c2g5 = 2*cg5*cg5-1.0;
s3g5 = s2g5*cg5+c2g5*sg5;
c3g5 = c2g5*cg5-s2g5*sg5;
s4g5 = 2*s2g5*c2g5;
c4g5 = 2*c2g5*c2g5-1.0;
/* s5g5 = s4g5*cg5+c4g5*sg5; */
/* c5g5 = c4g5*cg5-s4g5*sg5; */
s2g6 = 2*sg6*cg6;
c2g6 = 2*cg6*cg6-1.0;
s3g6 = s2g6*cg6+c2g6*sg6;
c3g6 = c2g6*cg6-s2g6*sg6;
s4g6 = 2*s2g6*c2g6;
c4g6 = 2*c2g6*c2g6-1.0;
s5g6 = s4g6*cg6+c4g6*sg6;
c5g6 = c4g6*cg6-s4g6*sg6;
s6g6 = 2*s3g6*c3g6;
c6g6 = 2*c3g6*c3g6-1.0;
s7g6 = s6g6*cg6+c6g6*sg6;
c7g6 = c6g6*cg6-s6g6*sg6;
s9g6 = s7g6*c2g6+c7g6*s2g6;
c9g6 = c7g6*c2g6-s7g6*s2g6;
s10g6 = 2*s5g6*c5g6;
c10g6 = 2*c5g6*c5g6-1.0;
s11g6 = s10g6*cg6+c10g6*sg6;
c11g6 = c10g6*cg6-s10g6*sg6;
s2g7 = 2*sg7*cg7;
c2g7 = 2*cg7*cg7-1.0;
s3g7 = s2g7*cg7+c2g7*sg7;
c3g7 = c2g7*cg7-s2g7*sg7;
s2gm4g = s2g5*c4g6-c2g5*s4g6;
c2gm4g = c2g5*c4g6+s2g5*s4g6;
s2gm5g = s2g5*c5g6-c2g5*s5g6;
c2gm5g = c2g5*c5g6+s2g5*s5g6;
s2gm6g = s2g5*c6g6-c2g5*s6g6;
c2gm6g = c2g5*c6g6+s2g5*s6g6;
l6 = varv180(g6 + 90.1109);
sl6 = sind(l6);
cl6 = cosd(l6);
s2l6 = 2*sl6*cl6;
c2l6 = 2*cl6*cl6-1.0;
s2gm7g = s2g5*c7g6-c2g5*s7g6;
c2gm7g = c2g5*c7g6+s2g5*s7g6;
lona =
(23045.-142.*t) * sg6
+ 5014.*t
- (2689.-60.*t) * c2gm5g
+ 2507.
+ (1177.+101.*t) * s2gm5g
- (826.-3.*t) * c2gm4g
+ (802.-11.*t) * s2g6
+ (425.+2.*t) * (sg5*c2g6-cg5*s2g6)
- (114.+229.*t) * cg6
- (153.-3.*t) * (c2g5*c6g6+s2g5*s6g6)
- (70.+3.*t) * c2l6
+ (67.-3.*t) * s2l6
+ (66.+6.*t) * s2gm6g
+ 41. * (sg5*c3g6-cg5*s3g6)
+ 39. * s3g6
+ 31. * (sg5*cg6-cg5*sg6)
+ 31. * (s2g5*c2g6-c2g5*s2g6)
- 29. * (c2g5*c3g6+s2g5*s3g6)
- 28. * (s2gm6g*c3g7+c2gm6g*s3g7)
+ 28. * (cg5*c3g6+sg5*s3g6)
- (12.-22.*t) * s2gm4g
- 22. * (sg6*c3g7-cg6*s3g7)
+ 20. * (s2g5*c3g6-c2g5*s3g6)
+ (20.+6.*t) * (c4g5*c10g6+s4g5*s10g6)
+ 19. * (c2g6*c3g7+s2g6*s3g7)
+ (19.-6.*t) * (s4g5*c10g6-c4g5*s10g6)
- (7.+17.*t) * c2g6
- 16. * (cg6*c3g7+sg6*s3g7);
lonb =
+ 12. * cg5
- 12. * (s2g6*c2g7-c2g6*s2g7)
- 11. * c2gm7g
+ 10. * (s2g6*c3g7-c2g6*s3g7)
+ 10. * (c2g5*c2g6+s2g5*s2g6)
+ 9. * (s4g5*c9g6-c4g5*s9g6)
- 8. * (sg6*c2g7-cg6*s2g7)
- 8. * (c2l6*cg6-s2l6*sg6)
+ 8. * (c2l6*cg6+s2l6*sg6)
+ 8. * (cg6*cg7+sg6*sg7)
- 8. * (s2l6*cg6-c2l6*sg6)
+ 7. * (s2l6*cg6+c2l6*sg6)
- (7.-4.*t) * (cg5*c2g6+sg5*s2g6)
- 5. * (s3g5*c7g6-c3g5*s7g6)
- 5. * (c3g5*c3g6+s3g5*s3g6)
- 5. * (c2g6*c2g7+s2g6*s2g7)
+ 5. * (s3g5*c4g6-c3g5*s4g6)
+ 5. * s2gm7g
+ 4. * (s3g5*c3g6-c3g5*s3g6)
+ 4. * (s3g5*c5g6-c3g5*s5g6)
+ 3. * (c2gm6g*c3g7-s2gm6g*s3g7)
+ 3. * (c3g5*c7g6+s3g5*s7g6)
+ 3. * (c4g5*c9g6+s4g5*s9g6)
+ 3. * (s3g5*c6g6-c3g5*s6g6)
+ 3. * (s2g5*cg6-c2g5*sg6)
+ 3. * (sg5*c4g6-cg5*s4g6)
+ 2. * (c3g6*c3g7+s3g6*s3g7)
+ 2. * s4g6
- 2. * (c3g5*c4g6+s3g5*s4g6)
- 2. * (c2g5*cg6+s2g5*sg6)
- 2. * (s2gm7g*c3g7+c2gm7g*s3g7)
+ 2. * (cg5*c4g6+sg5*s4g6)
+ 2. * (c4g5*c11g6+s4g5*s11g6)
- 2. * (sg6*cg7-cg6*sg7);
*lon = l6 + ( lona + lonb ) / 3600.0;
*lat = (
(8297.+18.*t) * sg6
- (3346.-79.*t) * cg6
+ (462.-4.*t) * s2g6
- 189. * c2g6
+ 185.-10.*t
- 71. * c2gm4g
+ 3.*t * s2gm4g
+ 46. * s2gm6g
- 45. * c2gm6g
+ 29. * s3g6
- 20. * (c2g5*c3g6+s2g5*s3g6)
- 14. * c2gm5g
- 11. * c3g6
+ 9. * (sg5*c3g6-cg5*s3g6)
+ 8. * (sg5*cg6-cg5*sg6)
- 6. * (s2g5*c3g6-c2g5*s3g6)
+ 5. * s2gm7g
- 5. * c2gm7g
+ 4. * s2gm5g
- 3. * (cg5*cg6+sg5*sg6)
+ 3. * (cg5*c3g6+sg5*s3g6)
+ 3. * (sg5*c2g6-cg5*s2g6)
+ 2. * s4g6
- 2. * (c2g5*c2g6+s2g5*s2g6)
) / 3600.0;
*r =
(9.55774-0.00028*t)
- (0.53252-0.00328*t) * cg6
- 0.01878 * s2gm4g
- 0.01482 * c2g6
+ 0.00817 * (sg5*cg6-cg5*sg6)
- 0.00539 * (cg5*c2g6+sg5*s2g6)
- (0.00225+0.00524*t) * sg6
+ 0.00349 * s2gm5g
+ 0.00347 * s2gm6g
+ 0.00149 * (c2g5*c6g6+s2g5*s6g6)
- 0.00126 * (c2g5*c2g6+s2g5*s2g6)
+ 0.00104 * (cg5*cg6+sg5*sg6)
+ 0.00101 * c2gm5g
+ 0.00098 * (cg5*c3g6+sg5*s3g6)
- 0.00073 * (c2g5*c3g6+s2g5*s3g6)
- 0.00062 * c3g6
+ 0.00042 * (s2g6*c3g7-c2g6*s3g7)
+ 0.00041 * (s2g5*c2g6-c2g5*s2g6)
- 0.00040 * (sg5*c3g6-cg5*s3g6)
+ 0.00040 * (c2g5*c4g6+s2g5*s4g6)
- 0.00023 * sg5
+ 0.00020 * (s2g5*c7g6-c2g5*s7g6);
} /* satpos_vf */
==================================cut===================================
/*
VF_URA.C - Uranus position
Referens: van Flandern & Pulkkinen
Astrophys. J. Suppl., 41:391-411, nov 1979
Paul Schlyter, 1987-06-14
*/
#include <math.h>
#include "trig.h"
static double frac( double x )
/****************************/
{
return( x - floor(x) );
} /* frac */
void urapos_vf( double jd, double *lon, double *lat, double *r )
/**************************************************************/
{
double dnr, t,
g5, g6, g7, g8,
sg5, cg5, sg6, cg6, sg7, cg7, sg8, cg8,
s2g5, c2g5,
s2g6, c2g6, s3g6, c3g6, s6g6, c6g6,
s2g7, c2g7, s3g7, c3g7,
s2g8, c2g8, s3g8, c3g8,
s2gm6g, c2gm6g,
l7, f7, sf7, cf7, s4g7, c4g7, s2f7, s4g8, c4g8;
dnr = jd - 2451545.0;
t = dnr/36525.0 + 1.0;
g5 = 360.0 * frac( 0.056531 + 0.00023080893 * dnr );
g6 = 360.0 * frac( 0.882987 + 0.00009294371 * dnr );
g7 = 360.0 * frac( 0.400589 + 0.00003269438 * dnr );
g8 = 360.0 * frac( 0.725368 + 0.00001672092 * dnr );
sg5 = sind(g5);
cg5 = cosd(g5);
sg6 = sind(g6);
cg6 = cosd(g6);
sg7 = sind(g7);
cg7 = cosd(g7);
sg8 = sind(g8);
cg8 = cosd(g8);
s2g5 = 2*sg5*cg5;
c2g5 = 2*cg5*cg5-1.0;
s2g6 = 2*sg6*cg6;
c2g6 = 2*cg6*cg6-1.0;
s3g6 = s2g6*cg6+c2g6*sg6;
c3g6 = c2g6*cg6-s2g6*sg6;
s6g6 = 2*s3g6*c3g6;
c6g6 = 2*c3g6*c3g6-1.0;
s2g7 = 2*sg7*cg7;
c2g7 = 2*cg7*cg7-1.0;
s3g7 = s2g7*cg7+c2g7*sg7;
c3g7 = c2g7*cg7-s2g7*sg7;
s2g8 = 2*sg8*cg8;
c2g8 = 2*cg8*cg8-1.0;
s3g8 = s2g8*cg8+c2g8*sg8;
c3g8 = c2g8*cg8-s2g8*sg8;
s2gm6g = s2g5*c6g6-c2g5*s6g6;
c2gm6g = c2g5*c6g6+s2g5*s6g6;
/* Uranus position: */
l7 = varv180(g7 + 169.0488);
f7 = 360.0 * frac( 0.664614 + 0.00003265562 * dnr );
sf7 = sind(f7);
cf7 = cosd(f7);
s4g7 = 2*s2g7*c2g7;
c4g7 = 2*c2g7*c2g7-1.0;
s2f7 = 2*sf7*cf7;
s4g8 = 2*s2g8*c2g8;
c4g8 = 2*c2g8*c2g8-1.0;
*lon = l7 + (
(19397.+110.*t-9.*t*t) * sg7
+ (570.+7.*t) * s2g7
- (12.+536.*t+12.*t*t) * cg7
+ 143. * (sg6*c2g7-cg6*s2g7)
+ (102.-7.*t) * (sg6*c3g7-cg6*s3g7)
+ (76.+7.*t) * (cg6*c3g7+sg6*s3g7)
- 49. * (sg5*cg7-cg5*sg7)
+ 32.*t*t
- 30.*t * c2g7
+ 29. * (s2gm6g*c3g7+c2gm6g*s3g7)
+ 29. * (c2g7*c2g8+s2g7*s2g8)
- 28. * (cg7*cg8+sg7*sg8)
+ 23. * s3g7
- 21. * (cg5*cg7+sg5*sg7)
+ 20. * (sg7*cg8-cg7*sg8)
+ (20.+8.*t) * (cg6*c2g7+sg6*s2g7)
- 19. * (cg6*cg7+sg6*sg7)
+ 17. * (s2g7*c3g8-c2g7*s3g8)
+ 14. * (s3g7*c3g8-c3g7*s3g8)
+ 13. * (sg6*cg7-cg6*sg7)
+ 10. * (s2g7*c2g8-c2g7*s2g8)
- 9. * s2f7
+ 9. * (c2g7*c3g8+s2g7*s3g8)
+ 6. * (s2gm6g*c2g7+c2gm6g*s2g7)
+ 6. * (c2gm6g*c2g7-s2gm6g*s2g7)
+ 5. * (sg6*c4g7-cg6*s4g7)
- 4. * (s3g7*c4g8-c3g7*s4g8)
+ 4. * (c3g7*c3g8+s3g7*s3g8)
- 3. * cg8
- 2. * sg8
) / 3600.0;
*lat = ( 2775. * sf7
+ 131. * (sg7*cf7-cg7*sf7)
+ 130. * (sg7*cf7+cg7*sf7)
) / 3600.0;
*r = 19.21216
- (0.90154+0.00508*t) * cg7
- 0.02488*t * sg7
- 0.02121 * c2g7
- 0.00585 * (cg6*c2g7+sg6*s2g7)
- 0.00451 * (cg5*cg7+sg5*sg7)
+ 0.00336 * (sg6*cg7-cg6*sg7)
+ 0.00198 * (sg5*cg7-cg5*sg7)
+ 0.00118 * (cg6*c3g7+sg6*s3g7)
+ 0.00107 * (sg6*c2g7-cg6*s2g7)
- 0.00103*t * s2g7
- 0.00081 * (c3g7*c3g8+s3g7*s3g8);
} /* urapos_vf */
==================================cut===================================
/*
VF_NEP.C - Neptunus position
Referens: van Flandern & Pulkkinen
Astrophys. J. Suppl., 41:391-411, nov 1979
Paul Schlyter, 1987-06-14
*/
#include <math.h>
#include "trig.h"
static double frac( double x )
/****************************/
{
return( x - floor(x) );
} /* frac */
void neppos_vf( double jd, double *lon, double *lat, double *r )
/**************************************************************/
{
double dnr, t,
g5, g6, g7, g8,
sg5, cg5, sg6, cg6, sg7, cg7, sg8, cg8,
s2g7, c2g7,
s2g8, c2g8, s3g8, c3g8,
l7,
l8, f8, sf8, cf8;
dnr = jd - 2451545.0;
t = dnr/36525.0 + 1.0;
g5 = 360.0 * frac( 0.056531 + 0.00023080893 * dnr );
g6 = 360.0 * frac( 0.882987 + 0.00009294371 * dnr );
g7 = 360.0 * frac( 0.400589 + 0.00003269438 * dnr );
g8 = 360.0 * frac( 0.725368 + 0.00001672092 * dnr );
sg5 = sind(g5);
cg5 = cosd(g5);
sg6 = sind(g6);
cg6 = cosd(g6);
sg7 = sind(g7);
cg7 = cosd(g7);
sg8 = sind(g8);
cg8 = cosd(g8);
s2g7 = 2*sg7*cg7;
c2g7 = 2*cg7*cg7-1.0;
s2g8 = 2*sg8*cg8;
c2g8 = 2*cg8*cg8-1.0;
s3g8 = s2g8*cg8+c2g8*sg8;
c3g8 = c2g8*cg8-s2g8*sg8;
/* Neptunus position: */
l7 = varv180(g7 + 169.0488);
l8 = varv180(g8 + 43.7558);
f8 = 360.0 * frac( 0.480856 + 0.00001663715 * dnr );
sf8 = sind(f8);
cf8 = cosd(f8);
*lon = l8 + (
(3523.-4.*t+4.*t*t) * sg8
- 50. * 2*sf8*cf8
- 43.*t * cg8
+ 29. * (sg5*cg8-cg5*sg8)
+ 19. * s2g8
- 18. * (cg5*cg8+sg5*sg8)
+ 13. * (cg6*cg8+sg6*sg8)
+ 13. * (sg6*cg8-cg6*sg8)
- 9. * (s2g7*c3g8-c2g7*s3g8)
+ 9. * (c2g7*c2g8+s2g7*s2g8)
- 5. * (c2g7*c3g8+s2g7*s3g8)
+ 4. * (cg7*c2g8+sg7*s2g8)
) / 3600.0;
*lat = ( (6404.-33.*t) * sf8
+ 55. * (sg8*cf8+cg8*sf8)
+ 55. * (sg8*cf8-cg8*sf8)
) / 3600.0;
*r = 30.07175
- 0.25701 * cg8
- 0.00787 * cosd(2*l7-g7-2*l8)
+ 0.00409 * (cg5*cg8+sg5*sg8)
- 0.00314*t * sg8
+ 0.00250 * (sg5*cg8-cg5*sg8)
- 0.00194 * (sg6*cg8-cg6*sg8)
+ 0.00185 * (cg6*cg8+sg6*sg8);
} /* neppos_vf */
==================================cut===================================
/*
PLTPOSVF.C
Plutos position
Referens: van Flandern & Pulkkinen
Astrophys. J. Suppl., 41:391-411, nov 1979
Paul Schlyter, 1987-06-14
*/
#include <math.h>
#include "trig.h"
static double frac( double x )
/****************************/
{
return( x - floor(x) );
} /* frac */
void plutpos_vf( double jd, double *lon, double *lat, double *r )
/***************************************************************/
{
double dnr, t,
l9, g9, f9,
sg9, cg9, sf9, cf9,
s2g9, c2g9, s3g9, c3g9, s4g9, c4g9,
s2f9, c2f9, s3f9, c3f9;
dnr = jd - 2451545.0;
t = dnr/36525.0 + 1.0;
l9 = 360.0 * frac( 0.663854 + 0.00001115482 * dnr );
g9 = 360.0 * frac( 0.041020 + 0.00001104864 * dnr );
f9 = varv180(g9 + 113.8806);
sg9 = sind(g9);
cg9 = cosd(g9);
sf9 = sind(f9);
cf9 = cosd(f9);
s2g9 = 2*sg9*cg9;
c2g9 = 2*cg9*cg9-1.0;
s3g9 = s2g9*cg9+c2g9*sg9;
c3g9 = c2g9*cg9-s2g9*sg9;
s4g9 = 2*s2g9*c2g9;
c4g9 = 2*c2g9*c2g9-1.0;
s2f9 = 2*sf9*cf9;
c2f9 = 2*cf9*cf9-1.0;
s3f9 = s2f9*cf9+c2f9*sf9;
c3f9 = c2f9*cf9-s2f9*sf9;
*lon = l9 + (
(101577.+200.*t+227.*t*t) * sg9
+ 15517. * s2g9
- 3593. * s2f9
+ 3414. * s3g9
- 2201. * (sg9*c2f9-cg9*s2f9)
- 1871. * (sg9*c2f9+cg9*s2f9)
+ 839. * s4g9
- 757. * (s2g9*c2f9+c2g9*s2f9)
- 285. * (s3g9*c2f9+c3g9*s2f9)
+ 218. * (s2g9*c2f9-c2g9*s2f9)
) / 3600.0;
*lat = (
57726. * sf9
+ 15257. * (sg9*cf9-cg9*sf9)
+ 14102. * (sg9*cf9+cg9*sf9)
+ 3870. * (s2g9*cf9+c2g9*sf9)
+ 1138. * (s3g9*cf9+c3g9*sf9)
+ 472. * (s2g9*cf9-c2g9*sf9)
+ 353. * (s4g9*cf9+c4g9*sf9)
- 144. * (sg9*c3f9-cg9*s3f9)
- 119. * s3f9
- 111. * (sg9*c3f9+cg9*s3f9)
) / 3600.0;
*r = 40.74638
- 9.58235 * cg9
- 1.16703 * c2g9
- 0.22649 * c3g9
- 0.04996 * c4g9;
} /* plutpos_vf */
==================================cut===================================
--
----------------------------------------------------------------
Paul Schlyter, Swedish Amateur Astronomer's Society (SAAF)
Grev Turegatan 40, S-114 38 Stockholm, SWEDEN
e-mail: pau...@saaf.se paul.s...@ausys.se pa...@inorbit.com
WWW: http://hotel04.ausys.se/pausch http://welcome.to/pausch
Joe Heafner
>
> I'm implementing and extending the methods found in "Low-precision
> Formulae for Planetary Positions" by T.C. Van Flandern and K. F.
> Pulkkinen (U.S. Naval Observatory, 1979) in Java. I've seen other works
> based on this paper, but none that implement the series to their full
> precision -- within one arc-second for +/- 300 years. It's a lot of
> tedious typing to get all of those series terms in that paper entered.
> So far I've got the data entered for the sun and moon, as well as
> Mercury, Venus, and Mars. I've been spending more time coding in Java
> than data entry, but I should soon have all of the planets out to Pluto
> entered.
>
Pulkkinen series on hand for years. It took some time to code them all
up, but I did it. I have a C version stashed away somewhere too.
> My goal is to create a flexible, extensible set of classes and methods
> that are easily incorporated into any Java program. Beyond computing the
> planetary elements themselves as given in the Van Flandern/Pulkkinen
> paper, I've developed methods for computing moon phase events, solstice
> and equinox events, rise and set times, and a number of utility
> functions. In the future, I'd like to add elongation and opposition
> events, Jovian moons, and any other good suggestions that come my way,
> that I feel capable of implementing. Once the code is further developed,
> I intend to make it all publicly available.
>
BEGIN SHAMELESS PLUG
You may be interested in my book, PRACTICAL EPHEMERIS COMPUTATIONS, to
be published in February or March by Willmann-Bell. I illustrate all the
various reduction algorithms needed to produce accurate planetary,
solar, and lunar ephemerides. Furthermore, I provide a software toolbox
in PowerBasic and C that allows the user to build powerful applications.
My software is designed to work with the JPL planetary ephemeris data
files (also available on CD-ROM from Willmann-Bell). While not written
in Java, you may get some ideas for algorithms from my book.
END SHAMELESS PLUG
Thanks,
Joe Heafner
Astronomy/Physics Instructor
Kerry Shetline
Paul Schlyter wrote:
> Perhaps I could save you some time. Below are Van Flandern-Pulkkinen's
> series, coded by me in C, back in 1987. I only coded the heliocentric
> series, since the geocentric series were both longer and less accurate,
> and one can compute the geocentric positions from the heliocentric
> positions anyway.
Rather than putting the series directly into code, the data I use are
external to my program so that they can easily be modified or replaced
if necessary. A section of the data file appears as follows:
MOON
60.40974
1
1
22640 0 S 2,1
-4586 0 S 2,1 4,-2
2370 0 S 4,2
769 0 S 2,2
-668 0 S 8,1
-412 0 S 3,2
-212 0 S 2,2 4,-2
-206 0 S 2,1 4,-2 8,1
There's an object called VPPlanet which contains the name of the planet,
a scaling factor (used for the U, V, and W series), the indices of the
mean elements to use for mean longitude and right ascension, and an
array of six VPSeries objects. Each VPSeries object corresponds to one
term in the series and forms a linked list to the subsequent terms. Each
term is expressed as a coefficient, the exponent x for T^x, a flag for
whether to compute the sine or the cosine, and a linked list of factors.
Each factor is a coefficient and an index into the mean elements.
>I ended up extensively using only the series for
> Jupiter and Saturn; for the other celestial bodies I picked the series
> by Jean Meeus in his "Astronomical formulae for Calculators" instead,
> which I found to be slightly more accurate.
I'll take it that if the formulae are meant for calculators, they're
probably much shorter and faster to compute. If on top of that the
formulae are actually more accurate, that will be quite a bonus.
By the way, I just checked amazon.com for the book. It's a special order
and will take 4-6 weeks, unfortunately. Guess it doesn't rank in sales
with the latest Stephen King :) By the way, there's a newer book by Jean
Meeus, "Astronomical Algorithms", that came out in 1998 -- that looks
like the way to go now.
> Then you may need to learn some FORTRAN, since there are implementations
> in FORTRAN of these series, to their full precision, on the Net. All
> the ASCII source code is some 35 MBytes in size....
I haven't used FORTRAN since the late 1970s, but I'm sure I can still
wade through it. I think I've found what may be the APIs for the VSOP87
code online, but not the code itself. Do you have a link for the actual
code that you could send me?
-Kerry
Paul Schlyter
Kerry Shetline <shet...@cyberzone.net> wrote:
> Thanks, Paul, for the fast and very helpful response.
>
> Paul Schlyter wrote:
>> series, coded by me in C, back in 1987. I only coded the heliocentric
>> series, since the geocentric series were both longer and less accurate,
>> and one can compute the geocentric positions from the heliocentric
>> positions anyway.
>
> external to my program so that they can easily be modified or replaced
> if necessary. A section of the data file appears as follows:
>
> MOON
> 60.40974
> 1
> 1
> 22640 0 S 2,1
> -4586 0 S 2,1 4,-2
> 2370 0 S 4,2
> 769 0 S 2,2
> -668 0 S 8,1
> -412 0 S 3,2
> -212 0 S 2,2 4,-2
> -206 0 S 2,1 4,-2 8,1
available, and one therefore can waste a lot of CPU cycles and still
get a decent performance. But back in 1983, when I first implemented
this (in FORTRAN btw; it got ported to Pascal in 1986 and to C in
1987), my target machine was an 8-bit micro: the Apple II running
CP/M-80 on an add-on board having a 2 MHz Z80. Floating-point
support was done in software, and to get a decent performance, I used
single precision as much as I could. For performance reasons I put
the series in the code, and I also wrote the code such that it called
as few trig functions as possible. I would probably not write the
code that way today -- but since the code already is available, I
still use it when I want lighing fast performance, and still decent
accuracy (ca 10 arcsec).
Using that code I have a program which computes the
rise/set/culmination times for the Sun and Moon daily, and the
planets from Mercury to Uranus weekly, and outputs a "calendar" with
these rise/set times. This program has survived through several
incarnations, and I still use it. Back in 1983, it took one full
hour for this program to compute one year of rise/set times of only
the Sun and Moon (no planets) on a 2 MHz Z80, and then I did as much
as I could in single precision. Today, the same program runs in
about one second: then it uses double precision everywhere, and it
also computes rise/set times for six planets once a week. A quite
dramatic increase in execution speed!
If I was to write that code today, I might try a different approach:
instead of reading and interpreting that external table at run time,
as you intend to do, I would want to write a "compiler" for the
table, which reads the series from the table and then outputs source
code for a function, in C or in Java, with the series directly in the
code. Then that function is compiled into the final program. This
would give both flexibility and speed - and it would probably be a
fun project too!
>> I ended up extensively using only the series for
>> Jupiter and Saturn; for the other celestial bodies I picked the series
>> by Jean Meeus in his "Astronomical formulae for Calculators" instead,
>> which I found to be slightly more accurate.
>
> probably much shorter and faster to compute. If on top of that the
> formulae are actually more accurate, that will be quite a bonus.
Flandern - Pulkkinen formulae, but Meeus followed another route:
instead of deriving his own new series, he instead used already
existing series (those by Simon Newcomb, some 100 years ago, which
were used in the Astronomical Ephemeris and Astronomical Almanac
until 1983), but truncated them to a more manageable size.
BTW the title "Astronomical Formulae for Calculators" is misleading,
since no pocket calculator available at the time of publishing of
this book had the capacity to implement the planetary position series
included in this book. With some effort, I managed to squeeze a
truncated version of the Van Flandern - Pulkkinen series (the same
one I have on my web site btw) into the HP-41C pocket calculator,
back in 1980. It occupied most of that calculator's available memory
(with all extra RAM modules inserted), or some 1.5 kBytes. Those
were the days....
Or, possible, when Meeus named this earlier book "Astronomical
Formulae for Calculators", he might not have referred to electronic
pocket calculators, but to human calculators, i.e. people like
himself!
> By the way, I just checked amazon.com for the book. It's a special order
> and will take 4-6 weeks, unfortunately. Guess it doesn't rank in sales
> with the latest Stephen King :) By the way, there's a newer book by Jean
> Meeus, "Astronomical Algorithms", that came out in 1998 -- that looks
> like the way to go now.
reasons. Everything in the older boook is available in the newer
book as well, but in an updated form. Some 5 years ago, the earlier
book was on sale at www.skypub.com, and was available immediately.
>> Then you may need to learn some FORTRAN, since there are implementations
>> in FORTRAN of these series, to their full precision, on the Net. All
>> the ASCII source code is some 35 MBytes in size....
>
> wade through it. I think I've found what may be the APIs for the VSOP87
> code online, but not the code itself. Do you have a link for the actual
Check the "How to compute planetary positions" link at my web page:
if you go to that link, and then to the bottom of that page, you'll
find several other links, and among those one link to the full
VSOP87 and ELP2000 theories.
Chris Marriott
Paul Schlyter wrote in message <86vipp$8ph$1...@merope.saaf.se>...
>Check the "How to compute planetary positions" link at my web page:
>if you go to that link, and then to the bottom of that page, you'll
>find several other links, and among those one link to the full
>VSOP87 and ELP2000 theories.
If you want ultra-high-precision planetary positions, the "Planetary Series
(1996)" theory is much, much more accurate than VSOP87, although only for
the restricted date range 1900-2100. VSOP87 has errors, compared with the
JPL numerical ephemerides, of a large fraction of an arc second, whereas
PS96 has errors of the order of milli-arcseconds.
The only programs I'm aware of which fully implement both VSOP87 and PS96 to
their full precision are my own SkyMap, and Bill Gray's "Guide", although
there may be others. That, I suspect, is why both SkyMap and Guide do so
well in reproducing the ultra-critical test of the occultation of 47 Sgr
(??) by Titan a few years ago which you yourself observed!
Chris
-----------------------------------------------------------------------
Chris Marriott, SkyMap Software, UK (ch...@skymap.com)
Visit our web site at http://www.skymap.com
Astronomy software written by astronomers, for astronomers
Paul Schlyter
Chris Marriott <ch...@NOSPAM.skymap.com> wrote:
> Paul Schlyter wrote in message <86vipp$8ph$1...@merope.saaf.se>...
>> if you go to that link, and then to the bottom of that page, you'll
>> find several other links, and among those one link to the full
>> VSOP87 and ELP2000 theories.
>
> (1996)" theory is much, much more accurate than VSOP87, although only for
> the restricted date range 1900-2100. VSOP87 has errors, compared with the
> JPL numerical ephemerides, of a large fraction of an arc second, whereas
> PS96 has errors of the order of milli-arcseconds.
by planets of other celestial objects. Or if you're JPL and want to
send spacecrafts to these planets. For other applications the VSOP87
accuracy is sufficient, and its larger valid time span is definitely
an advantage.
Also: remember that he wanted to do this in Java, a pseudo-
interpreted language with a performance penalty because of the
interpretation phase. You could use most accurate, and longest,
series in C, C++ or Fortran, but I think you should think twice
before doing it in pseudo-interpreted languages, such a Java or UCSD
Pascal (which too implemented the idea of a standard virtual machine,
some 20 years before Java did it!).
Of course, even the full VSOP87 is probably too much for a Java
implementation, but one can truncate the VSOP87 series at a suitable
point.
> The only programs I'm aware of which fully implement both VSOP87 and PS96 to
> their full precision are my own SkyMap, and Bill Gray's "Guide", although
> there may be others. That, I suspect, is why both SkyMap and Guide do so
> well in reproducing the ultra-critical test of the occultation of 47 Sgr
> (??)
28 Sgr
> by Titan a few years ago
11 years ago - in 1989
> which you yourself observed!
Yes! It was a very memorable event, perhaps even more exciting than
a total solar eclipse; much because of the rarity of such a phenomenon.
However, when observing such an event, you *also* need a good theory
for the motion of the satellite.
Chris Marriott
Paul Schlyter wrote in message <871e0c$p4c$1...@merope.saaf.se>...
>In article <949183275.12679.0...@news.demon.co.uk>,
>Chris Marriott <ch...@NOSPAM.skymap.com> wrote:
>> If you want ultra-high-precision planetary positions, the "Planetary
Series
>> (1996)" theory is much, much more accurate than VSOP87, although only for
>> the restricted date range 1900-2100. VSOP87 has errors, compared with the
>> JPL numerical ephemerides, of a large fraction of an arc second, whereas
>> PS96 has errors of the order of milli-arcseconds.
>
>This extra accuracy would be significant only for computing occultations
>by planets of other celestial objects. Or if you're JPL and want to
>send spacecrafts to these planets. For other applications the VSOP87
>accuracy is sufficient, and its larger valid time span is definitely
>an advantage.
True. One other benefit of PS96, however, is that it has considerably fewer
"terms" than VSOP87 and is much easier and quicker to compute. If you want
"current" planetary positions, it really is a much better choice, IMHO.
>> The only programs I'm aware of which fully implement both VSOP87 and PS96
to
>> their full precision are my own SkyMap, and Bill Gray's "Guide", although
>> there may be others. That, I suspect, is why both SkyMap and Guide do so
>> well in reproducing the ultra-critical test of the occultation of 47 Sgr
>> (??)
>
>28 Sgr
>
>> by Titan a few years ago
>
>11 years ago - in 1989
>
>> which you yourself observed!
>
>Yes! It was a very memorable event, perhaps even more exciting than
>a total solar eclipse; much because of the rarity of such a phenomenon.
>
>However, when observing such an event, you *also* need a good theory
>for the motion of the satellite.
You also need a third thing - the ability to correctly compute the apparent
position of the star. It's surprising how few programs do this
"rigourously"!
Joe Knapp
> >> If you want ultra-high-precision planetary positions, the "Planetary
> Series
> >> (1996)" theory is much, much more accurate than VSOP87, although only for
> >> the restricted date range 1900-2100. VSOP87 has errors, compared with the
> >> JPL numerical ephemerides, of a large fraction of an arc second, whereas
> >> PS96 has errors of the order of milli-arcseconds.
One alternative for somewhat lower precision would be Bretagnon and Simon's
"Planetary Tables & Programs -4000 to 2800." For example, they quote an accuracy
for the Mars position of ~24 arcsec in current times and no worse than ~36
arcsec over the entire date range.
http://www.amazon.com/exec/obidos/ASIN/0943396085/o/qid=949456440/sr=8-2/104-6135283-0402060
Joe
Paul Schlyter
Chris Marriott <ch...@NOSPAM.skymap.com> wrote:
> Paul Schlyter wrote in message <871e0c$p4c$1...@merope.saaf.se>...
>>In article <949183275.12679.0...@news.demon.co.uk>,
>>Chris Marriott <ch...@NOSPAM.skymap.com> wrote:
> Series
>>> (1996)" theory is much, much more accurate than VSOP87, although only for
>>> the restricted date range 1900-2100. VSOP87 has errors, compared with the
>>> JPL numerical ephemerides, of a large fraction of an arc second, whereas
>>> PS96 has errors of the order of milli-arcseconds.
>>
>>by planets of other celestial objects. Or if you're JPL and want to
>>send spacecrafts to these planets. For other applications the VSOP87
>>accuracy is sufficient, and its larger valid time span is definitely
>>an advantage.
>
> True. One other benefit of PS96, however, is that it has considerably fewer
> "terms" than VSOP87 and is much easier and quicker to compute. If you want
> "current" planetary positions, it really is a much better choice, IMHO.
use a truncated version of VSOP87 - for instance the version which was
included in Jean Meeus' "Astronomical Algorithms". That will still give
you an accuracy to (a farily large) fraction of an arcsecond, which should
be enough in almost all cases.
>>> The only programs I'm aware of which fully implement both VSOP87 and PS96
> to
>>> their full precision are my own SkyMap, and Bill Gray's "Guide", although
>>> there may be others. That, I suspect, is why both SkyMap and Guide do so
>>> well in reproducing the ultra-critical test of the occultation of 47 Sgr
>>> (??)
>>
>>28 Sgr
>>
>>> by Titan a few years ago
>>
>>11 years ago - in 1989
>>
>>> which you yourself observed!
>>
>>Yes! It was a very memorable event, perhaps even more exciting than
>>a total solar eclipse; much because of the rarity of such a phenomenon.
>>
>>However, when observing such an event, you *also* need a good theory
>>for the motion of the satellite.
>
> You also need a third thing - the ability to correctly compute the
> apparent position of the star.
that accuracy. Before Hipparcos, such a catalog didn't really exist.
> It's surprising how few programs do this "rigourously"!
by satellites of other planets. The only ones I know about so far is
Bill Gray's "Guide" and your "SkyMap". If all you want to do is to
plot a star background, you can certainly do quite well without
doing these computations rigourously.
Paul Schlyter
Joe Knapp <j...@worldnet.att.net> wrote:
>>>> If you want ultra-high-precision planetary positions, the "Planetary
>>>> Series
>>>> (1996)" theory is much, much more accurate than VSOP87, although only for
>>>> the restricted date range 1900-2100. VSOP87 has errors, compared with the
>>>> JPL numerical ephemerides, of a large fraction of an arc second, whereas
>>>> PS96 has errors of the order of milli-arcseconds.
>
> "Planetary Tables & Programs -4000 to 2800." For example, they quote an accuracy
> for the Mars position of ~24 arcsec in current times and no worse than ~36
> arcsec over the entire date range.
in Meeus' earlier book "Astronomical Formulae for Calculators" or Van
Flandern and Pulkkinen "Low Precision Planetary Positions" paper
(from which I wrote that C code I posted here a few days ago): both
yield typical errors of some 10 arc seconds. For epochs very far
into the past or future, Bretagnon-Simon's book may give better
accuracy though.
I didn't like the Bretagnon-Simon book for another reason: all the
terms of the series were give in the form:
V * (sin or cos) ( A * t + B )
where t=time and A and B just were some "magical numbers", different
for each term, with no obvious connection to one another. This makes it
quite hard to try to rewrite the series in order to not have to compute
as many trigs. I definitely prefer the variety where a few fundamental
arguments, A, B, C, D, are computed, and then each term is given as:
V * (sin or cos) ( i*A + j*B + k*C + m*D )
where i,j,k,m all are fairly small integer values, positive or negative,
sometimes zero.
> http://www.amazon.com/exec/obidos/ASIN/0943396085/o/qid=949456440/sr=8-2/104-6135283-0402060
Also available from http://www.skypub.com and http://www.willbell.com
Joe Knapp
Paul Schlyter <pau...@merope.saaf.se> wrote in message
news:878nad$ct$1...@merope.saaf.se...
Yep, magic number galore! Here's an example java implementation of the Mars
calc:
// return heliocentric coordinates for Mars at the current date/time
// input: u -- current date in unit of 10,000 Julian days from J2000
// returns: spherical and rectangular coordinate vectors:
// spherCoords -- lng (radians), lat (radians), radius (AU)
// rectCoords -- X, Y, Z -- all in AU
//
// For algorithm see Bretagnon and Simon, "Planetary Programs and Tables"
public void marsHelio(double[] spherCoords,double[] rectCoords,double u)
{
double mlng, mlat, mrad ; // spherical coords lng,lat,rad
double[] marL1 = {186563.7,18135.0,-1332.0,-704.0,-65.0,-89.0};
double[] marL2 =
{0.337967,33405.348759,0.031676,-0.007354,0.001143,-0.00029,-0.0001};
double[] marB1 =
{319714.0,-10277.0,24272.0,-2420.0,-10850.0,3880.0,5310.0,-1050.0};
double[] marB2 =
{5.339102,33407.21879,0.048,-0.04831,0.01402,0.029,-0.0073,-0.0112};
double[] marB3 = {29803.0,1904.0,1865.0,-60.0,-950.0,220.0,270.0};
double[] marB4 = {5.67694,66812.5668,0.0803,-0.0536,0.0147,0.028};
double[] marB5 = {3137.0,472.0,111.0,70.0};
double[] marB6 = {6.0173,100217.928,0.093,-0.086,0.037};
double[] marR1 = {141849.5,13651.8,-1230.0,-378.0,187.0,-153.0,-73.0};
double[] marR2 =
{3.479698,33405.34956,0.030669,-0.00909,0.00223,0.00083,-0.00048};
double[] marR3 = {6607.8,1272.8,-53.0,-46.0,14.0,-12.0,99.0};
double[] marR4 =
{3.81781,66810.6991,0.0613,-0.0182,0.0044,0.0012,0.002};
double[][] marL = {
{424067.0,4.725053, 1.599646},
{117053.0,0.92177 , 66810.71641 },
{ 31286.0,5.17451 , 66811.55202 },
{ 10248.0,1.3885 ,100216.0920 },
{ 4933.0,1.4302 , 66813.4604 },
{ 4130.0,5.4976 ,100216.8648 },
{ 2605.0,1.382 , -33.896 },
{ 1334.0,3.104 , -34.791 },
{ 1232.0,2.420 , 28109.218 },
{ 1180.0,2.148 ,133621.421 },
{ 959.0,4.404 , 22813.140 },
{ 827.0,6.053 ,133621.745 },
{ 778.0,4.910 , 56218.431 },
{ 692.0,1.881 ,100218.894 },
{ 667.0,1.267 , -3980.684 },
{ 416.0,1.799 , 29424.634 },
{ 354.0,3.146 , 25443.942 },
{ 341.0,1.005 , 66815.503 },
{ 339.0,1.420 , 33371.351 },
{ 310.0,0.916 ,-33439.371 },
{ 294.0,2.57 , 1915.95 },
{ 260.0,2.69 , -7961.45 },
{ 245.0,3.70 , 33370.81 },
{ 242.0,2.94 , 5.41 },
{ 198.0,2.78 , 5296.69 },
{ 182.0,3.09 ,-33440.09 },
{ 167.0,4.28 , 17516.29 },
{ 153.0,2.72 , 61514.60 },
{ 148.0,5.27 , 22812.71 },
{ 141.0,4.51 , 21463.18 },
{ 128.0,1.06 ,100220.45 },
{ 126.0,5.19 , 50921.57 },
{ 105.0,0.77 , 5296.06 },
{ 103.0,5.26 , 89623.78 },
{ 86.0,3.97 , 29140.97 },
{ 81.0,5.38 , 1559.53 },
{ 72.0,0.90 ,-37386.05 },
{ 71.0,5.23 , 10592.20 },
{ 66.0,2.06 , 31273.21 },
{ 64.0,4.49 , 84327.58 },
{ 63.0,1.43 ,167027.36 },
{ 60.0,1.77 ,133625.34 },
{ 55.0,4.02 , -7.37 },
{ 55.0,6.05 , 17482.57 },
{ 47.0,4.23 ,-11942.07 },
{ 45.0,1.23 , 1914.57 },
{ 43.0,6.24 , -3981.49 },
{ 43.0,2.12 , 62829.99 },
{ 42.0,2.07 , 2131.89 },
{ 38.0,4.90 , -7962.27 },
{ 36.0,0.02 , 17514.95 },
{ 33.0,2.37 ,167026.02 },
{ 31.0,4.65 , 25443.27 },
{ 31.0,2.57 , 66776.76 },
{ 30.0,1.19 ,-66844.78 },
{ 30.0,6.06 , 35321.38 },
{ 30.0,4.35 , 62546.28 },
{ 29.0,5.79 , 207.81 },
{ 28.0,1.17 , 13501.91 },
{ 28.0,2.09 ,-31489.48 }
};
double[][] marB = {
{ 32962.0,1.67255 , 1.77247 },
{ 9705.0,3.5531 , 0.9198 },
{ 367.0,0.353 ,133623.307 },
{ 101.0,4.74 ,133624.11 },
{ 44.0,0.94 , 33373.21 },
{ 44.0,2.83 ,-33441.23 },
{ 40.0,0.42 ,167028.67 }
};
double[][] marR = {
{ 27946.0,4.8846 , 0.4677 },
{ 5147.0,4.5968 ,100216.0535 },
{ 2196.0,2.568 ,100216.663 },
{ 811.0,5.560 , 28109.218 },
{ 749.0,1.772 , 56218.430 },
{ 559.0,5.112 ,133621.486 },
{ 503.0,1.272 , 22813.139 },
{ 332.0,2.701 ,133621.941 },
{ 258.0,4.56 , 33371.35 },
{ 248.0,4.93 , 29424.63 },
{ 236.0,0.92 ,-33439.37 },
{ 231.0,0.09 , 25443.93 },
{ 186.0,0.56 , 33370.81 },
{ 138.0,3.09 ,-33440.09 },
{ 117.0,2.11 , 50921.57 },
{ 110.0,1.27 , -3980.70 },
{ 99.0,5.84 , 61514.59 },
{ 90.0,4.41 , 5296.11 },
{ 90.0,5.64 ,167026.95 },
{ 81.0,2.10 , 10592.24 },
{ 80.0,2.83 , -7961.39 },
{ 74.0,1.50 , 21463.25 },
{ 73.0,1.25 , 84327.62 },
{ 71.0,2.86 ,167027.21 },
{ 69.0,2.10 , 22812.70 },
{ 69.0,2.13 , 89623.77 },
{ 63.0,1.29 , 17516.40 },
{ 57.0,0.83 , 29140.97 },
{ 53.0,0.90 ,-37386.04 }
};
mlng = reduce(6.2458611 +
33408.5620646 * u +
1E-6 * ps(marL1,6,u) *
sinred(ps(marL2,7,u)) +
1E-7* sinsum(marL,60,u)) ;
mlat = 1E-7*
(ps(marB1,8,u) *
sinred(ps(marB2,8,u)) +
ps(marB3,7,u) *
sinred(ps(marB4,6,u)) +
ps(marB5,4,u) *
sinred(ps(marB6,5,u)) +
sinsum(marB,7,u)) ;
mrad = 1.5298560 + 1E-6 * (ps(marR1,7,u) *
cosred(ps(marR2,7,u)) +
ps(marR3,7,u) *
cosred(ps(marR4,7,u))) +
1E-7* cossum(marR,29,u) ;
spherCoords[0] = mlng ;
spherCoords[1] = mlat ;
spherCoords[2] = mrad ;
// calc rectangular XYZ coords, all in AU
rectCoords[0] = mrad * Math.cos(mlat) * Math.cos(mlng);
rectCoords[1] = mrad * Math.cos(mlat) * Math.sin(mlng);
rectCoords[2] = mrad * Math.sin(mlat);
}
// Various math utilities used above:
// Given a 2d array 'tbl' of 'count' 3-element vectors, calculate summation
of a[0]*sin(a[1] + a[2]*u)
public double sinsum(double[][] tbl,int count, double u)
{
int i;
double result;
result = 0.0;
for(i = 0; i < count; ++i)
result += tbl[i][0] *
sinred(tbl[i][1] + tbl[i][2] * u);
return(result);
}
// calculate summation of a[0]*cos(a[1] + a[2]*u)
public double cossum(double[][] tbl,int count, double u)
{
int i;
double result;
result = 0.0;
for(i = 0; i < count; ++i)
result += tbl[i][0] *
cosred(tbl[i][1] + tbl[i][2] * u);
return(result);
}
public double pix2 = 6.28318530718 ; // 2pi
// do sin or cos of an angle after first reducing it to normal range 0-2pi
public double sinred(double angle)
{
return(Math.sin(fmod(angle,pix2)));
}
public double cosred(double angle)
{
return(Math.cos(fmod(angle,pix2)));
}
// reduce an angle to the range 0 - 2 Pi
public double reduce(double angle)
{
double temp;
temp = fmod(angle,pix2);
return((temp >= 0.0) ? temp : temp + pix2);
}
// calculate power series of an array
public double ps(double[] tbl,int count, double u)
{
int i;
double result;
result = tbl[count-1]; /* last factor */
for(i = count - 2; i >= 0; --i)
result = tbl[i] + result * u;
return(result);
}
// floating point modulo operation
public double fmod(double f,double div)
{
int idiv ;
idiv= (int)(f/div) ;
return f - idiv*div ;
}
Paul Schlyter
Joe Knapp <j...@worldnet.att.net> wrote:
> Yep, magic number galore! Here's an example java implementation of the Mars
> calc:
> // floating point modulo operation
> public double fmod(double f,double div)
> {
> int idiv ;
>
> idiv= (int)(f/div) ;
> return f - idiv*div ;
> }
// floating point modulo operation
public double fmod(double f,double div)
{
}
???
This version has two advantages:
1. There's no risk of integer overflow.
2. It works as expected also for negative values of f.
------------------------------------------------------------
Finally, one question: Why are you doing this in Java? Why not in
plain C instead? Your Java programming style appears very much C-like.
Joe Knapp
Paul Schlyter <pau...@merope.saaf.se> wrote > Why not instead:
>
> // floating point modulo operation
> public double fmod(double f,double div)
> {
> }
>
>
> This version has two advantages:
>
> 1. There's no risk of integer overflow.
>
> 2. It works as expected also for negative values of f.
>
> plain C instead? Your Java programming style appears very much C-like.
1. The program works as-is (e.g., no risk of integer overflow in this calc).
2. It's part of a web Java applet (the Mars calc ported from C code).
Joe
Joe Knapp
Paul Schlyter <pau...@merope.saaf.se> wrote
> > 1. The program works as-is (e.g., no risk of integer overflow in this
calc).
>
> do you have to NOT remove this potential bug?
Never, ever mess with working software! Also, life is too short to worry
about unpotentiated bugs...
Joe
Sam Midkiff
> <snip>
>Also: remember that he wanted to do this in Java, a pseudo-
>interpreted language with a performance penalty because of the
>interpretation phase. You could use most accurate, and longest,
>series in C, C++ or Fortran, but I think you should think twice
>before doing it in pseudo-interpreted languages, such a Java or UCSD
>Pascal (which too implemented the idea of a standard virtual machine,
>some 20 years before Java did it!).
>
>Of course, even the full VSOP87 is probably too much for a Java
>implementation, but one can truncate the VSOP87 series at a suitable
>point.
The research project I have been working on has shown that it is
possible to achieve from 60 - 100% of Fortran performance with Java.
There is not a lot that is intrinsic to the Java language that forces
bad performance on computationally intensive numeric codes. For more
details, see our project web page at
http://www.research.ibm.com/ninja/. Some of our software (the Array
package for Java) is available at
http://www.alphaworks.ibm.com/tech/ninja. Finally, an organization
devoted to Java for numeric computing is the Java Grande Forum, at
www.javagrande.org/.
Having said that, any commerically available off the shelf Java
implementation that i know of suffers from the problems you describe.
Sam Midkiff
slmidkiff connecting from att dot worldnet dot net
These are my opinions, and not the opinions of my employer.
Paul Schlyter
Sam Midkiff <mid...@watson.ibm.com> wrote:
> On 30 Jan 2000 14:23:24 +0100, Paul Schlyter wrote:
>
>> <snip>
>> Also: remember that he wanted to do this in Java, a pseudo-
>> interpreted language with a performance penalty because of the
>> interpretation phase. You could use most accurate, and longest,
>> series in C, C++ or Fortran, but I think you should think twice
>> before doing it in pseudo-interpreted languages, such a Java or UCSD
>> Pascal (which too implemented the idea of a standard virtual machine,
>> some 20 years before Java did it!).
>>
>> Of course, even the full VSOP87 is probably too much for a Java
>> implementation, but one can truncate the VSOP87 series at a suitable
>> point.
>
> The research project I have been working on has shown that it is
> possible to achieve from 60 - 100% of Fortran performance with Java.
too. I actually had one implementation for CP/M-80 many years
ago, called "Nevada Fortran", which compiled to pseudo-code which
then was interpreted.... just like Java!
> There is not a lot that is intrinsic to the Java language that forces
> bad performance on computationally intensive numeric codes.
> For more
> details, see our project web page at
> http://www.research.ibm.com/ninja/. Some of our software (the Array
> package for Java) is available at
> http://www.alphaworks.ibm.com/tech/ninja. Finally, an organization
> devoted to Java for numeric computing is the Java Grande Forum, at
> www.javagrande.org/.
>
> Having said that, any commerically available off the shelf Java
> implementation that i know of suffers from the problems you describe.
Not everyone can write their own efficient Java compilers...
Clark Lindsey
versions from Sun. "Just-in-time" compilers are now available
for most platforms. JIT compilers dynamically convert the Java
byte codes to native processor instructions. So the first pass
through, say a loop, is at the slower interpreted pace but from
then on you get the full processor speed.
Sun's new Hotspot compiler is even more sophisticated. It
monitors the program as it is running and then adaptively optimizes
those sections that are soaking up time.
If you really demand the highest speeds, then there are also
compilers that output directly to native executables. This loses the
Java virtue of portability but for some applications that may not
be important.
These advances explain why Java has gone from being 10-20
times slower than comparable C/C++ programs to just 10-20% slower.
In fact, in real life a Java program can just as easily end up
running faster than a similar C program because seldom
are programs optimized to their ultimate potential performance.
Finally, the faster development times with Java will for many people
more than make up for small performance penalties. I know this
is not just hype because I've had programmers tell me that
their development times with Java are as much as 50% faster
than comparable C++ projects they did previously.
- Clark Lindsey
clarkl...@hobbyspace.com
Paul Schlyter <pau...@merope.saaf.se> wrote in message
news:87hpi0$hkd$1...@merope.saaf.se...
Sam Midkiff
wrote:
>In article <389c1a9d...@netnews.worldnet.att.net>,
>Sam Midkiff <mid...@watson.ibm.com> wrote:
>
>> On 30 Jan 2000 14:23:24 +0100, Paul Schlyter wrote:
>>
>>
>> The research project I have been working on has shown that it is
>> possible to achieve from 60 - 100% of Fortran performance with Java.
>
>.....well, OK, there probably are some inefficient Fortran compilers
>too. I actually had one implementation for CP/M-80 many years
>ago, called "Nevada Fortran", which compiled to pseudo-code which
>then was interpreted.... just like Java!
The compiler we compared against was IBM's xlf compiler, which is one
of the better Fortran compilers (not having numbers at hand, I won't
say that it is the best.) It was the latest version, using high
optimization levels. Our kernel code was compared against hand tuned
libraries (IBM's ESSL). To use a 1980's compiler for comparison, or
a modern compiler crippled with low opt levels, would have been not
only dishonest, but bad science -- which I try not to do, as a matter
of prinicple. And, in the interest of accuracy, some of the 60-100%
numbers use the FMA instruction, which is not, strictly speaking,
legal Java (see a few paragraphs later). The papers have the
relevant details for individual benchmarks.
>
>> There is not a lot that is intrinsic to the Java language that forces
>> bad performance on computationally intensive numeric codes.
>
>Not even the interpretation of the J-code ????
> snip...
There is no such beast as J-code mentioned in the language spec. I
assume you are referring to byte code?
The real performance limiter on numeric codes is the fact that
specification for floating point arithmetic disallows the FMA
(floating multiply-add) instruction found on some machines. There is
a JSR (Java Specification Request) wending its way through the process
to fix this. A similar fix has already been made to accomodate the
way that X86 architectures handle floating point. The problem with
FMA's is that they are not only faster, but more precise, which messes
up bit-for-bit compatibility.
Nothing I see in the Java Language Specification requires
interpretation. I think the reference implementation does
interpretation, and the semantics of the language are defined in
terms of interpretation, but as long as an implemention of the
language behaves according to the spec, how it is actually implemented
is irrelevant. (More realistically, as long as you pass Sun's
verification process, it doesn't make any difference how it is
implemented.) Sun, IBM, and others sell non-interpreting Java Virtual
Machines.
>>
>>
>> Having said that, any commerically available off the shelf Java
>> implementation that i know of suffers from the problems you describe.
>
>Well, that's the kinds of Java's which are available to most of us.
>Not everyone can write their own efficient Java compilers...
>
My point, which I seem not to have made clearly enough, is that the
problems raised with Java are not intrinsic to the language, but
instead are artifacts of current implementations. Therefore, the
situation with respect to performance should improve dramatically over
time. When developing a piece of software with a significant
lifetime, consideration should be given not only to current
performance and availability of the implementation environment, but
also the expected performance over the lifetime of the project.
As well, current off-the-shelf version of Java can give 25% or so of
Fortran. Given current improvements in CPU speeds, you will have 2000
Fortran performance sometime in 2002 or so, without improvements to
the Java language processing. Because Java is a relatively immature
language, performance improvements should come much faster than for
more mature languages like C, C++ and Fortran.
Sam Midkiff
these are my opinions, and not the opinions of the IBM Corp.