Accurate and fast Sun/Moon ephemerides suitable for Android (and iOS) projects

NEW: The code in this page is also available at this respository, including all previous versions of the code with diff files. You can use it to track changes and update the code in your project. For instance, the diff file between the July 3, 2020 and the August 7, 2020 versions would be this link.
It is very common the need of ephemerides for the Sun and/or the Moon in other projects. For projects related with my institute I have developed from time to time some programs using JPARSEC, and often thought it would be better to have a more simple program completely independent from JPARSEC, which is not very handy for simple calculations. Recently I was asked about ephemerides for the Sun and the Moon for a commercial project in Android. I refused to give permission to use JPARSEC in a closed project, but suggested as an alternative the program I present in this post, very suitable for Android devices.
The class SunMoonCalculator is a good implementation of ephemerides considering such a little piece of code. The implementation takes into account precession, nutation, aberration, and also the difference between terrestrial time and universal time. Rise and set times are calculated considering the angular radius of the object and the refraction, but not the depression of the horizon (that requires knowledge of the elevation about sea level) since this is the standard way to give rise and set times. Everything is simplified but the corrections are there and they work as expected, giving accuracies around 0.001 degrees for the Sun and 0.01 degrees for the Moon, tested between 1800 and 2200. The program can calculate accurate rise, set, and transit times, measured from local horizon and also civil, nautical, and astronomical twilights. Accuracy in the times of rise and set are around <1s for the Sun and few seconds for the Moon.
To obtain the difference TTUT1 I have done a polynomial fitting of this difference in two ranges of years: from year 600 to 1600, and from 1600 to 2200. Values between 2020 and 2200 are extrapolated and only tentative. Outside this range a parabola is used following the work by Morrison et al. (2004). The effect of this correction in the output values is not very relevant (unless the error gets above the minute of time), since the main source of uncertainty is the reduced number of terms in the calculation of the ecliptic coordinates of the Moon. Anyway, the combination of accurate secular expansions by S. Moshier with a considerable number of periodic terms allows the program to reach an accuracy within the arcminute over many centuries or even several millenia, which is a remarkably good accuracy. You may notice the elevation above the sea level is not an input to the program, so the topocentric positions are not corrected for this. Since this correction can improve the output for the Moon, in the section about improving the position of the Moon I present a better topocentric correction.
The main ephemerides methods are getSun() and getMoon(), that returns the ecliptic longitude and latitude, the distance, and the angular radius of both objects. A common method doCalc(…) is used to transform these values into the final coordinates azimuth and elevation, although right ascension (RA) and declination (DEC) are also calculated. Elevation is corrected for refraction (if it is above 3 degrees), but RA and DEC are not (another 'standard' convention when calculating ephemerides). Then the method obtainAccurateRiseSetTransit(…) is used to correct the rise, set, and transit times taking into account the movement of the Sun or Moon between the calculation time and the rise/set/transit time. This is required for any moving object, for stars it isn't necessary. The transit time also includes calculation of the maximum elevation, also corrected for refraction. The topocentric distance to the Sun/Moon is also provided, as well as the age of the Moon (approximate number of days since last new Moon). The distances are topocentric, you can use the geocentric flag to get the geocentric distances. For static objects I describe how to compute their positions starting from J2000 catalog coordinates at the end of the discussion (2020/09/01).
The main program contains a simple test case. The method getDateAsString(…) is used to show the date for a given rise/set/transit time. It considers the possible output of 1 in these values in case the Sun/Moon has no rise/set/transit for a given observer, for instance one close to North of South pole (or for geocentric calculations). Input and output times are in Universal Time (UT). The method calcSunAndMoon is used to compute everything at once comfortably, you may need to call only the Sun or the Moon (doCalc(getSun()) for the Sun) but accurate rise/set/transit times require iterations. For an accurate Moon age in days you should compute the Sun position before the Moon.
My tests show this program is very accurate during many centuries around year 2000. In case you would like to use it in a commercial program, please let me know, but no explicit permission is required. In case you will be using this code for some specific interesting purpose related to space, please let me know, I could help improving the accuracy or adapting it to other needs.
The program has been extended from its original version to support additional computations. The program listed below can also compute lunar phases and equinoxes/solstices. Accuracy of these results is around one minute for lunar phases and 23 minutes for the seasons. To reach this accuracy I have changed the algorithm to compute the position of the Sun to use the expansion in “Planetary Programs and Tables”, by Pierre Bretagnon and JeanLouis Simon, WillmanBell, 1986. I have also extended it to support ephemerides of planets with a very good accuracy over a millenia around year 2000. The code was also rewritten to improve its quality and readability, so it is quite different from the first original version.
/** * A very simple yet accurate Sun/Moon calculator without using JPARSEC library. * @author T. Alonso Albi  OAN (Spain), email t.alonso@oan.es * @version January 19, 2021 (added TTUT in ancient times, mean elements up to t^7 in getMoon to improve accuracy in ancient times) * @version January 11, 2021 (fixed an error with Julian day computation, refraction iteration improved) * @version November 10, 2020 (better refraction correction, still with Bennet's formula) * @version August 7, 2020 (forgot to add equation of equinoxes to lst, affects 1s at most to rise/set/transit times) * @version July 1, 2020 (improved aberration in getSun, illumination phase for planets fixed, refraction with ambient P/T and extensible to radio wavelengths) * @version June 15, 2020 (more terms for Sun, equinoxes/solstices methods revised to improve accuracy and performance) * @version June 12, 2020 (added planetary ephemerides from another class, better performance) * @version June 3, 2020 (fix rise/set problem for high latitudes, return azimuth and elevation even in geocentric computations) * @version November 26, 2018 (two new methods getCulminationTime and getAzimuthTime) * @version November 6, 2018 (better accuracy for Moon, angular radius in ephemeris, cosmetic improvements) * @version July 24, 2018 (new class to hold results, illumination phase, moon phases, equinoxes and solstices) * @version May 25, 2017 (fixed nutation correction and moon age, better accuracy in Moon) */ public class SunMoonCalculator { /** Radians to degrees. */ public static final double RAD_TO_DEG = 180.0 / Math.PI; /** Degrees to radians. */ public static final double DEG_TO_RAD = 1.0 / RAD_TO_DEG; /* Arcseconds to radians */ public static final double ARCSEC_TO_RAD = (DEG_TO_RAD / 3600.0); /** Astronomical Unit in km. As defined by JPL. */ public static final double AU = 149597870.691; /** Earth equatorial radius in km. IERS 2003 Conventions. */ public static final double EARTH_RADIUS = 6378.1366; /** Two times Pi. */ public static final double TWO_PI = 2.0 * Math.PI; /** Pi divided by two. */ public static final double PI_OVER_TWO = Math.PI / 2.0; /** Julian century conversion constant = 100 * days per year. */ public static final double JULIAN_DAYS_PER_CENTURY = 36525.0; /** Seconds in one day. */ public static final double SECONDS_PER_DAY = 86400; /** Light time in days for 1 AU. DE405 definition. */ public static final double LIGHT_TIME_DAYS_PER_AU = 0.00577551833109; /** Our default epoch. The Julian Day which represents noon on 20000101. */ public static final double J2000 = 2451545.0; /** The set of twilights to calculate (types of rise/set events). */ public static enum TWILIGHT { /** * Event ID for calculation of rising and setting times for astronomical * twilight. In this case, the calculated time will be the time when the * center of the object is at 18 degrees of geometrical elevation below the * astronomical horizon. At this time astronomical observations are possible * because the sky is dark enough. */ ASTRONOMICAL, /** * Event ID for calculation of rising and setting times for nautical * twilight. In this case, the calculated time will be the time when the * center of the object is at 12 degrees of geometric elevation below the * astronomical horizon. */ NAUTICAL, /** * Event ID for calculation of rising and setting times for civil twilight. * In this case, the calculated time will be the time when the center of the * object is at 6 degrees of geometric elevation below the astronomical * horizon. */ CIVIL, /** * The standard value of 34' for the refraction at the local horizon. */ HORIZON_34arcmin }; /** The set of events to calculate (rise/set/transit events). */ public static enum EVENT { /** Rise. */ RISE, /** Set. */ SET, /** Transit. */ TRANSIT } /** The set of phases to compute the moon phases. */ public static enum MOONPHASE { /** New Moon phase. */ NEW_MOON ("New Moon: ", 0), /** Crescent quarter phase. */ CRESCENT_QUARTER ("Crescent quarter:", 0.25), /** Full Moon phase. */ FULL_MOON ("Full Moon: ", 0.5), /** Descent quarter phase. */ DESCENT_QUARTER ("Descent quarter: ", 0.75); /** Phase name. */ public String phaseName; /** Phase value. */ public double phase; private MOONPHASE(String name, double ph) { phaseName = name; phase = ph; } } /** The set of bodies to compute ephemerides. */ public static enum BODY { MERCURY (0, 2439.7), VENUS (1, 6051.8), MARS (3, 3396.19), JUPITER (4, 71492), SATURN (5, 60268), URANUS (6, 25559), NEPTUNE (7, 24764), Moon (2, 1737.4), Sun (1, 696000), EMB (2, 0); /** Equatorial radius in km. */ public double eqRadius; /** Body index for computing the position. */ public int index; private BODY (int i, double r) { index = i; eqRadius = r; } } /** Input values and nutation/obliquity parameters only calculated once. */ private double obsLon = 0, obsLat = 0; private TWILIGHT twilight = TWILIGHT.HORIZON_34arcmin; protected double jd_UT = 0; protected double t = 0; protected double nutLon = 0, nutObl = 0; protected double meanObliquity = 0; protected double TTminusUT = 0; protected double lst = 0; /** * Class to hold the results of ephemerides. * @author T. Alonso Albi  OAN (Spain) */ public class Ephemeris { private Ephemeris(double azi, double alt, double rise2, double set2, double transit2, double transit_alt, double ra, double dec, double dist, double eclLon, double eclLat, double angR) { azimuth = azi; elevation = alt; rise = rise2; set = set2; transit = transit2; transitElevation = transit_alt; rightAscension = ra; declination = dec; distance = dist; illuminationPhase = 100; eclipticLongitude = eclLon; eclipticLatitude = eclLat; angularRadius = angR; } /** Values for azimuth, elevation, rise, set, and transit for the Sun. Angles in radians, rise ... * as Julian days in UT. Distance in AU. */ public double azimuth, elevation, rise, set, transit, transitElevation, distance, rightAscension, declination, illuminationPhase, eclipticLongitude, eclipticLatitude, angularRadius; } /** Ephemeris for the Sun and Moon bodies. */ public Ephemeris sun, moon; /** Moon's age in days as an independent variable. */ public double moonAge; /** * Main constructor for Sun/Moon calculations. Time should be given in * Universal Time (UT), observer angles in radians. * @param year The year. * @param month The month. * @param day The day. * @param h The hour. * @param m Minute. * @param s Second. * @param obsLon Longitude for the observer. * @param obsLat Latitude for the observer. * @throws Exception If the date does not exists. */ public SunMoonCalculator(int year, int month, int day, int h, int m, int s, double obsLon, double obsLat) throws Exception { double jd = toJulianDay(year, month, day, h, m, s); double x = year + (month  1 + day / 30.0) / 12.0; double x2 = x * x, x3 = x2 * x, x4 = x3 * x; if (year < 600  year >= 2200) { double u = (x  1820) * 0.01; TTminusUT = 20 + 32 * u * u; } else { if (year < 1600) { TTminusUT = 10535.328003  9.9952386275 * x + 0.00306730763 * x2  7.7634069836E6 * x3 + 3.1331045394E9 * x4 + 8.2255308544E12 * x2 * x3  7.4861647156E15 * x4 * x2 + 1.936246155E18 * x4 * x3  8.4892249378E23 * x4 * x4; } else { TTminusUT = 1027175.34776 + 2523.2566254 * x  1.8856868491 * x2 + 5.8692462279E5 * x3 + 3.3379295816E7 * x4 + 1.7758961671E10 * x2 * x3  2.7889902806E13 * x2 * x4 + 1.0224295822E16 * x3 * x4  1.2528102371E20 * x4 * x4; } } double c = 0.000012932 * Math.pow(x  1955.0, 2); if (year < 1955  year > 2005) TTminusUT += c; this.obsLon = obsLon; this.obsLat = obsLat; setUTDate(jd); } /** * Sets the rise/set times to return. Default is for the local horizon. * @param t The Twilight. */ public void setTwilight(TWILIGHT t) { this.twilight = t; } /** * Sets the UT date from the provided Julian day and computes the nutation, obliquity, and * sidereal time. TT minuts UT1 is not updated since it changes very slowly with time. * Use this only to update the computation time for the same year as the one used when the * instance was created. * @param jd The new Julian day in UT. */ protected void setUTDate(double jd) { this.jd_UT = jd; this.t = (jd + TTminusUT / SECONDS_PER_DAY  J2000) / JULIAN_DAYS_PER_CENTURY; // Compute nutation double M1 = (124.90  1934.134 * t + 0.002063 * t * t) * DEG_TO_RAD; double M2 = (201.11 + 72001.5377 * t + 0.00057 * t * t) * DEG_TO_RAD; nutLon = (0.0047785 * Math.sin(M1)  0.0003667 * Math.sin(M2)) * DEG_TO_RAD; nutObl = (0.002558 * Math.cos(M1)  0.00015339 * Math.cos(M2)) * DEG_TO_RAD; // Compute mean obliquity double t2 = this.t / 100.0; double tmp = t2 * (27.87 + t2 * (5.79 + t2 * 2.45)); tmp = t2 * (249.67 + t2 * (39.05 + t2 * (7.12 + tmp))); tmp = t2 * (1.55 + t2 * (1999.25 + t2 * (51.38 + tmp))); tmp = (t2 * (4680.93 + tmp)) / 3600.0; meanObliquity = (23.4392911111111 + tmp) * DEG_TO_RAD; // Obtain local apparent sidereal time double jd0 = Math.floor(jd_UT  0.5) + 0.5; double T0 = (jd0  J2000) / JULIAN_DAYS_PER_CENTURY; double secs = (jd_UT  jd0) * SECONDS_PER_DAY; double gmst = (((((6.2e6 * T0) + 9.3104e2) * T0) + 8640184.812866) * T0) + 24110.54841; double msday = 1.0 + (((((1.86e5 * T0) + 0.186208) * T0) + 8640184.812866) / (SECONDS_PER_DAY * JULIAN_DAYS_PER_CENTURY)); gmst = (gmst + msday * secs) * (15.0 / 3600.0) * DEG_TO_RAD; lst = normalizeRadians(gmst + obsLon + nutLon * Math.cos(meanObliquity + nutObl)); } /** Calculates everything for the Sun and the Moon. */ public void calcSunAndMoon() { double jd = this.jd_UT; // First the Sun sun = doCalc(getSun()); int niter = 15; // Number of iterations to get accurate rise/set/transit times sun.rise = obtainAccurateRiseSetTransit(sun.rise, EVENT.RISE, niter, true); sun.set = obtainAccurateRiseSetTransit(sun.set, EVENT.SET, niter, true); sun.transit = obtainAccurateRiseSetTransit(sun.transit, EVENT.TRANSIT, niter, true); if (sun.transit == 1) { sun.transitElevation = 0; } else { // Update Sun's maximum elevation setUTDate(sun.transit); sun.transitElevation = doCalc(getSun()).transitElevation; } // Now Moon setUTDate(jd); moon = doCalc(getMoon()); double ma = moonAge; //niter = 15; // Number of iterations to get accurate rise/set/transit times moon.rise = obtainAccurateRiseSetTransit(moon.rise, EVENT.RISE, niter, false); moon.set = obtainAccurateRiseSetTransit(moon.set, EVENT.SET, niter, false); moon.transit = obtainAccurateRiseSetTransit(moon.transit, EVENT.TRANSIT, niter, false); if (moon.transit == 1) { moon.transitElevation = 0; } else { // Update Moon's maximum elevation setUTDate(moon.transit); moon.transitElevation = doCalc(getMoon()).transitElevation; } setUTDate(jd); moonAge = ma; // Compute illumination phase percentage for the Moon getIlluminationPhase(moon); } // Sun data from the expansion "Planetary Programs // and Tables" by Pierre Bretagnon and JeanLouis // Simon, WillmanBell, 1986 private static double sun_elements[][] = { new double[] { 403406.0, 0.0, 4.721964, 1.621043 }, new double[] { 195207.0, 97597.0, 5.937458, 62830.348067 }, new double[] { 119433.0, 59715.0, 1.115589, 62830.821524 }, new double[] { 112392.0, 56188.0, 5.781616, 62829.634302 }, new double[] { 3891.0, 1556.0, 5.5474, 125660.5691 }, new double[] { 2819.0, 1126.0, 1.512, 125660.9845 }, new double[] { 1721.0, 861.0, 4.1897, 62832.4766 }, new double[] { 0.0, 941.0, 1.163, .813 }, new double[] { 660.0, 264.0, 5.415, 125659.31 }, new double[] { 350.0, 163.0, 4.315, 57533.85 }, new double[] { 334.0, 0.0, 4.553, 33.931 }, new double[] { 314.0, 309.0, 5.198, 777137.715 }, new double[] { 268.0, 158.0, 5.989, 78604.191 }, new double[] { 242.0, 0.0, 2.911, 5.412 }, new double[] { 234.0, 54.0, 1.423, 39302.098 }, new double[] { 158.0, 0.0, .061, 34.861 }, new double[] { 132.0, 93.0, 2.317, 115067.698 }, new double[] { 129.0, 20.0, 3.193, 15774.337 }, new double[] { 114.0, 0.0, 2.828, 5296.67 }, new double[] { 99.0, 47.0, .52, 58849.27 }, new double[] { 93.0, 0.0, 4.65, 5296.11 }, new double[] { 86.0, 0.0, 4.35, 3980.7 }, new double[] { 78.0, 33.0, 2.75, 52237.69 }, new double[] { 72.0, 32.0, 4.5, 55076.47 }, new double[] { 68.0, 0.0, 3.23, 261.08 }, new double[] { 64.0, 10.0, 1.22, 15773.85 }, new double[] { 46.0, 16.0, .14, 188491.03 }, new double[] { 38.0, 0.0, 3.44, 7756.55 }, new double[] { 37.0, 0.0, 4.37, 264.89 }, new double[] { 32.0, 24.0, 1.14, 117906.27 }, new double[] { 29.0, 13.0, 2.84, 55075.75 }, new double[] { 28.0, 0.0, 5.96, 7961.39 }, new double[] { 27.0, 9.0, 5.09, 188489.81 }, new double[] { 27.0, 0.0, 1.72, 2132.19 }, new double[] { 25.0, 17.0, 2.56, 109771.03 }, new double[] { 24.0, 11.0, 1.92, 54868.56 }, new double[] { 21.0, 0.0, .09, 25443.93 }, new double[] { 21.0, 31.0, 5.98, 55731.43 }, new double[] { 20.0, 10.0, 4.03, 60697.74 }, new double[] { 18.0, 0.0, 4.27, 2132.79 }, new double[] { 17.0, 12.0, .79, 109771.63 }, new double[] { 14.0, 0.0, 4.24, 7752.82 }, new double[] { 13.0, 5.0, 2.01, 188491.91 }, new double[] { 13.0, 0.0, 2.65, 207.81 }, new double[] { 13.0, 0.0, 4.98, 29424.63 }, new double[] { 12.0, 0.0, .93, 7.99 }, new double[] { 10.0, 0.0, 2.21, 46941.14 }, new double[] { 10.0, 0.0, 3.59, 68.29 }, new double[] { 10.0, 0.0, 1.5, 21463.25 }, new double[] { 10.0, 9.0, 2.55, 157208.4 } }; protected double[] getSun() { double L = 0.0, R = 0.0, t2 = t * 0.01; double Lp = 0.0, deltat = 0.5, t2p = (t + deltat / JULIAN_DAYS_PER_CENTURY) * 0.01; for (int i = 0; i < sun_elements.length; i++) { double v = sun_elements[i][2] + sun_elements[i][3] * t2; double u = normalizeRadians(v); L = L + sun_elements[i][0] * Math.sin(u); R = R + sun_elements[i][1] * Math.cos(u); double vp = sun_elements[i][2] + sun_elements[i][3] * t2p; double up = normalizeRadians(vp); Lp = Lp + sun_elements[i][0] * Math.sin(up); } double lon = normalizeRadians(4.9353929 + normalizeRadians(62833.196168 * t2) + L / 10000000.0) * RAD_TO_DEG; double sdistance = 1.0001026 + R / 10000000.0; // Now subtract aberration double dlon = ((Lp  L) / 10000000.0 + 62833.196168 * (t2p  t2)) / deltat; double aberration = dlon * sdistance * LIGHT_TIME_DAYS_PER_AU; lon = aberration * RAD_TO_DEG; double slongitude = lon * DEG_TO_RAD; // apparent longitude (error<0.001 deg) double slatitude = 0; // Sun's ecliptic latitude is always negligible return new double[] {slongitude, slatitude, sdistance, Math.atan(BODY.Sun.eqRadius / (AU * sdistance))}; } protected double[] getMoon() { // These expansions up to t^7 for the mean elements are taken from S. L. Moshier, see program cmoon /* Mean elongation of moon = D */ double x = (1.6029616009939659e+09 * t + 1.0722612202445078e+06); x += (((((3.207663637426e013 * t + 2.555243317839e011) * t + 2.560078201452e009) * t  3.702060118571e005) * t + 6.9492746836058421e03) * t /* D, t^3 */  6.7352202374457519e+00) * t * t; /* D, t^2 */ double phase = ARCSEC_TO_RAD * x; /* Mean distance of moon from its ascending node = F */ x = (1.7395272628437717e+09 * t + 3.3577951412884740e+05); x += (((((4.474984866301e013 * t + 4.189032191814e011) * t  2.790392351314e009) * t  2.165750777942e006) * t  7.5311878482337989e04) * t /* F, t^3 */  1.3117809789650071e+01) * t * t; /* F, t^2 */ double node = ARCSEC_TO_RAD * x; /* Mean anomaly of sun = l' (J. Laskar) */ x = (1.2959658102304320e+08 * t + 1.2871027407441526e+06); x += ((((((((1.62e20 * t  1.0390e17) * t  3.83508e15) * t + 4.237343e13) * t + 8.8555011e11) * t  4.77258489e8) * t  1.1297037031e5) * t + 8.7473717367324703e05) * t  5.5281306421783094e01) * t * t; double sanomaly = ARCSEC_TO_RAD * x; /* Mean anomaly of moon = l */ x = (1.7179159228846793e+09 * t + 4.8586817465825332e+05); x += (((((1.755312760154e012 * t + 3.452144225877e011) * t  2.506365935364e008) * t  2.536291235258e004) * t + 5.2099641302735818e02) * t /* l, t^3 */ + 3.1501359071894147e+01) * t * t; /* l, t^2 */ double anomaly = ARCSEC_TO_RAD * x; /* Mean longitude of moon, re mean ecliptic and equinox of date = L */ x = (1.7325643720442266e+09 * t + 7.8593980921052420e+05); x += (((((7.200592540556e014 * t + 2.235210987108e010) * t  1.024222633731e008) * t  6.073960534117e005) * t + 6.9017248528380490e03) * t /* L, t^3 */  5.6550460027471399e+00) * t * t; /* L, t^2 */ double l = ARCSEC_TO_RAD * x * RAD_TO_DEG; // Now longitude, with the three main correcting terms of evection, // variation, and equation of year, plus other terms (error<0.01 deg) // P. Duffet's MOON program taken as reference for the periodic terms double E = 1.0  (.002495 + 7.52E06 * (t + 1.0)) * (t + 1.0); l += 6.28875 * Math.sin(anomaly) + 1.274018 * Math.sin(2 * phase  anomaly) + .658309 * Math.sin(2 * phase); l += 0.213616 * Math.sin(2 * anomaly)  E * .185596 * Math.sin(sanomaly)  0.114336 * Math.sin(2 * node); l += .058793 * Math.sin(2 * phase  2 * anomaly) + .057212 * E * Math.sin(2 * phase  anomaly  sanomaly) + .05332 * Math.sin(2 * phase + anomaly); l += .045874 * E * Math.sin(2 * phase  sanomaly) + .041024 * E * Math.sin(anomaly  sanomaly)  .034718 * Math.sin(phase)  E * .030465 * Math.sin(sanomaly + anomaly); l += .015326 * Math.sin(2 * (phase  node))  .012528 * Math.sin(2 * node + anomaly)  .01098 * Math.sin(2 * node  anomaly) + .010674 * Math.sin(4 * phase  anomaly); l += .010034 * Math.sin(3 * anomaly) + .008548 * Math.sin(4 * phase  2 * anomaly); l += E * .00791 * Math.sin(sanomaly  anomaly + 2 * phase)  E * .006783 * Math.sin(2 * phase + sanomaly) + .005162 * Math.sin(anomaly  phase) + E * .005 * Math.sin(sanomaly + phase); l += .003862 * Math.sin(4 * phase) + E * .004049 * Math.sin(anomaly  sanomaly + 2 * phase) + .003996 * Math.sin(2 * (anomaly + phase)) + .003665 * Math.sin(2 * phase  3 * anomaly); l += E * 2.695E3 * Math.sin(2 * anomaly  sanomaly) + 2.602E3 * Math.sin(anomaly  2*(node+phase)); l += E * 2.396E3 * Math.sin(2*(phase  anomaly)  sanomaly)  2.349E3 * Math.sin(anomaly+phase); l += E * E * 2.249E3 * Math.sin(2*(phasesanomaly))  E * 2.125E3 * Math.sin(2*anomaly+sanomaly); l += E * E * 2.079E3 * Math.sin(2*sanomaly) + E * E * 2.059E3 * Math.sin(2*(phasesanomaly)anomaly); l += 1.773E3 * Math.sin(anomaly+2*(phasenode))  1.595E3 * Math.sin(2*(node+phase)); l += E * 1.22E3 * Math.sin(4*phasesanomalyanomaly)  1.11E3 * Math.sin(2*(anomaly+node)); double longitude = l * DEG_TO_RAD; double Psin = 29.530588853; if (sun != null) { // Get Moon age, more accurate than 'phase' but we need the Sun position moonAge = normalizeRadians(longitude  sun.eclipticLongitude) * Psin / TWO_PI; } else { // Use the phase variable as estimate, less accurate but this is used only when we don't need an accurate value moonAge = phase * Psin / TWO_PI; } // Now Moon parallax double parallax = .950724 + .051818 * Math.cos(anomaly) + .009531 * Math.cos(2 * phase  anomaly); parallax += .007843 * Math.cos(2 * phase) + .002824 * Math.cos(2 * anomaly); parallax += 0.000857 * Math.cos(2 * phase + anomaly) + E * .000533 * Math.cos(2 * phase  sanomaly); parallax += E * .000401 * Math.cos(2 * phase  anomaly  sanomaly) + E * .00032 * Math.cos(anomaly  sanomaly)  .000271 * Math.cos(phase); parallax += E * .000264 * Math.cos(sanomaly + anomaly)  .000198 * Math.cos(2 * node  anomaly); parallax += 1.73E4 * Math.cos(3 * anomaly) + 1.67E4 * Math.cos(4*phaseanomaly); // So Moon distance in Earth radii is, more or less, double distance = 1.0 / Math.sin(parallax * DEG_TO_RAD); // Ecliptic latitude with nodal phase (error<0.01 deg) l = 5.128189 * Math.sin(node) + 0.280606 * Math.sin(node + anomaly) + 0.277693 * Math.sin(anomaly  node); l += .173238 * Math.sin(2 * phase  node) + .055413 * Math.sin(2 * phase + node  anomaly); l += .046272 * Math.sin(2 * phase  node  anomaly) + .032573 * Math.sin(2 * phase + node); l += .017198 * Math.sin(2 * anomaly + node) + .009267 * Math.sin(2 * phase + anomaly  node); l += .008823 * Math.sin(2 * anomaly  node) + E * .008247 * Math.sin(2 * phase  sanomaly  node) + .004323 * Math.sin(2 * (phase  anomaly)  node); l += .0042 * Math.sin(2 * phase + node + anomaly) + E * .003372 * Math.sin(node  sanomaly  2 * phase); l += E * 2.472E3 * Math.sin(2 * phase + node  sanomaly  anomaly); l += E * 2.222E3 * Math.sin(2 * phase + node  sanomaly); l += E * 2.072E3 * Math.sin(2 * phase  node  sanomaly  anomaly); double latitude = l * DEG_TO_RAD; return new double[] {longitude, latitude, distance * EARTH_RADIUS / AU, Math.atan(BODY.Moon.eqRadius / (distance * EARTH_RADIUS))}; } /** * Compute the topocentric position of the body. * @param pos Values for the ecliptic longitude, latitude, distance and so on from previous methods for the specific body. * @return The ephemeris object with the output position */ protected Ephemeris doCalc(double[] pos) { return doCalc(pos, false); } /** * Compute the position of the body. * @param pos Values for the ecliptic longitude, latitude, distance and so on from previous methods for the specific body. * @param geocentric True to return geocentric position. Set this to false generally. * @return The ephemeris object with the output position */ protected Ephemeris doCalc(double[] pos, boolean geocentric) { // Correct for nutation in longitude and obliquity pos[0] = pos[0] + nutLon; pos[1] = pos[1] + nutObl; // Ecliptic to equatorial coordinates double cl = Math.cos(pos[1]); double x = pos[2] * Math.cos(pos[0]) * cl; double y = pos[2] * Math.sin(pos[0]) * cl; double z = pos[2] * Math.sin(pos[1]); double sinEcl = Math.sin(meanObliquity), cosEcl = Math.cos(meanObliquity); double tmp = y * cosEcl  z * sinEcl; z = y * sinEcl + z * cosEcl; y = tmp; // Obtain topocentric rectangular coordinates double sinLat = Math.sin(obsLat); double cosLat = Math.cos(obsLat); double radiusAU = geocentric ? 0 : EARTH_RADIUS / AU; double correction[] = new double[] { radiusAU * cosLat * Math.cos(lst), radiusAU * cosLat * Math.sin(lst), radiusAU * sinLat}; double xtopo = x  correction[0]; double ytopo = y  correction[1]; double ztopo = z  correction[2]; // Obtain topocentric equatorial coordinates double ra = 0.0; double dec = PI_OVER_TWO; if (ztopo < 0.0) dec = dec; if (ytopo != 0.0  xtopo != 0.0) { ra = Math.atan2(ytopo, xtopo); dec = Math.atan2(ztopo / Math.hypot(xtopo, ytopo), 1.0); } double dist = Math.sqrt(xtopo * xtopo + ytopo * ytopo + ztopo * ztopo); // Hour angle double angh = lst  ra; // Obtain azimuth and geometric alt double sinDec = Math.sin(dec), cosDec = Math.cos(dec); double h = sinLat * sinDec + cosLat * cosDec * Math.cos(angh); double alt = Math.asin(h); double azy = Math.sin(angh); double azx = Math.cos(angh) * sinLat  sinDec * cosLat / cosDec; double azi = Math.PI + Math.atan2(azy, azx); // 0 = north if (geocentric) return new Ephemeris(azi, alt, 1, 1, 1, 1, normalizeRadians(ra), dec, dist, pos[0], pos[1], pos[3]); // Get apparent elevation alt = refraction(alt); switch (twilight) { case HORIZON_34arcmin: // Rise, set, transit times, taking into account Sun/Moon angular radius (pos[3]). // The 34' factor is the standard refraction at horizon. // Removing angular radius will do calculations for the center of the disk instead // of the upper limb. tmp = (34.0 / 60.0) * DEG_TO_RAD  pos[3]; break; case CIVIL: tmp = 6 * DEG_TO_RAD; break; case NAUTICAL: tmp = 12 * DEG_TO_RAD; break; case ASTRONOMICAL: tmp = 18 * DEG_TO_RAD; break; } // Compute cosine of hour angle tmp = (Math.sin(tmp)  sinLat * sinDec) / (cosLat * cosDec); /** Length of a sidereal day in days according to IERS Conventions. */ double siderealDayLength = 1.00273781191135448; double celestialHoursToEarthTime = 1.0 / (siderealDayLength * TWO_PI); // Make calculations for the meridian double transit_time1 = celestialHoursToEarthTime * normalizeRadians(ra  lst); double transit_time2 = celestialHoursToEarthTime * (normalizeRadians(ra  lst)  TWO_PI); double transit_alt = Math.asin(sinDec * sinLat + cosDec * cosLat); transit_alt = refraction(transit_alt); // Obtain the current event in time double transit_time = transit_time1; double jdToday = Math.floor(jd_UT  0.5) + 0.5; double transitToday2 = Math.floor(jd_UT + transit_time2  0.5) + 0.5; // Obtain the transit time. Preference should be given to the closest event // in time to the current calculation time if (jdToday == transitToday2 && Math.abs(transit_time2) < Math.abs(transit_time1)) transit_time = transit_time2; double transit = jd_UT + transit_time; // Make calculations for rise and set double rise = 1, set = 1; if (Math.abs(tmp) <= 1.0) { double ang_hor = Math.abs(Math.acos(tmp)); double rise_time1 = celestialHoursToEarthTime * normalizeRadians(ra  ang_hor  lst); double set_time1 = celestialHoursToEarthTime * normalizeRadians(ra + ang_hor  lst); double rise_time2 = celestialHoursToEarthTime * (normalizeRadians(ra  ang_hor  lst)  TWO_PI); double set_time2 = celestialHoursToEarthTime * (normalizeRadians(ra + ang_hor  lst)  TWO_PI); // Obtain the current events in time. Preference should be given to the closest event // in time to the current calculation time (so that iteration in other method will converge) double rise_time = rise_time1; double riseToday2 = Math.floor(jd_UT + rise_time2  0.5) + 0.5; if (jdToday == riseToday2 && Math.abs(rise_time2) < Math.abs(rise_time1)) rise_time = rise_time2; double set_time = set_time1; double setToday2 = Math.floor(jd_UT + set_time2  0.5) + 0.5; if (jdToday == setToday2 && Math.abs(set_time2) < Math.abs(set_time1)) set_time = set_time2; rise = jd_UT + rise_time; set = jd_UT + set_time; } Ephemeris out = new Ephemeris(azi, alt, rise, set, transit, transit_alt, normalizeRadians(ra), dec, dist, pos[0], pos[1], pos[3]); return out; } /** * Corrects geometric elevation for refraction if it is greater than 3 degrees. * @param alt Geometric elevation in radians. * @return Apparent elevation. */ private double refraction(double alt) { if (alt <= 3 * DEG_TO_RAD) return alt; double altIn = alt, prevAlt = alt; int niter = 0; do { double altOut = computeGeometricElevation(alt); alt = altIn  (altOutalt); niter ++; if (Math.abs(prevAltalt) < 0.001 * DEG_TO_RAD) break; prevAlt = alt; } while (niter < 8); return alt; } /** * Compute geometric elevation from apparent elevation. Note ephemerides * calculates geometric elevation, so an inversion is required, something * achieved in method {@linkplain #refraction(double)} by iteration. * @param alt Apparent elevation in radians. * @return Geometric elevation in radians. */ private double computeGeometricElevation(double alt) { double Ps = 1010; // Pressure in mb double Ts = 10 + 273.15; // Temperature in K double altDeg = alt * RAD_TO_DEG; // Bennet 1982 formulae for optical wavelengths, do the job but not accurate close to horizon // Yan 1996 formulae would be better but with much more lines of code double r = DEG_TO_RAD * Math.abs(Math.tan(PI_OVER_TWO  (altDeg + 7.31 / (altDeg + 4.4)) * DEG_TO_RAD)) / 60.0; double refr = r * (0.28 * Ps / Ts); return Math.min(alt  refr, PI_OVER_TWO); /* // Bennet formulae adapted to radio wavelenths. Use this for position in radio wavelengths // Reference for some values: http://ictsyebes.oan.es/reports/doc/ITOAN20032.pdf (Yebes 40m radiotelescope) double Hs = 20; // Humidity % // Water vapor saturation pressure following Crane (1976), as in the ALMA memorandum double esat = 6.105 * Math.exp(25.22 * (Ts  273.15) / Ts) + Math.pow(Ts / 273.15, 5.31); double Pw = Hs * esat / 100.0; double R0 = (16.01 / Ts) * (Ps  0.072 * Pw + 4831 * Pw / Ts) * ARCSEC_TO_RAD; double refr = R0 * Math.abs(Math.tan(PI_OVER_TWO  (altDeg + 5.9 / (altDeg + 2.5)) * DEG_TO_RAD)); return Math.min(alt  refr, PI_OVER_TWO); */ } /** * Sets the illumination phase field for the provided body. * Sun position must be computed before calling this method. * @param body The ephemeris object for this body. */ protected void getIlluminationPhase(Ephemeris body) { double dlon = body.rightAscension  sun.rightAscension; double cosElong = (Math.sin(sun.declination) * Math.sin(body.declination) + Math.cos(sun.declination) * Math.cos(body.declination) * Math.cos(dlon)); double RE = sun.distance; double RO = body.distance; // Use elongation cosine as trick to solve the rectangle and get RP (distance body  sun) double RP = Math.sqrt((cosElong * 2.0 * RE * RO  RE * RE  RO * RO)); double DPH = ((RP * RP + RO * RO  RE * RE) / (2.0 * RP * RO)); body.illuminationPhase = 100 * (1.0 + DPH) * 0.5; } /** * Transforms a common date into a Julian day number (counting days from Jan 1, 4713 B.C. at noon). * Dates before October, 15, 1582 are assumed to be in the Julian calendar, after that the Gregorian one is used. * @param year Year. * @param month Month. * @param day Day. * @param h Hour. * @param m Minute. * @param s Second. * @return Julian day number. * @throws Exception For an inexistent date. */ protected double toJulianDay(int year, int month, int day, int h, int m, int s) throws Exception { // The conversion formulas are from Meeus, chapter 7. boolean julian = false; // Use Gregorian calendar if (year < 1582  (year == 1582 && month < 10)  (year == 1582 && month == 10 && day < 15)) julian = true; int D = day; int M = month; int Y = year; if (M < 3) { Y; M += 12; } int A = Y / 100; int B = julian ? 0 : 2  A + A / 4; double dayFraction = (h + (m + (s / 60.0)) / 60.0) / 24.0; double jd = dayFraction + (int) (365.25D * (Y + 4716)) + (int) (30.6001 * (M + 1)) + D + B  1524.5; return jd; } /** * Transforms a Julian day (rise/set/transit fields) to a common date. * @param jd The Julian day. * @return A set of integers: year, month, day, hour, minute, second. * @throws Exception If the input date does not exists. */ protected static int[] getDate(double jd) throws Exception { // The conversion formulas are from Meeus, // Chapter 7. double Z = Math.floor(jd + 0.5); double F = jd + 0.5  Z; double A = Z; if (Z >= 2299161D) { int a = (int) ((Z  1867216.25) / 36524.25); A += 1 + a  a / 4; } double B = A + 1524; int C = (int) ((B  122.1) / 365.25); int D = (int) (C * 365.25); int E = (int) ((B  D) / 30.6001); double exactDay = F + B  D  (int) (30.6001 * E); int day = (int) exactDay; int month = (E < 14) ? E  1 : E  13; int year = C  4715; if (month > 2) year; double h = ((exactDay  day) * SECONDS_PER_DAY) / 3600.0; int hour = (int) h; double m = (h  hour) * 60.0; int minute = (int) m; int second = (int) ((m  minute) * 60.0); return new int[] {year, month, day, hour, minute, second}; } /** * Returns a date as a string. * @param jd The Julian day. * @return The String. * @throws Exception If the date does not exists. */ public static String getDateAsString(double jd) throws Exception { if (jd == 1) return "NO RISE/SET/TRANSIT FOR THIS OBSERVER/DATE"; int date[] = getDate(jd); String zyr = "", zmo = "", zh = "", zm = "", zs = ""; if (date[1] < 10) zyr = "0"; if (date[2] < 10) zmo = "0"; if (date[3] < 10) zh = "0"; if (date[4] < 10) zm = "0"; if (date[5] < 10) zs = "0"; return date[0]+"/"+zyr+date[1]+"/"+zmo+date[2]+" "+zh+date[3]+":"+zm+date[4]+":"+zs+date[5]+" UT"; } /** * Reduce an angle in radians to the range (0  2 Pi). * @param r Value in radians. * @return The reduced radians value. */ protected static double normalizeRadians(double r) { if (r < 0 && r >= TWO_PI) return r + TWO_PI; if (r >= TWO_PI && r < 2*TWO_PI) return r  TWO_PI; if (r >= 0 && r < TWO_PI) return r; r = TWO_PI * Math.floor(r / TWO_PI); if (r < 0.) r += TWO_PI; return r; } /** * Computes an accurate rise/set/transit time for a moving object. * @param riseSetJD Start date for the event. * @param index Event identifier. * @param niter Maximum number of iterations. * @param sun True for the Sun. * @return The Julian day in UT for the event, 1s accuracy. */ private double obtainAccurateRiseSetTransit(double riseSetJD, EVENT index, int niter, boolean sun) { double step = 1; for (int i = 0; i< niter; i++) { if (riseSetJD == 1) return riseSetJD; // 1 means no rise/set from that location setUTDate(riseSetJD); Ephemeris out = null; if (sun) { out = doCalc(getSun()); } else { out = doCalc(getMoon()); } double val = out.rise; if (index == EVENT.SET) val = out.set; if (index == EVENT.TRANSIT) val = out.transit; step = Math.abs(riseSetJD  val); riseSetJD = val; if (step <= 1.0 / SECONDS_PER_DAY) break; // convergency reached } if (step > 1.0 / SECONDS_PER_DAY) return 1; // did not converge => without rise/set/transit in this date return riseSetJD; } /** * Main test program. * @param args Not used. */ public static void main(String[] args) { System.out.println("SunMoonCalculator test run"); try { int year = 2019, month = 4, day = 1, h = 3, m = 28, s = 52; // in UT !!! double obsLon = 40.64941 * DEG_TO_RAD, obsLat = 54.68653 * DEG_TO_RAD; // lon is negative to the west SunMoonCalculator smc = new SunMoonCalculator(year, month, day, h, m, s, obsLon, obsLat); smc.calcSunAndMoon(); String degSymbol = "\u00b0"; System.out.println("Sun"); System.out.println(" Az: "+(float) (smc.sun.azimuth * RAD_TO_DEG)+degSymbol); System.out.println(" El: "+(float) (smc.sun.elevation * RAD_TO_DEG)+degSymbol); System.out.println(" Dist: "+(float) (smc.sun.distance)+" AU"); System.out.println(" RA: "+(float) (smc.sun.rightAscension * RAD_TO_DEG)+degSymbol); System.out.println(" DEC: "+(float) (smc.sun.declination * RAD_TO_DEG)+degSymbol); System.out.println(" Ill: "+(float) (smc.sun.illuminationPhase)+"%"); System.out.println(" ang.R: "+(float) (smc.sun.angularRadius * RAD_TO_DEG)+degSymbol); System.out.println(" Rise: "+SunMoonCalculator.getDateAsString(smc.sun.rise)); System.out.println(" Set: "+SunMoonCalculator.getDateAsString(smc.sun.set)); System.out.println(" Transit: "+SunMoonCalculator.getDateAsString(smc.sun.transit)+" (elev. "+(float) (smc.sun.transitElevation * RAD_TO_DEG)+degSymbol+")"); /* System.out.println(" Az=+angR: "+SunMoonCalculator.getDateAsString(smc.getAzimuthTime(true, Math.PI+smc.sun.angularRadius))); System.out.println(" Max Elev: "+SunMoonCalculator.getDateAsString(smc.getCulminationTime(true, false))); System.out.println(" Az=0: "+SunMoonCalculator.getDateAsString(smc.getAzimuthTime(true, 0))); System.out.println(" Min Elev: "+SunMoonCalculator.getDateAsString(smc.getCulminationTime(true, true))); */ System.out.println("Moon"); System.out.println(" Az: "+(float) (smc.moon.azimuth * RAD_TO_DEG)+degSymbol); System.out.println(" El: "+(float) (smc.moon.elevation * RAD_TO_DEG)+degSymbol); System.out.println(" Dist: "+(float) (smc.moon.distance * AU)+" km"); System.out.println(" RA: "+(float) (smc.moon.rightAscension * RAD_TO_DEG)+degSymbol); System.out.println(" DEC: "+(float) (smc.moon.declination * RAD_TO_DEG)+degSymbol); System.out.println(" Ill: "+(float) (smc.moon.illuminationPhase)+"%"); System.out.println(" ang.R: "+(float) (smc.moon.angularRadius * RAD_TO_DEG)+degSymbol); System.out.println(" Age: "+(float) (smc.moonAge)+" days"); System.out.println(" Rise: "+SunMoonCalculator.getDateAsString(smc.moon.rise)); System.out.println(" Set: "+SunMoonCalculator.getDateAsString(smc.moon.set)); System.out.println(" Transit: "+SunMoonCalculator.getDateAsString(smc.moon.transit)+" (elev. "+(float) (smc.moon.transitElevation * RAD_TO_DEG)+degSymbol+")"); /* System.out.println(" Az=+angR: "+SunMoonCalculator.getDateAsString(smc.getAzimuthTime(false, Math.PI+smc.moon.angularRadius))); System.out.println(" Max Elev: "+SunMoonCalculator.getDateAsString(smc.getCulminationTime(false, false))); System.out.println(" Az=0: "+SunMoonCalculator.getDateAsString(smc.getAzimuthTime(false, 0))); System.out.println(" Min Elev: "+SunMoonCalculator.getDateAsString(smc.getCulminationTime(false, true))); */ smc.setTwilight(TWILIGHT.ASTRONOMICAL); smc.calcSunAndMoon(); System.out.println(""); System.out.println("Astronomical twilights:"); System.out.println("Sun"); System.out.println(" Rise: "+SunMoonCalculator.getDateAsString(smc.sun.rise)); System.out.println(" Set: "+SunMoonCalculator.getDateAsString(smc.sun.set)); System.out.println("Moon"); System.out.println(" Rise: "+SunMoonCalculator.getDateAsString(smc.moon.rise)); System.out.println(" Set: "+SunMoonCalculator.getDateAsString(smc.moon.set)); /* System.out.println(""); System.out.println("Closest Moon phases:"); for (int i=0; i<MOONPHASE.values().length; i++) { MOONPHASE mp = MOONPHASE.values()[i]; System.out.println(" "+mp.phaseName+" "+SunMoonCalculator.getDateAsString(smc.getMoonPhaseTime(mp))); } double equinox[] = smc.getEquinoxes(); double solstices[] = smc.getSolstices(); System.out.println(""); System.out.println("Equinoxes and solstices:"); System.out.println(" Spring equinox: "+SunMoonCalculator.getDateAsString(equinox[0])); System.out.println(" Autumn equinox: "+SunMoonCalculator.getDateAsString(equinox[1])); System.out.println(" Summer solstice: "+SunMoonCalculator.getDateAsString(solstices[0])); System.out.println(" Winter solstice: "+SunMoonCalculator.getDateAsString(solstices[1])); */ // Expected accuracy over 1800  2200: //  Sun: 0.001 deg in RA/DEC, 0.003 deg or 10 arcsec in Az/El. // <1s in rise/set/transit times. 1 min in Equinoxes/Solstices // Can be used 6 millenia around year 2000 with similar accuracy //  Mon: 0.03 deg or better. // 10s or better in rise/set/transit times. 2 minutes in lunar phases. // In most cases the actual accuracy in the Moon will be better, but it is not guaranteed. } catch (Exception exc) { exc.printStackTrace(); } } }
You may want to draw a figure with the appearance of the Moon disk as seen by an observer on Earth, considering the illumination percentage, the axis position angle, the pole position angle, or even the longitude of the central meridian to show coordinates on the disk. In the following algorithm the first two values calculated are the optical librations lp and bp (lunar coordinates of the center of the disk, or in other words, longitude of central meridian and position angle of pole), then the position angle of axis p, the bright limb angle bl (angle respect north direction from where the sunlight is coming), and the paralactic angle par. In case your chart will be oriented respect celestial north, par can be neglected, but in case you prefer to show the lunar disk as seen by the observer, the apparent inclination of the figure will be (ppar) instead of p (and, of course, par should be added also to the apparent bright limb angle). The illumination fraction is computed in the main program, just take into account the region of the disk illuminated changes after the full Moon phase.
In a first approximation for a drawing you can consider only the illumination fraction, position angle of axis, and paralactic angle. The other angles produce a noticeable but little wobble effect. As an example of that effect, rendered with my JPARSEC library as the rest, you can watch the following 2015 video from minute 01:00: http://www.oan.es/servidorEfem/video/2015.mp4.
You may need to see some charts in Internet to fully understand the angles and to check you apply the rotations correctly, since charting requires some care. Computations are based on the method by Eckhardt, simplified to neglect the physical librations of the Moon. Accuracy is better than 0.5 degrees.
/** * Returns the orientation angles of the lunar disk figure. Illumination fraction * is returned in the main program. Simplification of the method presented by * Eckhardt, D. H., "Theory of the Libration of the Moon", Moon and planets 25, 3 * (1981), without the physical librations of the Moon. Accuracy around 0.5 deg * for each value. * Moon and Sun positions must be computed before calling this method. * @return Optical libration in longitude, latitude, position angle of * axis, bright limb angle, and paralactic angle. */ public double[] getMoonDiskOrientationAngles() { double moonLon = moon.eclipticLongitude, moonLat = moon.eclipticLatitude, moonRA = moon.rightAscension, moonDEC = moon.declination; double sunRA = sun.rightAscension, sunDEC = sun.declination; // Obliquity of ecliptic double eps = meanObliquity + nutObl; // Moon's argument of latitude double F = (93.2720993 + 483202.0175273 * t  0.0034029 * t * t  t * t * t / 3526000.0 + t * t * t * t / 863310000.0) * DEG_TO_RAD; // Moon's inclination double I = 1.54242 * DEG_TO_RAD; // Moon's mean ascending node longitude double omega = (125.0445550  1934.1361849 * t + 0.0020762 * t * t + t * t * t / 467410.0  t * t * t * t / 18999000.0) * DEG_TO_RAD; double cosI = Math.cos(I), sinI = Math.sin(I); double cosMoonLat = Math.cos(moonLat), sinMoonLat = Math.sin(moonLat); double cosMoonDec = Math.cos(moonDEC), sinMoonDec = Math.sin(moonDEC); // Obtain optical librations lp and bp double W = moonLon  omega; double sinA = Math.sin(W) * cosMoonLat * cosI  sinMoonLat * sinI; double cosA = Math.cos(W) * cosMoonLat; double A = Math.atan2(sinA, cosA); double lp = normalizeRadians(A  F); double sinbp = Math.sin(W) * cosMoonLat * sinI  sinMoonLat * cosI; double bp = Math.asin(sinbp); // Obtain position angle of axis p double x = sinI * Math.sin(omega); double y = sinI * Math.cos(omega) * Math.cos(eps)  cosI * Math.sin(eps); double w = Math.atan2(x, y); double sinp = Math.hypot(x, y) * Math.cos(moonRA  w) / Math.cos(bp); double p = Math.asin(sinp); // Compute bright limb angle bl double bl = (Math.PI + Math.atan2(Math.cos(sunDEC) * Math.sin(moonRA  sunRA), Math.cos(sunDEC) * sinMoonDec * Math.cos(moonRA  sunRA)  Math.sin(sunDEC) * cosMoonDec)); // Paralactic angle par y = Math.sin(lst  moonRA); x = Math.tan(obsLat) * cosMoonDec  sinMoonDec * Math.cos(lst  moonRA); double par = 0.0; if (x != 0.0) { par = Math.atan2(y, x); } else { par = (y / Math.abs(y)) * PI_OVER_TWO; } return new double[] {lp, bp, p, bl, par}; }
An old version of the previous code to compute the position of the Moon and the Sun, and the code to compute the orientation angles of the Moon was used by Deep Pradhan to create a very nice demostration program for Java desktop. The project including the files for Eclipse IDE can be downloaded as a zip file. The program shows basic information about the Sun and the Moon and a visual representation of the Moon, allowing to change the position of the observer and the date and time. The date can be changed in steps of 1 day so that the changes in the lunar illumination and the disk orientation during the lunar phase appears very clearly. The program uses precomputed images to show the lunar disk, so the librations angles are not considered.
Deep Pradhan created this demostration program as a first step towards integrating the code presented here in his free Android app Deesha. There are a few changes in this project respect the code presented above, mainly to better integrate the computation of the orientation angles and the graphical representation with the ephemerides of the Sun and the Moon. The code is very clean so you should have no problems testing it. The program is released as free software, so you can use it in your own desktop/Android project.
Many thanks to Deep Pradhan for his great contribution!
NOTE: this program uses the old version of the Sun/Moon calculator.
Deep Pradhan is also working in a iOS version of his Android app. The most recent and updated code is available from this GitHub repository, where the main code follows the Swift code conventions and a little test program is also provided. This is the recommended source, although here you can still download his first Swift v4 code of the Sun/Moon calculator, in which he preserved the code to be as closely as possible to the Java version.
I'm very grateful to him for providing a code that will be very useful also for the community of iOS developers, and for contributing so much to make this post really successful and interesting for all developers. I would love to port my own Android planetarium ClearSky to iOS, but I think is too much hard work.
I would like to emphasize the ports to other languages created by the community. Some of them may come from old versions of the calculator, while others could be from recent versions of the Java code. Please email me any other implementation you may create or find, especially if belongs to the latest code presented above, so I can list it here. If you are a developer of any of the following ports you can check my Bitbucket repository to update your code to the latest version. Diff files between versions can be shown from the source. For instance, the diff file between the July 3, 2020 and the August 7, 2020 versions would be this link.
Swift (iOS):
https://github.com/kanchudeep/SunMoonCalculator
C++ (for Arduino):
https://github.com/ThingPulse/esp8266weatherstation/blob/master/src/SunMoonCalc.cpp
Kotlin (includes planetary positions using a different code from the one presented below):
Python:
There are many different kind of computations that can be solved using iterative and searchlike algorithms. For Androidlike devices it is important to implement them wisely to perform the search as fast as possible. Here I present some example methods that can be useful to obtain the instant of moon phases, seasons, the culmination time, and the time the center of the Sun or the Moon reaches a given azimuth. The culmination time is the instant of the maximum elevation of the Sun/Moon for a given day, which is slightly different from the transit time (instant when the Sun/Mon crosses the meridian exactly towards South, with an azimuth of 180 degrees). The last mentioned method will return the transit time if the azimuth selected is 180 degrees, but can also obtain the inferior meridian transit if you set azimuth to 0 (North), which will be close to the minimum elevation. Please note these methods are not threadsafe, any computation of the Sun/Moon position during the iterations will return wrong results since these iterations change the time. Previous time is recovered at the end before leaving the methods.
With the same logic behind but other considerations it is possible to predict eclipses, for instance. The Saros class in JPARSEC shows an implementation of this idea.
/** * Returns the instant of a given moon phase. * @param phase The phase. * @return The instant of that phase, accuracy around 1 minute or better. */ public double getMoonPhaseTime(MOONPHASE phase) { double accuracy = 10 / (30 * SECONDS_PER_DAY); // 10s / lunar cycle length in s => 10s accuracy double refPhase = phase.phase; double oldJD = jd_UT, oldMoonAge = moonAge; while (true) { double age = normalizeRadians((getMoon()[0]  getSun()[0])) / TWO_PI  refPhase; if (age < 0) age += 1; if (age < accuracy  age > 1  accuracy) break; if (age < 0.5) { jd_UT = age; } else { jd_UT += 1age; } setUTDate(jd_UT); } double out = jd_UT; setUTDate(oldJD); moonAge = oldMoonAge; return out; } /** * Returns the dates of the official (geocentric) Spring and Autumn equinoxes. * @return Dates of equinoxes, accuracy around 1 minute. * @throws Exception If an error occurs. */ public double[] getEquinoxes() throws Exception { double jdOld = jd_UT; double out[] = new double[2]; double prec = 1.0 / 86400.0; // Output precision 1s, accuracy around 1 minute int year = getDate(jd_UT)[0]; for (int i=0; i<2; i++) { int month = 3, day = 18; if (i == 1) month = 9; double jd = toJulianDay(year, month, day, 0, 0, 0); setUTDate(jd); double min = 1, minT = 1; double stepDays = 0.25, lastDec = 1; while (true) { double decAbs = Math.abs(doCalc(getSun(), true).declination); if (decAbs < min  min == 1) { min = decAbs; minT = jd; } if (decAbs > lastDec && lastDec >= 0) { if (Math.abs(stepDays) < prec) { out[i] = minT; break; } stepDays = stepDays / 2; } lastDec = decAbs; jd += stepDays; setUTDate(jd); } } setUTDate(jdOld); return out; } /** * Returns the dates of the official (geocentric) Summer and Winter solstices. * @return Dates of solstices, accuracy around 1 minute. * @throws Exception If an error occurs. */ public double[] getSolstices() throws Exception { double jdOld = jd_UT; double out[] = new double[2]; double prec = 1.0 / 86400.0; // Output precision 1s, accuracy around 1 minute int year = getDate(jd_UT)[0]; for (int i=0; i<2; i++) { int month = 6, day = 18; if (i == 1) month = 12; double jd = toJulianDay(year, month, day, 0, 0, 0); setUTDate(jd); double max = 1, maxT = 1; double stepDays = 0.25, lastDec = 1; while (true) { double decAbs = Math.abs(doCalc(getSun(), true).declination); if (decAbs > max  max == 1) { max = decAbs; maxT = jd; } if (decAbs < lastDec && lastDec >= 0) { if (Math.abs(stepDays) < prec) { out[i] = maxT; break; } stepDays = stepDays / 2; } lastDec = decAbs; jd += stepDays; setUTDate(jd); } } setUTDate(jdOld); return out; } /** * Returns the maximum/minimum elevation time for the Sun or the Moon. * @param forSun True for the Sun, false for the Moon. * @param inferior True to get the minimum elevation time. * @return The Julian day of the culmination instant, which is * only slightly different than the transit. */ public double getCulminationTime(boolean forSun, boolean inferior) { double jdOld = jd_UT; double jd = forSun ? sun.transit : moon.transit; if (inferior) jd += 0.5 * (jd > jdOld ? 1 : 1); double startPrecSec = 20.0 / SECONDS_PER_DAY, endPrecSec = 0.25 / SECONDS_PER_DAY; setUTDate(jd); Ephemeris ephem = forSun ? doCalc(getSun(), false) : doCalc(getMoon(), false); double refAlt = ephem.elevation; while (Math.abs(startPrecSec) > endPrecSec) { jd += startPrecSec; setUTDate(jd); ephem = forSun ? doCalc(getSun(), false) : doCalc(getMoon(), false); if (ephem.elevation < refAlt && !inferior) startPrecSec *= 0.25; if (ephem.elevation > refAlt && inferior) startPrecSec *= 0.25; refAlt = ephem.elevation; } setUTDate(jdOld); return jd; } /** * Returns the instant when the Sun or the Moon reaches a given azimuth. * @param forSun True for the Sun, false for the Moon. * @param azimuth The azimuth value to search for. * @return The Julian day of the azimuth instant. */ public double getAzimuthTime(boolean forSun, double azimuth) { double jdOld = jd_UT; double jd = forSun ? sun.transit : moon.transit; double startPrecSec = 500.0 / SECONDS_PER_DAY, endPrecSec = 0.25 / SECONDS_PER_DAY; setUTDate(jd); Ephemeris ephem = forSun ? doCalc(getSun(), false) : doCalc(getMoon(), false); double refDif = Math.abs(ephem.azimuth  azimuth); while (Math.abs(startPrecSec) > endPrecSec) { jd += startPrecSec; setUTDate(jd); ephem = forSun ? doCalc(getSun(), false) : doCalc(getMoon(), false); double dif = Math.abs(ephem.azimuth  azimuth); if (dif == 0) break; if (dif > refDif) startPrecSec *= 0.25; refDif = dif; } setUTDate(jdOld); return jd; }
The implementation presented above for the position of the Moon can be improved by adding all the terms from the MOON program by Peter DuffetSmith. Instead of using his implementation for the mean elements (secular terms), I use the expansions by Moshier because they are the best ones available, providing excellent results for the ancient past. Here we have the problem of accuracy versus speed, if you prefer accuracy you may want to replace the getMoon method with the code below, but if you are going to use the code in a low powerful device like an Arduino you will prefer the previous code. Since trigonometric functions are slow and down here there are a lot of them you may need some speed tests to decide, especially in case you will use some iterative algorithms. I still recommend using the previous code since the code below provides only a marginal improvement at the cost of many more little periodic terms. However, it is more interesting to replace the code to compute better topocentric positions.
protected double[] getMoon() { // These expansions up to t^7 for the mean elements are taken from S. L. Moshier /* Mean elongation of moon = D */ double x = (1.6029616009939659e+09 * t + 1.0722612202445078e+06); x += (((((3.207663637426e013 * t + 2.555243317839e011) * t + 2.560078201452e009) * t  3.702060118571e005) * t + 6.9492746836058421e03) * t /* D, t^3 */  6.7352202374457519e+00) * t * t; /* D, t^2 */ double phase = ARCSEC_TO_RAD * x; /* Mean distance of moon from its ascending node = F */ x = (1.7395272628437717e+09 * t + 3.3577951412884740e+05); x += (((((4.474984866301e013 * t + 4.189032191814e011) * t  2.790392351314e009) * t  2.165750777942e006) * t  7.5311878482337989e04) * t /* F, t^3 */  1.3117809789650071e+01) * t * t; /* F, t^2 */ double node = ARCSEC_TO_RAD * x; /* Mean anomaly of sun = l' (J. Laskar) */ x = (1.2959658102304320e+08 * t + 1.2871027407441526e+06); x += ((((((((1.62e20 * t  1.0390e17) * t  3.83508e15) * t + 4.237343e13) * t + 8.8555011e11) * t  4.77258489e8) * t  1.1297037031e5) * t + 8.7473717367324703e05) * t  5.5281306421783094e01) * t * t; double sanomaly = ARCSEC_TO_RAD * x; /* Mean anomaly of moon = l */ x = (1.7179159228846793e+09 * t + 4.8586817465825332e+05); x += (((((1.755312760154e012 * t + 3.452144225877e011) * t  2.506365935364e008) * t  2.536291235258e004) * t + 5.2099641302735818e02) * t /* l, t^3 */ + 3.1501359071894147e+01) * t * t; /* l, t^2 */ double anomaly = ARCSEC_TO_RAD * x; /* Mean longitude of moon, re mean ecliptic and equinox of date = L */ x = (1.7325643720442266e+09 * t + 7.8593980921052420e+05); x += (((((7.200592540556e014 * t + 2.235210987108e010) * t  1.024222633731e008) * t  6.073960534117e005) * t + 6.9017248528380490e03) * t /* L, t^3 */  5.6550460027471399e+00) * t * t; /* L, t^2 */ double l = ARCSEC_TO_RAD * x * RAD_TO_DEG; // Now longitude, with the three main correcting terms of evection, // variation, and equation of year, plus other terms (error<0.01 deg) // P. Duffet's MOON program taken as reference for the periodic terms double E = 1.0  (.002495 + 7.52E06 * (t + 1.0)) * (t + 1.0), E2 = E * E; double td = t + 1, td2 = t * t; double M6 = td * JULIAN_DAYS_PER_CENTURY * 360.0 / 6.798363307E3; double NA = normalizeRadians((2.59183275E2  M6 + (2.078E3 + 2.2E6 * td) * td2) * DEG_TO_RAD); double C = NA + DEG_TO_RAD * (275.05  2.3 * td); l += 6.28875 * Math.sin(anomaly) + 1.274018 * Math.sin(2 * phase  anomaly) + .658309 * Math.sin(2 * phase); l += 0.213616 * Math.sin(2 * anomaly)  E * .185596 * Math.sin(sanomaly)  0.114336 * Math.sin(2 * node); l += .058793 * Math.sin(2 * phase  2 * anomaly) + .057212 * E * Math.sin(2 * phase  anomaly  sanomaly) + .05332 * Math.sin(2 * phase + anomaly); l += .045874 * E * Math.sin(2 * phase  sanomaly) + .041024 * E * Math.sin(anomaly  sanomaly)  .034718 * Math.sin(phase)  E * .030465 * Math.sin(sanomaly + anomaly); l += .015326 * Math.sin(2 * (phase  node))  .012528 * Math.sin(2 * node + anomaly)  .01098 * Math.sin(2 * node  anomaly) + .010674 * Math.sin(4 * phase  anomaly); l += .010034 * Math.sin(3 * anomaly) + .008548 * Math.sin(4 * phase  2 * anomaly); l += E * .00791 * Math.sin(sanomaly  anomaly + 2 * phase)  E * .006783 * Math.sin(2 * phase + sanomaly) + .005162 * Math.sin(anomaly  phase) + E * .005 * Math.sin(sanomaly + phase); l += .003862 * Math.sin(4 * phase) + E * .004049 * Math.sin(anomaly  sanomaly + 2 * phase) + .003996 * Math.sin(2 * (anomaly + phase)) + .003665 * Math.sin(2 * phase  3 * anomaly); l += E * 2.695E3 * Math.sin(2 * anomaly  sanomaly) + 2.602E3 * Math.sin(anomaly  2*(node+phase)); l += E * 2.396E3 * Math.sin(2*(phase  anomaly)  sanomaly)  2.349E3 * Math.sin(anomaly+phase); l += E * E * 2.249E3 * Math.sin(2*(phasesanomaly))  E * 2.125E3 * Math.sin(2*anomaly+sanomaly); l += E * E * 2.079E3 * Math.sin(2*sanomaly) + E * E * 2.059E3 * Math.sin(2*(phasesanomaly)anomaly); l += 1.773E3 * Math.sin(anomaly+2*(phasenode))  1.595E3 * Math.sin(2*(node+phase)); l += E * 1.22E3 * Math.sin(4*phasesanomalyanomaly)  1.11E3 * Math.sin(2*(anomaly+node)); l += 8.92E4 * Math.sin(anomaly  3 * phase)  E * 8.11E4 * Math.sin(sanomaly + anomaly + 2 * phase); l += E * 7.61E4 * Math.sin(4 * phase  sanomaly  2 * anomaly); l += E2 * 7.04E4 * Math.sin(anomaly  2 * (sanomaly + phase)); l += E * 6.93E4 * Math.sin(sanomaly  2 * (anomaly  phase)); l += E * 5.98E4 * Math.sin(2 * (phase  node)  sanomaly); l += 5.5E4 * Math.sin(anomaly + 4 * phase) + 5.38E4 * Math.sin(4 * anomaly); l += E * 5.21E4 * Math.sin(4 * phase  sanomaly) + 4.86E4 * Math.sin(2 * anomaly  phase); l += E2 * 7.17E4 * Math.sin(anomaly  2 * sanomaly); double longitude = l * DEG_TO_RAD; double Psin = 29.530588853; if (sun != null) { // Get Moon age, more accurate than 'phase' but we need the Sun position moonAge = normalizeRadians(longitude  sun.eclipticLongitude) * Psin / TWO_PI; } else { // Use the phase variable as estimate, less accurate but this is used only when we don't need an accurate value moonAge = phase * Psin / TWO_PI; } // Now Moon parallax double p = .950724 + .051818 * Math.cos(anomaly) + .009531 * Math.cos(2 * phase  anomaly); p += .007843 * Math.cos(2 * phase) + .002824 * Math.cos(2 * anomaly); p += 0.000857 * Math.cos(2 * phase + anomaly) + E * .000533 * Math.cos(2 * phase  sanomaly); p += E * .000401 * Math.cos(2 * phase  anomaly  sanomaly) + E * .00032 * Math.cos(anomaly  sanomaly)  .000271 * Math.cos(phase); p += E * .000264 * Math.cos(sanomaly + anomaly)  .000198 * Math.cos(2 * node  anomaly); p += 1.73E4 * Math.cos(3 * anomaly) + 1.67E4 * Math.cos(4*phaseanomaly); p += E * 1.11E4 * Math.cos(sanomaly) + 1.03E4 * Math.cos(4 * phase  2 * anomaly); p += 8.4E5 * Math.cos(2 * anomaly  2 * phase)  E * 8.3E5 * Math.cos(2 * phase + sanomaly); p += 7.9E5 * Math.cos(2 * phase + 2 * anomaly) + 7.2E5 * Math.cos(4 * phase); p += E * 6.4E5 * Math.cos(2 * phase  sanomaly + anomaly)  E * 6.3E5 * Math.cos(2 * phase + sanomaly  anomaly); p += E * 4.1E5 * Math.cos(sanomaly + phase) + E * 3.5E5 * Math.cos(2 * anomaly  sanomaly); p += 3.3E5 * Math.cos(3 * anomaly  2 * phase)  3E5 * Math.cos(anomaly + phase); p += 2.9E5 * Math.cos(2 * (node  phase))  E * 2.9E5 * Math.cos(2 * anomaly + sanomaly); p += E2 * 2.6E5 * Math.cos(2 * (phase  sanomaly))  2.3E5 * Math.cos(2 * (node  phase) + anomaly); p += E * 1.9E5 * Math.cos(4 * phase  sanomaly  anomaly); // So Moon distance in Earth radii is, more or less, double distance = 1.0 / Math.sin(p * DEG_TO_RAD); // Ecliptic latitude with nodal phase (error<0.01 deg) l = 5.128189 * Math.sin(node) + 0.280606 * Math.sin(node + anomaly) + 0.277693 * Math.sin(anomaly  node); l += .173238 * Math.sin(2 * phase  node) + .055413 * Math.sin(2 * phase + node  anomaly); l += .046272 * Math.sin(2 * phase  node  anomaly) + .032573 * Math.sin(2 * phase + node); l += .017198 * Math.sin(2 * anomaly + node) + .009267 * Math.sin(2 * phase + anomaly  node); l += .008823 * Math.sin(2 * anomaly  node) + E * .008247 * Math.sin(2 * phase  sanomaly  node) + .004323 * Math.sin(2 * (phase  anomaly)  node); l += .0042 * Math.sin(2 * phase + node + anomaly) + E * .003372 * Math.sin(node  sanomaly  2 * phase); l += E * 2.472E3 * Math.sin(2 * phase + node  sanomaly  anomaly); l += E * 2.222E3 * Math.sin(2 * phase + node  sanomaly); l += E * 2.072E3 * Math.sin(2 * phase  node  sanomaly  anomaly); l += E * 1.877E3 * Math.sin(node  sanomaly + anomaly) + 1.828E3 * Math.sin(4 * phase  node  anomaly); l += E * 1.803E3 * Math.sin(node + sanomaly)  1.75E3 * Math.sin(3 * node); l += E * 1.57E3 * Math.sin(anomaly  sanomaly  node)  1.487E3 * Math.sin(node + phase); l += E * 1.481E3 * Math.sin(node + sanomaly + anomaly) + E * 1.417E3 * Math.sin(node  sanomaly  anomaly); l += E * 1.35E3 * Math.sin(node  sanomaly) + 1.33E3 * Math.sin(node  phase); l += 1.106E3 * Math.sin(node + 3 * anomaly) + 1.02E3 * Math.sin(4 * phase  node); l += 8.33E4 * Math.sin(node + 4 * phase  anomaly) + 7.81E4 * Math.sin(anomaly  3 * node); l += 6.7E4 * Math.sin(node + 4 * phase  2 * anomaly) + 6.06E4 * Math.sin(2 * phase  3 * node); l += 5.97E4 * Math.sin(2 * (phase + anomaly)  node); l += E * 4.92E4 * Math.sin(2 * phase + anomaly  sanomaly  node) + 4.5E4 * Math.sin(2 * (anomaly  phase)  node); l += 4.39E4 * Math.sin(3 * anomaly  node) + 4.23E4 * Math.sin(node + 2 * (phase + anomaly)); l += 4.22E4 * Math.sin(2 * phase  node  3 * anomaly)  E * 3.67E4 * Math.sin(sanomaly + node + 2 * phase  anomaly); l += E * 3.53E4 * Math.sin(sanomaly + node + 2 * phase) + 3.31E4 * Math.sin(node + 4 * phase); l += E * 3.17E4 * Math.sin(2 * phase + node  sanomaly + anomaly); l += E2 * 3.06E4 * Math.sin(2 * (phase  sanomaly)  node)  2.83E4 * Math.sin(anomaly + 3 * node); double W1 = 4.664E4 * Math.cos(NA); double W2 = 7.54E5 * Math.cos(C); double latitude = l * DEG_TO_RAD * (1.0  W1  W2); return new double[] {longitude, latitude, distance * EARTH_RADIUS / AU, Math.atan(BODY.Moon.eqRadius / (distance * EARTH_RADIUS))}; }
Another improvement is to compute an accurate topocentric correction to refer the position of the bodies to the location of the observer instead of the center of the Earth. In doCalc() method you will see the section 'Obtain topocentric rectangular coordinates'. The previous correction did not take into account the Earth ellipsoid neither the altitude of the observer, so you can replace that section with the code below for an accurate topocentric position.
// Obtain topocentric rectangular coordinates double xtopo = x, ytopo = y, ztopo = z; if (!geocentric) { double geocLat = (obsLat  .1925 * Math.sin(2 * obsLat) * DEG_TO_RAD); double sinLat = Math.sin(geocLat); double cosLat = Math.cos(geocLat); double geocR = 1.0  Math.pow(Math.sin(obsLat), 2) / 298.257; double radiusAU = (geocR * EARTH_RADIUS + obsAlt * 0.001) / AU; double correction[] = new double[] { radiusAU * cosLat * Math.cos(lst), radiusAU * cosLat * Math.sin(lst), radiusAU * sinLat}; xtopo = correction[0]; ytopo = correction[1]; ztopo = correction[2]; } double sinLat = Math.sin(obsLat); double cosLat = Math.cos(obsLat);
Where obsAlt is a new variable representing the observer altitude above the sea level in meters. You will have to add it to the program along with obsLon and obsLat, and to the constructor. The sinLat and cosLat variables at the end are used later in doCalc to compute azimuth and elevation, using the geodetic latitude (the input latitude) instead of the geocentric latitude geocLat computed in the topocentric correction. Note the refraction correction applied to the elevation uses as input the ambient pressure in mb and temperature in C. In case you want more control you can add these variables to the class, since they depend on the altitude/weather, and in case the position to compute is for an observer with a very high altitude (for instance when using this code with stratospheric balloons) it could be good to completely disable the refraction correction. That is another problem when you look for a greater accuracy: you have to add more and more terms, variables, or considerations. In this code I did not consider the depression of the horizon, there is a code to compute it inside the discussion section below.
With these changes the accuracy in the position of the Moon will be boosted from 0.015 to 0.005 degrees (just a few arcseconds) in current dates, which means around 12 seconds in the rise, set, and transit times. Moon distances will be accurate to 30 km or better. As you will see below, even around year 2000 B.C. the accuracy will be similar.
I have done some accuracy tests to show the behavior of the code over several millenia. To compute the position of the Sun and the Moon I use the periodic terms included in some old algorithms, and I wanted to know how they perform against the DE404 ephemerides integration by JPL, which is also old but still quite adequate to compute ephemerides for the Sun and the Moon over a long time span. For these tests I used the complete algorithm by Moshier, which is a fit to DE404 with residuals well below the arcsecond.
One of the points was to see if I could identify some additional periodic terms and approximately fit them to improve the ephemerides even more, and it was the case. The image below show the residuals in the ecliptic longitude, latitude, and distance between this program and DE404, computed from year 2000 to year +4000. The interval I used was 100 days or 22000 points.
As you can see, the ecliptic longitude (blue) of the Sun shows a secular deviation from DE404 (or even DE200, since the deviation also appears against the VSOP algorithm). The ecliptic longitude of the Moon shows a periodic deviation that can be at least partially corrected. The discrepancy peaks in the ecliptic longitude is around 36”, but only 15” in the ecliptic latitude (red). Distance errors (black) are around 30 km or better, with peaks of 40 km. This comparison is against the complete MOON program by Peter Duffet described above, for the version with less terms the errors are greater.
The right panels show the results after adding some terms to improve the ephemerides even more. Now the ecliptic longitude of the Sun departs from DE404 less than 1” almost always. Close to year 3500 the comparison cannot longer be done because the algorithm by Moshier is not valid in that epoch. The ecliptic latitude of the Sun is very close to 0 always, so the bad behavior of the red line shown is due to that algorithm. With VSOP the longitude ends with a discrepancy of 3.5”. Note also the distance errors, with typical values of 500 km and peaks around 1000 km. Respect the Moon, the discrepancy peaks respect DE404 can be reduced to 24” in the ecliptic longitude with just one term, or even to 18” with a second one for 700 years around year 2000. With both terms the discrepancy peaks in longitude and latitude are similar, as it should be, and the typical position errors are below 15” respect DE404. Even with a much larger sampling the peak errors are always below 30”. This is well beyond my initial intentions with this simple code.
If you are interested, here are the additional terms I used:
For the Sun, add to the sun_elements array
// Additional terms to 'fit' DE404, to better than 1" almost always new double[] { 4.0 * 48.5, 0.0, 3.3, 8 }, new double[] { 1.0 * 48.5, 0.0, 1.65, 15 }
For the Moon, add to the ecliptic longitude
// Additional terms to add beyond the original program to better 'fit' DE404 l += 14.1983 * Math.cos(2.3232 * t + 0.4964) / 3600.0; if (Math.abs(t) < 7) l += 7.2167 * Math.cos(33.8624 * t  0.5582) / 3600.0;
As reference, the discrepancy respect the most recent integration DE430 is around 150” for year 2000, ten times less than in a previous version of the program when the secular terms used in getMoon method were different and less adequate than those by Moshier for the ancient past or future. This discrepancy may increase with future numerical integrations, but probably not to alarming levels, so the current version of the program can be used several millenia far from year 2000.
I present below a class called PlanetCalculator that you can copy to the same package in which you have the latest version of the SunMoonCalculator class to extend it to support planetary ephemerides. It is not tested in depth but seems to work within the expected accuracy, which is remarkably good for the terrestrial planets and for Neptune. It can be used between years 10003000 A.D., but will give better results between 18002100. The expected accuracy in the rise/set/transit times is of a few seconds in the worst case.
The method is based on the work by J. L. Simon et al. published in A&A 282, 663 (1994). It is based on solving the Kepler equation for elliptic orbits using orbital elements that are modified with time to follow the perturbations between the planets. I took the code from SOFA (in C, easy to adapt to Java) to avoid having to type so many numbers myself. I have used many approximate algorithms over the years, never this one, but I think is a good option since the amount of code lines it requires is limited. The original code rotates the coordinates at the end to equatorial J2000, but since I use ecliptic ones in doCalc() I prefer to skip that and directly account for the precession using a code that works directly in ecliptic coordinates.
To obtain the geocentric vector to the target planet I use the geocentric position of the Sun (with minus sign is the heliocentric position of Earth) and then I add it to the heliocentric position of the planet. Using the EarthMoon barycenter as the position of the Earth (with only the new code) I would get an error of about 5000 km in the output positions, and since I have a nice geocentric getSun() method it is easy to correct it. The planetary aberration or lighttime are also corrected to first order (getPlanet method), something important for the external planets.
This new class also accepts the Sun and Moon as targets, but returns exactly the same output as the Sun/Moon calculator.
/** * A planetary ephemeris calculator to complete the Sun/Moon calculator, also without * using JPARSEC library or any other dependency. * @author T. Alonso Albi  OAN (Spain), email t.alonso@oan.es * @version July 1, 2020 (method calcEarth removed (getSun used), precession in ecliptic system) * @version June 10, 2020 (only basic test done) */ public class PlanetCalculator extends SunMoonCalculator { /** * Main constructor for Sun/Moon calculations. Time should be given in * Universal Time (UT), observer angles in radians. * @param year The year. * @param month The month. * @param day The day. * @param h The hour. * @param m Minute. * @param s Second. * @param obsLon Longitude for the observer. * @param obsLat Latitude for the observer. * @throws Exception If the date does not exists. */ public PlanetCalculator(int year, int month, int day, int h, int m, int s, double obsLon, double obsLat) throws Exception { super(year, month, day, h, m, s, obsLon, obsLat); } /* Planetary inverse masses */ private static double amas[] = { 6023600.0, /* Mercury */ 408523.5, /* Venus */ 328900.5, /* EMB */ 3098710.0, /* Mars */ 1047.355, /* Jupiter */ 3498.5, /* Saturn */ 22869.0, /* Uranus */ 19314.0 /* Neptune */ }; /* * Tables giving the mean Keplerian elements, limited to t^2 terms: * a semimajor axis (au) dlm mean longitude (deg and arcsec) * e eccentricity pi longitude of the perihelion (deg and arcsec) * dinc inclination (deg and arcsec) omega longitude of the ascending node (deg and arcsec) */ private static double a[][] = { new double[] { 0.3870983098, 0.0, 0.0 }, /* Mercury */ new double[] { 0.7233298200, 0.0, 0.0 }, /* Venus */ new double[] { 1.0000010178, 0.0, 0.0 }, /* EMB */ new double[] { 1.5236793419, 3e10, 0.0 }, /* Mars */ new double[] { 5.2026032092, 19132e10, 39e10 }, /* Jupiter */ new double[] { 9.5549091915, 0.0000213896, 444e10 }, /* Saturn */ new double[] { 19.2184460618, 3716e10, 979e10 }, /* Uranus */ new double[] { 30.1103868694, 16635e10, 686e10 } /* Neptune */ }; private static double dlm[][] = { new double[] { 252.25090552, 5381016286.88982, 1.92789 }, new double[] { 181.97980085, 2106641364.33548, 0.59381 }, new double[] { 100.46645683, 1295977422.83429, 2.04411 }, new double[] { 355.43299958, 689050774.93988, 0.94264 }, new double[] { 34.35151874, 109256603.77991, 30.60378 }, new double[] { 50.07744430, 43996098.55732, 75.61614 }, new double[] { 314.05500511, 15424811.93933, 1.75083 }, new double[] { 304.34866548, 7865503.20744, 0.21103 } }; private static double e[][] = { new double[] { 0.2056317526, 0.0002040653, 28349e10 }, new double[] { 0.0067719164, 0.0004776521, 98127e10 }, new double[] { 0.0167086342, 0.0004203654, 0.0000126734 }, new double[] { 0.0934006477, 0.0009048438, 80641e10 }, new double[] { 0.0484979255, 0.0016322542, 0.0000471366 }, new double[] { 0.0555481426, 0.0034664062, 0.0000643639 }, new double[] { 0.0463812221, 0.0002729293, 0.0000078913 }, new double[] { 0.0094557470, 0.0000603263, 0.0 } }; private static double pi[][] = { new double[] { 77.45611904, 5719.11590, 4.83016 }, new double[] { 131.56370300, 175.48640, 498.48184 }, new double[] { 102.93734808, 11612.35290, 53.27577 }, new double[] { 336.06023395, 15980.45908, 62.32800 }, new double[] { 14.33120687, 7758.75163, 259.95938 }, new double[] { 93.05723748, 20395.49439, 190.25952 }, new double[] { 173.00529106, 3215.56238, 34.09288 }, new double[] { 48.12027554, 1050.71912, 27.39717 } }; private static double dinc[][] = { new double[] { 7.00498625, 214.25629, 0.28977 }, new double[] { 3.39466189, 30.84437, 11.67836 }, new double[] { 0.0, 469.97289, 3.35053 }, new double[] { 1.84972648, 293.31722, 8.11830 }, new double[] { 1.30326698, 71.55890, 11.95297 }, new double[] { 2.48887878, 91.85195, 17.66225 }, new double[] { 0.77319689, 60.72723, 1.25759 }, new double[] { 1.76995259, 8.12333, 0.08135 } }; private static double omega[][] = { new double[] { 48.33089304, 4515.21727, 31.79892 }, new double[] { 76.67992019, 10008.48154, 51.32614 }, new double[] { 174.87317577, 8679.27034, 15.34191 }, new double[] { 49.55809321, 10620.90088, 230.57416 }, new double[] { 100.46440702, 6362.03561, 326.52178 }, new double[] { 113.66550252, 9240.19942, 66.23743 }, new double[] { 74.00595701, 2669.15033, 145.93964 }, new double[] { 131.78405702, 221.94322, 0.78728 } }; /* Tables for trigonometric terms to be added to the mean elements of */ /* the semimajor axes */ private static double kp[][] = { new double[] { 69613, 75645, 88306, 59899, 15746, 71087, 142173, 3086, 0 }, new double[] { 21863, 32794, 26934, 10931, 26250, 43725, 53867, 28939, 0 }, new double[] { 16002, 21863, 32004, 10931, 14529, 16368, 15318, 32794, 0 }, new double[] { 6345, 7818, 15636, 7077, 8184, 14163, 1107, 4872, 0 }, new double[] { 1760, 1454, 1167, 880, 287, 2640, 19, 2047, 1454 }, new double[] { 574, 0, 880, 287, 19, 1760, 1167, 306, 574 }, new double[] { 204, 0, 177, 1265, 4, 385, 200, 208, 204 }, new double[] { 0, 102, 106, 4, 98, 1367, 487, 204, 0 } }; private static double ca[][] = { new double[] { 4, 13, 11, 9, 9, 3, 1, 4, 0 }, new double[] { 156, 59, 42, 6, 19, 20, 10, 12, 0 }, new double[] { 64, 152, 62, 8, 32, 41, 19, 11, 0 }, new double[] { 124, 621, 145, 208, 54, 57, 30, 15, 0 }, new double[] { 23437, 2634, 6601, 6259, 1507,1821, 2620, 2115, 1489 }, new double[] { 62911,119919, 79336,17814,24241,12068, 8306, 4893, 8902 }, new double[] { 389061,262125,44088, 8387,22976,2093, 615, 9720, 6633 }, new double[] { 412235,157046,31430,37817, 9740, 13, 7449, 9644, 0 } }; private static double sa[][] = { new double[] { 29, 1, 9, 6, 6, 5, 4, 0, 0 }, new double[] { 48, 125, 26, 37, 18, 13, 20, 2, 0 }, new double[] { 150, 46, 68, 54, 14, 24, 28, 22, 0 }, new double[] { 621, 532, 694, 20, 192, 94, 71, 73, 0 }, new double[] { 14614,19828, 5869, 1881, 4372, 2255, 782, 930, 913 }, new double[] { 139737, 0, 24667, 51123, 5102, 7429, 4095, 1976, 9566 }, new double[] { 138081, 0, 37205,49039,41901,33872,27037,12474, 18797 }, new double[] { 0, 28492,133236, 69654, 52322,49577,26430, 3593, 0 } }; /* Tables giving the trigonometric terms to be added to the mean */ /* elements of the mean longitudes */ private static double kq[][] = { new double[] { 3086,15746,69613,59899,75645,88306, 12661, 2658, 0, 0 }, new double[] { 21863,32794,10931, 73, 4387,26934, 1473, 2157, 0, 0 }, new double[] { 10,16002,21863,10931, 1473,32004, 4387, 73, 0, 0 }, new double[] { 10, 6345, 7818, 1107,15636, 7077, 8184, 532, 10, 0 }, new double[] { 19, 1760, 1454, 287, 1167, 880, 574, 2640, 19, 1454 }, new double[] { 19, 574, 287, 306, 1760, 12, 31, 38, 19, 574 }, new double[] { 4, 204, 177, 8, 31, 200, 1265, 102, 4, 204 }, new double[] { 4, 102, 106, 8, 98, 1367, 487, 204, 4, 102 } }; private static double cl[][] = { new double[] { 21, 95, 157, 41, 5, 42, 23, 30, 0, 0 }, new double[] { 160, 313, 235, 60, 74, 76, 27, 34, 0, 0 }, new double[] { 325, 322, 79, 232, 52, 97, 55, 41, 0, 0 }, new double[] { 2268, 979, 802, 602, 668, 33, 345, 201, 55, 0 }, new double[] { 7610, 4997,7689,5841,2617, 1115,748,607, 6074, 354 }, new double[] { 18549, 30125,20012, 730, 824, 23,1289,352, 14767, 2062 }, new double[] { 135245,14594, 4197,4030,5630,2898,2540,306, 2939, 1986 }, new double[] { 89948, 2103, 8963, 2695, 3682, 1648, 866,154, 1963, 283 } }; private static double sl[][] = { new double[] { 342, 136, 23, 62, 66, 52, 33, 17, 0, 0 }, new double[] { 524, 149, 35, 117, 151, 122, 71, 62, 0, 0 }, new double[] { 105, 137, 258, 35, 116, 88,112, 80, 0, 0 }, new double[] { 854, 205, 936, 240, 140, 341, 97, 232, 536, 0 }, new double[] { 56980, 8016, 1012, 1448,3024,3710, 318, 503, 3767, 577 }, new double[] { 138606,13478,4964, 1441,1319,1482, 427, 1236, 9167, 1918 }, new double[] { 71234,41116, 5334,4935,1848, 66, 434, 1748, 3780, 701 }, new double[] { 47645, 11647, 2166, 3194, 679, 0,244, 419, 2531, 48 } }; /** * Computes the position of a given planet. * @param body The body to compute. * @return Ecliptic longitude, latitude (degrees, of date), distance (AU), and apparent radius (radians), as * given also in other methods for the Sun and the Moon. */ private double[] getPlanet(BODY body) { if (body == BODY.Moon) return super.getMoon(); if (body == BODY.Sun) return super.getSun(); // Not required, but faster // Compute Earth position in the Solar System (= geocentric position of sun) double sunNow[] = getSun(); double eX1 = sunNow[2] * Math.cos(sunNow[0]) * Math.cos(sunNow[1]); double eY1 = sunNow[2] * Math.sin(sunNow[0]) * Math.cos(sunNow[1]); double eZ1 = sunNow[2] * Math.sin(sunNow[1]); double earth[] = new double[] {eX1, eY1, eZ1, 0, 0, 0}; // Using vel = 0 since aberration is already corrected in getSun() // Compute position and velocity of target in the SS double target[] = calcBody(body); // Compute geocentric geometric position double px = target[0]  earth[0]; double py = target[1]  earth[1]; double pz = target[2]  earth[2]; double vx = target[3]  earth[3]; double vy = target[4]  earth[4]; double vz = target[5]  earth[5]; // First order lighttime correction double dist = Math.sqrt(px * px + py * py + pz * pz) * LIGHT_TIME_DAYS_PER_AU; px = px  vx * dist; py = py  vy * dist; pz = pz  vz * dist; // Output geocentric apparent ecliptic position of date double radius = body.eqRadius / AU; double distance = Math.sqrt(px * px + py * py + pz * pz); double longitude = Math.atan2(py, px); double latitude = Math.atan2(pz / Math.hypot(px, py), 1.0); return new double[] {longitude, latitude, distance, Math.atan(radius / distance)}; } /** * Computes the position and velocity of a given body. Based on iauPlan94 program as published * by SOFA (Standards Of Fundamental Astronomy) software collection. Original source: Simon et al. * Fortran code from Simon, J.L., Bretagnon, P., Chapront, J., ChaprontTouze, M., Francou, G., and * Laskar, J., Astronomy and Astrophysics 282, 663 (1994). * * Expected maximum error over the interval 18002100 A.D.: * better than 110" for Jupiter, Saturn, and Uranus. 40" for Mars, <=20" for the others. * Over 10003000 A.D: about 1.5 times the previous errors. Note 150" means around 10s of error in rise/set times, * so the expected error is below 10s for all bodies even in this wide interval. * * @param body The body. * @return Position and velocity, ecliptic of date. */ private double[] calcBody(BODY body) { if (body == BODY.Sun) return new double[] {0, 0, 0, 0, 0, 0}; // For the Sun, located at barycenter int np = body.index; /* Gaussian constant */ double GK = 0.017202098950; /* Maximum number of iterations allowed to solve Kepler's equation */ int KMAX = 10; /* Time: Julian millennia since J2000.0 */ double t = super.t / 10.0; /* Compute the mean elements */ double da = a[np][0] + (a[np][1] + a[np][2] * t) * t; double dl = (3600.0 * dlm[np][0] + (dlm[np][1] + dlm[np][2] * t) * t) * ARCSEC_TO_RAD; double de = e[np][0] + (e[np][1] + e[np][2] * t) * t; double dp = normalizeRadians((3600.0 * pi[np][0] + (pi[np][1] + pi[np][2] * t) * t) * ARCSEC_TO_RAD); double di = (3600.0 * dinc[np][0] + (dinc[np][1] + dinc[np][2] * t) * t) * ARCSEC_TO_RAD; double dom = normalizeRadians((3600.0 * omega[np][0] + (omega[np][1] + omega[np][2] * t) * t) * ARCSEC_TO_RAD); /* Apply the trigonometric terms */ double dmu = 0.35953620 * t, arga = 0, argl = 0; for (int k = 0; k < 8; k++) { arga = kp[np][k] * dmu; argl = kq[np][k] * dmu; da += (ca[np][k] * Math.cos(arga) + sa[np][k] * Math.sin(arga)) * 1e7; dl += (cl[np][k] * Math.cos(argl) + sl[np][k] * Math.sin(argl)) * 1e7; } arga = kp[np][8] * dmu; da += t * (ca[np][8] * Math.cos(arga) + sa[np][8] * Math.sin(arga)) * 1e7; for (int k = 8; k < 10; k++) { argl = kq[np][k] * dmu; dl += t * (cl[np][k] * Math.cos(argl) + sl[np][k] * Math.sin(argl)) * 1e7; } dl = normalizeRadians(dl); /* Iterative soln. of Kepler's equation to get eccentric anomaly */ double am = dl  dp; double ae = am + de * Math.sin(am); int k = 0; double dae = 1.0; while (k < KMAX && Math.abs(dae) > 1e12) { dae = (am  ae + de * Math.sin(ae)) / (1.0  de * Math.cos(ae)); ae += dae; k++; } /* True anomaly */ double ae2 = ae / 2.0; double at = 2.0 * Math.atan2(Math.sqrt((1.0 + de) / (1.0  de)) * Math.sin(ae2), Math.cos(ae2)); /* Distance (au) and speed (radians per day) */ double r = da * (1.0  de * Math.cos(ae)); double v = GK * Math.sqrt((1.0 + 1.0 / amas[np]) / (da * da * da)); double si2 = Math.sin(di / 2.0); double xq = si2 * Math.cos(dom); double xp = si2 * Math.sin(dom); double tl = at + dp; double xsw = Math.sin(tl); double xcw = Math.cos(tl); double xm2 = 2.0 * (xp * xcw  xq * xsw); double xf = da / Math.sqrt(1  de * de); double ci2 = Math.cos(di / 2.0); double xms = (de * Math.sin(dp) + xsw) * xf; double xmc = (de * Math.cos(dp) + xcw) * xf; double xpxq2 = 2 * xp * xq; /* Position (J2000.0 ecliptic x,y,z in au) */ double x = r * (xcw  xm2 * xp); double y = r * (xsw + xm2 * xq); double z = r * (xm2 * ci2); double pos[] = new double[] {x, y, z}; /* Velocity (J2000.0 ecliptic xdot,ydot,zdot in au/d) */ x = v * (( 1.0 + 2.0 * xp * xp) * xms + xpxq2 * xmc); y = v * (( 1.0  2.0 * xq * xq) * xmc  xpxq2 * xms); z = v * (2.0 * ci2 * (xp * xms + xq * xmc)); double vel[] = new double[] {x, y, z}; // Additional computations needed outside the referenced program // Precession directly to ecliptic of date in ecliptic system, without rotation to equatorial pos = precessionFromJ2000(pos); vel = precessionFromJ2000(vel); return new double[] {pos[0], pos[1], pos[2], vel[0], vel[1], vel[2]}; } /** * Precess coordinates from J2000 (t0 = 0) to a given date using the classical method * by Lieske 1979, but in the ecliptic system. * @param p The x, y, z ecliptic positions. * @return Output x, y, z ecliptic positions for the final date. */ private double[] precessionFromJ2000(double p[]) { double t0 = 0, dt = t  t0; double lon = 174.876383888889 + (((3289.4789 + .60622 * t0) * t0) + ((869.8089  .50491 * t0) + .03536 * dt) * dt) / 3600; double ang = ((47.0029  (0.06603  .000598 * t0) * t0) + ((.03302 + .000598 * t0) + .00006 * dt) * dt) * dt / 3600; double prec = ((5029.0966 + (2.22226  .000042 * t0) * t0) + ((1.11113  .000042 * t0)  .000006 * dt) * dt) * dt / 3600; lon *= DEG_TO_RAD; ang *= DEG_TO_RAD; prec *= DEG_TO_RAD; double c1 = Math.cos(lon + prec), c2 = Math.cos(ang), c3 = Math.cos(lon); double s1 = Math.sin(lon + prec), s2 = Math.sin(ang), s3 = Math.sin(lon); double mat1[] = new double[] {c1 * c3 + s1 * c2 * s3, c1 * s3  s1 * c2 * c3, s1 * s2}; double mat2[] = new double[] {s1 * c3  c1 * c2 * s3, s1 * s3 + c1 * c2 * c3, c1 * s2}; double mat3[] = new double[] {s2 * s3, s2 * c3, c2}; double x[] = new double[] { p[0] * mat1[0] + p[1] * mat1[1] + p[2] * mat1[2], p[0] * mat2[0] + p[1] * mat2[1] + p[2] * mat2[2], p[0] * mat3[0] + p[1] * mat3[1] + p[2] * mat3[2] }; return x; } /** * Computes an accurate rise/set/transit time for a moving object. * @param riseSetJD Start date for the event. * @param index Event identifier. * @param niter Maximum number of iterations. * @param body The body. * @return The Julian day in UT for the event, 1s accuracy. */ private double obtainAccurateRiseSetTransit(double riseSetJD, EVENT index, int niter, BODY body) { double step = 1; for (int i = 0; i< niter; i++) { if (riseSetJD == 1) return riseSetJD; // 1 means no rise/set from that location setUTDate(riseSetJD); Ephemeris out = doCalc(getPlanet(body)); double val = out.rise; if (index == EVENT.SET) val = out.set; if (index == EVENT.TRANSIT) val = out.transit; step = Math.abs(riseSetJD  val); riseSetJD = val; if (step <= 1.0 / SECONDS_PER_DAY) break; // convergency reached } if (step > 1.0 / SECONDS_PER_DAY) return 1; // did not converge => without rise/set/transit in this date return riseSetJD; } /** * Calculates ephemeris for a planet. * @param body The target body. * @return The ephemeris for the body. */ public Ephemeris calcPlanet(BODY body) { double jd = jd_UT; // First the Sun (needed for illumination phase) sun = doCalc(getSun()); Ephemeris planet = doCalc(getPlanet(body)); int niter = 15; // Number of iterations to get accurate rise/set/transit times planet.rise = obtainAccurateRiseSetTransit(planet.rise, EVENT.RISE, niter, body); planet.set = obtainAccurateRiseSetTransit(planet.set, EVENT.SET, niter, body); planet.transit = obtainAccurateRiseSetTransit(planet.transit, EVENT.TRANSIT, niter, body); if (planet.transit == 1) { planet.transitElevation = 0; } else { // Update maximum elevation setUTDate(planet.transit); planet.transitElevation = doCalc(getPlanet(body)).transitElevation; } // Compute illumination phase percentage for the planet if (body != BODY.Sun) getIlluminationPhase(planet); setUTDate(jd); return planet; } /** * Main test program. * @param args Not used. */ public static void main(String args[]) { System.out.println("PlanetCalculator test run"); try { int year = 1820, month = 6, day = 9, h = 20, m = 0, s = 0; // UT double obsLon = 4 * DEG_TO_RAD, obsLat = 40 * DEG_TO_RAD; BODY body = BODY.Sun; PlanetCalculator pc = new PlanetCalculator(year, month, day, h, m, s, obsLon, obsLat); Ephemeris planet = pc.calcPlanet(body); String degSymbol = "\u00b0"; System.out.println(body.name()); System.out.println(" Az: "+(float) (planet.azimuth * RAD_TO_DEG)+degSymbol); System.out.println(" El: "+(float) (planet.elevation * RAD_TO_DEG)+degSymbol); System.out.println(" Dist: "+(float) (planet.distance)+" AU"); System.out.println(" RA: "+(float) (planet.rightAscension * RAD_TO_DEG)+degSymbol); System.out.println(" DEC: "+(float) (planet.declination * RAD_TO_DEG)+degSymbol); System.out.println(" Ill: "+(float) (planet.illuminationPhase)+"%"); System.out.println(" ang.R: "+(float) (planet.angularRadius * RAD_TO_DEG * 3600)+"\""); System.out.println(" Rise: "+SunMoonCalculator.getDateAsString(planet.rise)); System.out.println(" Set: "+SunMoonCalculator.getDateAsString(planet.set)); System.out.println(" Transit: "+SunMoonCalculator.getDateAsString(planet.transit)+" (elev. "+(float) (planet.transitElevation * RAD_TO_DEG)+degSymbol+")"); // Can be used in the year interval 10003000, but 18002100 is recommended // Expected accuracy (maximum errors) 18002100: // Jupiter, Saturn, and Uranus: 2 arcminutes, 8s in rise/set/transit times // Mars: 0.7 arcminutes, 3s in rise/set/transit times // Other planets: <0.5 arcminutes, 1s in rise/set/transit times // Sun/Moon: same results } catch (Exception exc) { exc.printStackTrace(); } } }
The history of all code changes to the SunMoonCalculator and PlanetCalculator classes is available at my repository at Bitbucket. The latest code appear at the top in the most recent commit. You can use this repository to check for the latest and most complete version of the source code presented in this page, or to check code changes between your version and the latest one.
Overall the code presented in this post allows to compute the position of the Sun, the Moon, and the planets. The accuracy is around 0.001 degrees for the Sun over several millennia, 0.005 for the Moon also over several millenia (respect the numerical integration DE404), and for the planets 0.05 degrees or better up to 1000 year far from year 2000. The program is quite complete since it provides some physical ephemerides for the bodies, and the rise and set times can be computed with an accuracy of a few seconds in the worst case. I show how to to compute moon phases or the seasons using iterative algorithms. I gave priority to speed over accuracy for the Moon, but including almost all possible effects and corrections with an approximate but still quite accurate approach, and showing how to improve the accuracy even more. The only weak part is the refraction correction, because an accurate approach would be CPU expensive. As a compromise the formulae by Yan (1996) could be an option.
With the hard astronomical part done you can further extend this program to compute other events such as planetary conjunctions or specific configurations from Earth, eclipses, or even occultations of planets by the Moon. Feel free to use this code for any purpose. I only ask you to send me a message so I can track where this code ends, especially if you plan to use this in the space industry.
In case this accuracy level is below your requirements my recent post about numerical integration of planetary orbits could be used to compute ephemerides with a typical accuracy of 0.01” or better respect the planetary ephemerides DE430. A auxiliary file with initial conditions for the integration is needed and sample code is given to create and read one. A 1.5 MB file can provide ephemerides for 3 millennia. The method is less suitable for low powerful devices, but extremely more accurate than this post for the planets. Have a look at the final section of that post. In case of using the methods described in this post to reduce the orbital positions and get the apparent location of the bodies as seen from the Earth the accuracy is limited to 1”, since much more sophisticated reduction algorithms are needed to preserve the nominal accuracy of <0.01”. But that is the purpose of my heavy JPARSEC library.
Discussion
Hi,
I want to calculate time of rise and set at local time. i.e GMT +5:30 how should i do this?
Second i want to obtain moon Nakshatra which is based on moon longitude. How could i obtain that?
Thanks
Hello. Rise/set time in the program above is given in UT (GMT), so you only need to add 5h 30m to the result of the program after you set the correct position for the observer.
Respect moon Nakshatra, I know nothing about Hindu astrology or astrology in general. You need the moon ecliptic longitude, something calculated in the variable longitude in the getMoon method. Seems you need to substract then the ecliptic longitude of the star Spica and do some more things after (according to Wikipedia), but that's something you should look into and do yourself.
Hey Thanks,
How come the Logitude is in thousands e.g. 67444.4444 ? Moon moves from one Nakshatra to other in about 1 Day+ as the total Nakshtras are 27. I just want to map moon's progress across each of them. Just like Moon age calculated in the class.
Longitude is in degrees without 0  360 reduction.
Suppose you have to do the same as in the moon age calculation, something like int nakshtras = (int) (normalizeRadians(longitude * DEG_TO_RAD) * 27 / TWO_PI (+ some constant))
HI The last logic you have shared seems to be working. Will need to do a lot of testing though. If anything needed I will surely write back.
Thank you very much for help.
Am I correct in assuming that for Sun rise/set/transit and Moon rise/set/transit/age calculations the hour, minute, and seconds parameters in “SunMoonCalculator smc = new SunMoonCalculator(year, month, day, h, m, s, obsLon, obsLat);” can be 0 with no adverse affects on the calculations?
Also will entering the obsLon, obsLat parameters with more than 2 decimal places have any effect?
Thank you for this program, it is very useful, and for your help.
Hi Bryan,
In case you mean adverse = wrong results, setting them as 0 is no problem… But the program returns the closest events in time to the input calculation time, so to obtain results referred to a given date, set minute = second = 0, hour = 12  time zone (input date is in UT, you need to set as hour the midday for your location in UT), or maybe get the results for consecutive days and remove those for the dates you don't need. Since the Moon moves considerably within a day, even with the midday trick you could obtain a rise/set time referred to a day before or after 'today', and also, due to that movement, sometimes the Moon has no rise or set in a given date (actually, this happens once per month). Anyway, you can forget all this in case you don't care about getting as rise or set values dates referred to the day before or after the calculation time.
Respect second point, the maximum possible error in the rise/set times with 2 decimal places is around 2.4 seconds, which is within the accuracy of this program. No problem with that.
Regards, Tomas.
So you are saying that for sun rise/set/transit calculations, it is OK to set the hour to (12 + 5) (I am in the USA central time zone which is (UTC  5)) and the min. and sec. to 0. And for the moon rise/set/transit calculations the same is true but I will not get times for events that happen on the previous/next day (I want to only display events for the current day (or a specified day) and if there is no moonrise that day I will display something like “none” or “no moonrise today”, if I understand you correctly, setting the time as described above will work just fine for my purposes.
Regarding the lat/lon decimal places, the reason I was asking is because I will be feeding output from the device's GPS to the equation and wondered if I should trim the decimal places, or just leave it as is (as long as it will fit into a double of course).
It is fine to do so as a start, but there is no guarrantee that the output rise/set will be always for the current date, considering also that you have to correct back from UT to the USA time zone. You will have to repeat calculations for +/ 1 day (maybe better a little less for the Moon, say +/ 20 hrs) and select from the three rise/set times the one corresponding for the specified day, if any.
This is a fairly simple project and I don't need really super accuracy (but I don't want it very far off either, say +/ 5 min. or so (but I would not reject something because it was to accurate, just if it got to complicated)), and I am going to ignore user altitude (although I could get it from the GPS (I am planning to use the location services, I wonder if I can get alt. from the network?)).
Would it be better to give the full current local date and time, or would I still have to repeat the calculations for yesterday, today, and tomorrow (and +/ 20 hrs for the moon)?
would the check for “today” be something like: run the 3 calculations and convert them all to local timezone and check for: [not a negative hour OR hour not greater than 24]?
Just curious, could the calculations be improved even more by doing something like: today, today +/ 10 (or say 14 for the sun) hours, and today +/ 20(or 24 for the sun) hours for a total of 5 calculations, then finding the ones for “today” (especially for the moon calculations)?
I forgot to mention that I do understand about converting back to local time from UTC.
Which brings up a question, would it be better to convert local time to UTC before sending the date/time to the method?
 User altitude is not important, I have used interpolation in an Earth's bump map, but it is complicated and unreliable.
 Use always +/ 20 hrs for Moon and +/ 1 day for Sun. To use the midday trick I think is less important, just test yourself.
 Rise and set times are given as Julian days, in units of days. Add time zone to those variables to get the Julian day in local time, and convert it back to year/month/day (getDate() method).
 No! The program automatically considers the movement of the Sun and the Moon, so it already does that. The real accuracy is a few seconds, in case you change the hour the rise/set times will be exactly the same (an accurate), or different only in case the closest event is for a different day.
 Well, convert local time to UTC before using the method is in fact what you must do to obtain the correct position of the body for that specific time. In case you only want rise/set you should also do that to have more chances of getting the rise/set for the specified date. Since the rise/set are also UTCs you must convert them back to local.
– I am not sure that I have the questions and answers in the right order, so here is what I think you are saying, if you would please arrange the questions with the correct answers I would appreciate it. Sorry I am kind of new at this.
This is a fairly simple project and I don't need really super accuracy (but I don't want it very far off either, say +/ 5 min. or so (but I would not reject something because it was to accurate, just if it got to complicated)), and I am going to ignore user altitude (although I could get it from the GPS (I am planning to use the location services, I wonder if I can get alt. from the network?)).
 User altitude is not important, I have used interpolation in an Earth's bump map, but it is complicated and unreliable.
– I am pretty sure I have this one in the right place, I was just mentioning that I was going to ignore the alt. and I still plan to do so, the rest was just conjecture or thoughts about possible future changes.
Would it be better to give the full current local date and time, or would I still have to repeat the calculations for yesterday, today, and tomorrow (and +/ 20 hrs for the moon)?
 Use always +/ 20 hrs for Moon and +/ 1 day for Sun. To use the midday trick I think is less important, just test yourself.
– I am pretty sure I have this one in the right place to.  I do plan to do this.
would the check for “today” be something like: run the 3 calculations and convert them all to local timezone and check for: [not a negative hour OR hour not greater than 24]?
– I don't see an answer to this one, or I completely don't understand the answer (this is one of my most important questions, how do I determine which of the 3 calculations is the right one to use?).
Just curious, could the calculations be improved even more by doing something like: today, today +/ 10 (or say 14 for the sun) hours, and today +/ 20(or 24 for the sun) hours for a total of 5 calculations, then finding the ones for “today” (especially for the moon calculations)?
– By now I am completely lost.
Then there is the comment I sent later:
I forgot to mention that I do understand about converting back to local time from UTC. Which brings up a question, would it be better to convert local time to UTC before sending the date/time to the method?
– I don't think you said anything about this one, and I don't think it is necessary to do so, it was just for your information.
Just FYI I am allowing the user to use a date picker to choose dates in the future and past.
← Would the check for “today” be something like: run the 3 calculations and convert them all to local timezone and check for: [not a negative hour OR hour not greater than 24]?
→ Rise and set times are given as Julian days, in units of days. Add time zone to those variables to get the Julian day in local time, and convert it back to year/month/day (let getDate() method do it for you). getDate() will never return a negative hour or a value greater than 24, so use getDate(smc.sunRise ÷ timeZoneInHours / 24.0) to get the local date and time of the sun rise instead of UT. Use ÷/ 1 day and check all 3 values of the returned day for the day you want, if any is present. In the main() I simply use getDate() to output the rise time in UT, just do this to output the rise time in local time.
← Just curious, could the calculations be improved even more by doing something like: today, today ÷/ 10 (or say 14 for the sun) hours, and today ÷/ 20(or 24 for the sun) hours for a total of 5 calculations, then finding the ones for “today” (especially for the moon calculations)?
→ The program automatically considers the movement of the Sun and the Moon, so it will returns the same rise/set times when you change the input time, unless the closest event for that time is for a different day. Hence with one 'try' with ÷/ 1 day for the Sun or ÷/ 20 hours for the Moon is enough.
← I forgot to mention that I do understand about converting back to local time from UTC. Which brings up a question, would it be better to convert local time to UTC before sending the date/time to the method?
→ (If by method you mean this program) This program only use UTC, so to convert local time to UTC before using it is in fact what you must do to obtain the correct position of the body for that specific time. In case you only want rise/set you should also do that to have more chances of getting the rise/set for the specified date after converting them back to local time. Although the calculations for ÷/ 1 day is the way to be sure.
–> getDate() will never return a negative hour or a value greater than 24, so use getDate(smc.sunRise ÷ timeZoneInHours / 24.0) to get the local date and time of the sun rise instead of UT.
I realize that getDate won't return a negative number, but when you apply the timezone offset to it you can get a negative number, for instance my timezone offset is 5 (CST 6 + daylight savings time offset 1 = 5) and getDate returns a UTC hour of 1 you get 1 + 5 = 4.
getDate(julian_day in UT) → year, month, day in UT
getDate(julian_day in LT) → year, month, day in LT
No need to do yourself any time zone correction to the output of getDate, do it in the input Julian day.
my Main code is pretty much just ripped from your Main method. Here is what I am doing (I have left out a lot of code that is not related for brevity). It builds and the app runs on an android phone, everything works, but it won't give me sun/moon rise/set times that even come close to the U.S. Naval Observatory Sun and Moon Data for One Day. (for some reason the sun transits are fairly close, I can't remember about the moon, I think it is close also). Just fyi my app uses a datpicker to select new dates after it starts up, and uses location services to change the obsLat and obsLon, these both work as expected.
I appreciate all of your help.
public class MainActivity extends AppCompatActivity implements GoogleApiClient.ConnectionCallbacks, GoogleApiClient.OnConnectionFailedListener {
… the rest of your setup code if I left any thing out
… more of the same as above
… other unrelated code that uses moonAge to change the picture displayed to the correct
some other things I have tried sunset
I said it works on an android phone, although I have run it on a phone once, I am normally running it on an android studio virtual device (Genymotion, because it allows me to change the GPS coordinates) for testing.
I have also used getDateAsString() in the same configurations as above (and probably some more). And I have used my own getTimeAsString() which is just a copy of your getDateAsString() that I modified to only return the hour min. and sec.
Seems you don't catch the thing, I guess some mistakes like
date2 = SunMoonCalculator.getDate(smc.sunRise +/ 24);
instead of
date2 = SunMoonCalculator.getDate(smc.sunRise +/ 1);
But even that code makes no sense without calling the entire instance to do the calculations. Calendar, AM/PM, and DST things should also be done with extreme care. Here is my code, tested against my main library (8s of error in the Moon, 1s in the Sun) and against some web sources.
public static void main(String args[]) {
}
At the end of the forceResultsForToday block I would add this to also solve the no Moon event problem once a month
Ok thanks, I will give it a try and let you know.
That helped a lot, but I am still getting 1+ to 2+ hours off on the sun/moon set/rise/transit times, and I am not sure about the moon age but it seems to be close enough that my program was giving the right moon phase (at least for today).
some possible causes: you wrote, ” Check time zone is in hours \n int timeZone = 5 + 1; (+1 since it is summer)”, although you have the right number here, the way you arrived at it is incorrect, the time zone here is 6 (CDT), because it is summer here (daylight savings time) we add 1, 6 + 1 = 5 (CDT + daylight savings time = 5) this in itself is not a problem, but you are using lat lon from Washington DC ” Check these values are in radians, they are for Washington DC \n double obsLon = 77.03 * DEG_TO_RAD, obsLat = 38.9 * DEG_TO_RAD;” which is EDT or 5 (with out daylight savings time), or currently, with daylight savings time accounted for it is 4. I don't know if this is all of the problem but it must have an effect of about an hour or so. I am doing * DEG_TO_RAD on my lat and lon's.
I test against: http://aa.usno.navy.mil/data/docs/RS_OneDay.php  the U.S. Naval Observatory
Build finally finished, my results compared to http://aa.usno.navy.mil/data/docs/RS_OneDay.php  the U.S. Naval Observatory are:
lat  N41° 4 lon W94° 22' May 04 2016 CST + DST = 5
my data
sun 6:41:59 14:14:27 21:18:27 moon 5:32:32 12:22:33 19:22:36 moon age 25.49805876518279 (waning crescent)
USNO data
sun 6:11 am 1:14 pm 8:18 pm moon 5:00 am 11:22 am 5:53 pm moon age unavailible (waning crescent)
Doing that myself with USNO I find the same:
Sun
Begin civil twilight 05:41
Sunrise 06:11
Sun transit 13:14
Sunset 20:18
End civil twilight 20:48
Moon
Moonrise 05:00
Moon transit 11:22
Moonset 17:53
And the results with the program presented here (unmodified), with the previous main method, only changing location and with TZ+DST = 6+1 = 5 for the new place:
Sun
Rise: 2016/5/4 6:10:56
Transit: 2016/5/4 13:14:11
Set: 2016/5/4 20:18:6
Moon
Rise: 2016/5/4 4:59:50
Transit: 2016/5/4 11:21:45
Set: 2016/5/4 17:52:51
Which are compatible with USNO up to the accuracy stated here: 12s with Sun, around 10s with Moon.
You mention you modified something to simplify the program. Check your program, you must be doing something wrong.
The changes I have made are:
 added the datepicker and use the calendar function rather than date  it displays the date correctly on the screen, so it seems to be working correctly. Just one question, the datepicker and the “static Calendar calendar = Calendar.getInstance();” gives the month in an odd way, it gives January = 0, February = 1 etc. Will you r code handle this or should I correct for it? Like month = month + 1.
 added location services to get the lat/lon (I do apply the DEG_TO_RAD to them)  this is also correctly displayed on the screen, lat/lon are given with quit a few decimal places (it fits into a double).
the above additions only change the values in the variables used by the calculation.
 added code (in my version of your Main method) to use moonAge to change the picture displayed to the correct moon phase picture, and display the appropriate text (full moon, first crescent, etc…)  this also displays correctly on the screen. here is the code:
added a new method to your calculation code to return only the hours, minuets, and seconds, I just copied your getDateAsString method, renamed it, and changed the return statment. here is the code:
I think it is the problem with the incorrect month, but I have not had a chance to test it yet. Thank you.
The sunrise and moonrise are both off by about 1 hour and 29/30 minutes, everything else is good to within a minute or less. for today (May, 6 2016) my sunrise was 5:29:35  USNO is 6:10, and my moonrise is 5:51:33  USNO is 6:21. I am adding 1 to the month before calculating. I will try it the other way and get back to you.
By not adding 1 to the month, before calculating, the sunrise is off by about 14 min. and the sun transit is good, sun set is off by about an hour, moonrise is only off by about 1 min., moon transit is off by about an hour, and moon set is off by about an hour. Looks to me like I should add 1 to the month before calculating.
When I add 1 to the “month” before calculating, the sunrise I am getting seems to be “Begin civil twilight” (within a minute or so) and the moonrise seems to be something similar, like a “moon Begin civil twilight” of sorts.
Everything else is within a minute or so, good enough for me.
I can't figure out what is up with the sun/moon rises though, I believe that setting “smc.setTwilight(SunMoonCalculator.TWILIGHT.TWILIGHT_CIVIL);” and “boolean forceResultsForToday = true;” and adding 1 to the “month” before calculating, are the settings I should be using, is this correct?
By the way the GPS gives lat/lon out to 14 decimal places (with 2 places, or 2 places + sign (+/), before the decimal).
You should use HORIZON_34arcmin as twilight, the default value (just do not set it). CIVIL one is when the Sun is 6 degrees below the horizon (documented at the top). Then check the (moon) sunrise/sunset with the values I found for your location lat  N41° 4 lon W94° 22' May 04 2016 CST + DST = 5.
Seem's everything else is fine!
Except for the moonage that fixed it, if you round the seconds it is to the minute (the ones that are a minute less than the online source have a seconds less than 30, the ones that the minutes match have a seconds greater than 30). Good enough for my purposes.
The moonage is off by 2 days, online sources say it should be 4.122 days, but I get 2.5325…. days, for May 10 2016 17:57 CDT (GMT  5). It does not make much difference if I do: “double moonAge = smc.moonAge;” or “double moonAge = smc.moonAge + timeZone / 24.0;”, I have no idea what is causing this.
I have found that my emulator was set to EDT (GMT  4) so you DO have to do something like: “static int daylightSavings = (calendar.get(Calendar.DST_OFFSET) / (1000 * 60 * 60) % 24); static int timeZone = (calendar.get(Calendar.ZONE_OFFSET) / (1000 * 60 * 60) % 24) + daylightSavings; total time zone offset”. Just FYI I have edited my “getTimeAsString” method to pad the times with zeros, 01:03:05 instead of 1:3:5. here is the code: “public static String getTimeAsString(double jd) throws Exception { if (jd == 1) return “NO RISE/SET/TRANSIT FOR THIS OBSERVER/DATE”; int date[] = SunMoonCalculator.getDate(jd); return String.format(”%02d”, date[3])+”:”+String.format(”%02d”, date[4])+”:”+String.format(”%02d”, date[5]);”
You should test the maths in an external and simple desktop program. You would then see the moon age returned is 4.66 days, not 2.53. Maths are maths, with no bugs you must obtain the same result as in the desktop program, not one minute or 2 days off.
4.66 for May 10 or May 11?
Are you using a specific time of day?
I am getting 4.681706 at around 1:00pm on May 11. Internet sources are giving me 5.20 and 4.895.
If I understand correctly, since the moonage is just a count of days since the last new moon, no timezone correction is necessary, is this correct?.
re: You should test the maths in an external and simple desktop program. You would then see the moon age returned is 4.66 days, not 2.53. Maths are maths, with no bugs you must obtain the same result as in the desktop program, not one minute or 2 days off.
I use your original, unaltered code and compile and run it in regular java, besides using internet sources, to check my output, that is why I know I am off and have a problem with my code.
I understand that the math is the math, and that your program works, and that I should not be off by any amount, but I was just hoping that by seeing how far off I was it would give you a clue as to what I was doing wrong.
right now I seem to be about 5 hours off, I get 5.06222 and the internet source says I should get 5.285.
All I am doing is getting the age, using it to determine the moon phase (without altering it in any way) (BTW there are some errors in that code above, I can repost if you want me to), and printing it to the screen:
double moonAge = smc.moonAge; … String moonAgeString = ”” + (float) (moonAge); … moonAgeText.setText(moonAgeString); this is just the android version of a print() or println() I REALLY appreaciate ALL of your help, thank you.
For May 11 2016 22:25:00
Running your unaltered program in java gives 5.6601515.
Internet source gives 5.329.
Using “double moonAge = smc.moonAge  timeZone / 24.0;” (Instinct tells me this is not the way to do this, it was just an experiment because I seem to be off by about 5 hours)
I get 5.3201737, this would be good enough for me if the internet source is correct.
4.66 was for the date and time you wrote, May 10 2016 17:57 CDT (GMT  5), which is May 10, 22:57 UT in this program. It needs no correction. Time since previous new moon was 4.14 days at that time.
There are different methods to compute moon age (or several ways to define it), giving differences of a few hours. You may be comparing with an internet source in which you provide CDT when it expects GMT, or in which a different computing method is used.
But shouldn't your code in java give something pretty close to the same code in android studio? It does not give the same answer even when I add 5 hours like above (I effectively used: “double moonAge = smc.moonAge + 5 / 24.0;”).
I have to ADD 5 hours to even come close, with out adding 5 hours it would be something like 5.1118404 while your unaltered program in java (my simple desktop program) gives 5.6601515. shouldn't the same method give similar results without having to adjust it?
I don't know what happened, but after I shut everything down last night, and restarted this morning, I restarted android studio and removed the +5 hours adjustment, and rebuilt the app, and now things seem to be working right ????
one internet source (http://cycletourist.com/moon/ (this is the one I think is the most accurate)) gives 5.894
your code on the desktop java gives 6.189334
android on the emulator gives 6.188877
android on the actual phone gives 6.1887045
the other internet source (http://www.maxx.moongiant.com/phase/05/12/2016) gives 6.20
your java code (desktop program) and both of the android codes (same code just running on different devices) are reasonably the same, and we are between the two internet sources, I think it is working now.
UPDATE: refreshed the date on the phone and emulator and now I get 5.6696715, rebuilt the app and I get about the same thing. I don't know what is going on, think I will just give up on it.
UPDATE 2: the only thing I can think of is the first test was at 11:50 am and the last test was at 12:something but that should not make any difference, the desktop program still gives a good number 6.2151933, I still don't know what is going on!!!
Thank you for ALL of your help.
Should I be inputting the time in 24 hour format? It is currently in am/pm format ie 1:00 pm is 1:00 not 13:00.
That was it, java.util.calendar returns 0 for both midnight and noon so
int am_pm = calendar.get(Calendar.AM_PM);
…
if (am_pm == 1) { h += 12; }
fixed it.
Both android apps match the desktop program.
Correction:
int h = calendar.get(Calendar.HOUR); gives the hour in 12 hour format (0  11) so the above correction is needed. int h = calendar.get(Calendar.HOUR_OF_DAY); gives the hour in 24 hour format (0  23) so the above correction is NOT needed.
Hi Tomás, First of all I would like to thank you for Sun Moon Calculator, very useful indeed. I am writing my own Android app, for my own private use and have implemented it quite successfully. I'm learning how to program in Java/Android and using the app writing process as a self learning tool, I'm slowly getting there :) I have also implemented your extra code to draw a representation of the Moon's disk, but I have a question, probably a silly one, if it is I'm sorry and apologise, like I say I'm still learning. Anyway, the question is,,, are the results returned for p and par in radians? As I am trying to accurately display the parallactic angle of the Moon's disk. Thanks again for this code it is fantastic, and also sorry again if the question I asked is considered to be silly,,, please forgive this old man trying to learn to programme and keep his brain active :) Regards, Tim.
Thanks a lot Tim, programming in Java was indeed also very slow for me when I started, but suddenly, after a year, I started to progress a lot faster. With Android is the same, but the problem here are the amount of 'tricks' or workarounds needed in this platform to solved the problems. As with many things, you will master it if you spend enough time.
And yes, all the angles returned in that code are in radians.
Regards, Tomás
Thank you for your response Tomás, Yes I found Java development hard going to start with, but I'm slowly getting the hang of it, mainly by examining other people's code, yours included, and experimentation. In radians, that's great, thanks Regards, Tim.
P.s. I'm also going to experiment with your full Jparsec for Android in the near future, I would like to eventually try and write a full blown planetarium, again just for my own use. Thanks again for your work, it is greatly appreciated. Regards, Tim.
I know you can get age of the moon. Is it possible to calculate the phase and the Illumination (12%, 98% etc.) as well?
Age of the Moon is usually related to difference between the ecliptic longitude of Sun and Moon, and it isn't strickly equivalent to phase and illumination.
Illumination phase can be computed from moon elongation (angular distance Sun  Moon):
Angular distance can be approximated from RA/DEC of both objects (you need to compute positions of both objects):
although you could use azimuth/elevation, its better to use RA/DEC because elevation could be affected by refraction.
Thanks!
Where can I get the values of RA1,RA2,DEC1,DEC2? Thanks!
The labels 1 and 2 are for the Sun and Moon. Just define global double variables sunRA, sunDEC, moonRA, moonDEC (as already done with azimuth and elevation) and, inside the method calcSunAndMoon(), after the call to method doCalc() for the Sun and Moon sections, add in each sections:
as you can see, RA and DEC are the elements 6 and 7 in the output array 'out', which I decided to ignore.
Thanks! did your code. Looks like the result matches internet!
The azimuth and elevation seem to be way off. Although other results are correct. You can use this US military site to compare your Azimuth: http://aa.usno.navy.mil/data/docs/AltAz.php. Regards.
Not possible, please give a full test example for that page, but please check first that date, location coordinates and time zone are the same in the program and that page.
Sorry. I was using local PC time which caused incorrect results. Thanks!
Hi, I just want to calculate next moon phase exact time, like for next event is new moon at lon 79.0368, lat 21.1039, 15 May 5.18 pm and so on. Can you help me with calculations?
A simple loop using the Moon age should be fine. Assuming the instance is called smc, desired accuracy is 10s, and you want to compute the event closest in time, something like:
where phase is the phase to compute (03).
Thanks to Tomás Alonso Albi, I got the idea now.
Hi, I found a bug in moon Azimuth angle.
time 1: 30 May 2018 11:55:01 = Azimuth 260 time 2: 30 May 2018 12:16:01 = Azimuth 258 time 3: 30 May 2018 12:55:01 = Azimuth 140 time 4: 30 May 2018 13:55:01 = Azimuth 98 time 5: 30 May 2018 14:55:01 = Azimuth 97 time 6: 30 May 2018 16:55:01 = Azimuth 101
Could you look into it?
I have checked the program carefully and others have also thought they found some kind of errors. Please check all inputs you use, specially the observer location and units (radians) and the time zone you use, if any. The program uses Universal Time and likely you are comparing results with other sources without using on both the same location for the observer and computation time (see the answer by Bruce Jin right before you wrote).
In case you still have problems, let me know all the input data you use.
I used same Demonstration program in this blog as source for my computation. You can use same time cases as I used in previous comment to reproduce on your side to check. It may reproduce on any locations.
Sorry but cannot find any problem from my side. You should first explain exactly why you think there is a bug in the azimuth. Please take into account azimuth depends on the observer position, what am I supposed to reproduce? I have checked against my own JPARSEC library, copying the code as it is posted, and everything is ok.
I ported this Java class to Cplusplus for use on an Arduino. With just a few tweaks it would also work in a plain C/Cplusplus environment.
We made it available as part of our ThingPulse Weather Station library at https://github.com/ThingPulse/esp8266weatherstation.
Proper credit was given in three places: the README, the release notes, and of course the code
Glad you find the class useful, and thanks for the detailed credits
Hi Tomás, I was wondering if it is possible to use your class to work out when solstices and equinoxes occur. If it is possible would you be kind enough to show an example piece of code please? Regards Tim
I have updated the program, with some reorganization of the code, to support those kind of calculations and to add also the lunar phases.
Did you ever think about hosting this code in a public repository like GitHub? It would be much easier to reason about, provide patches (if necessary) or allow people to contribute new features.
I didn't expect this post to become popular, but I agree it would be the best approach now. I have many other things to do, so with more free time I will consider the idea in case things here continue growing.
Thank you, much appreciated. Regards, Tim
It's great. Thank you very much. It's very useful for creating Android project and save a lot of development time.
Hello. Im try to use moon.eclipticLongitude But get some unnormal value. I wait degree 0360, but get 90714.26460707808 What I do wrong? Thanks a lot
–smc.moon.eclipticLongitude * RAD_TO_DEG : 90714.26460707808 —Transit: 2018/10/21 19:18:12 UT
I didn't apply the reduction 02PI on that value. Just call the included function normalizeRadians on that value before converting radians to degrees.
Hello on the page you published the SunMoonCalculator.java very nice, But I need javascript code for html. I'm very glad you send this code to my mail address. thanks
Sorry this is Java code, you will have to port it to javascript. You can partially save work using http://www.jsweet.org/.
I want to get longitude of sun and moon. How can i get ?
Can't understand what you ask for. Celestial longitude = right ascension, suncentered longitude = ecliptic longitude, and the longitude respect the horizon of the observer = azimuth. The program returns all that in the sun/moon objects, just use it …
Hi!
I have some problem with this particular algorithm above.
A year ago I did made some electronic device for home use (sun, moon, weather and etc). Took the algo, rewrote it in C++ and then in C and it work absolutely fine.
Recently I've decided to make an Android program: nothing special, just simple clock with moonphase, sunrise and sunset, but now in Java. And it didn't work as expected.
For example take my location: lat = 54.68653, lon = 40.64941, timezone = +3, year = 2019, month = 4, hour = 63, min = 28, sec = 52, put it in online java compile (jdoodle.com in my case) and get some extraordinary results for sun: Sun Rise: 2019/04/01 02:50:46 UT Set: 2019/04/01 15:53:11 UT Transit: 2019/04/01 09:21:22 UT (elev. 39.824867?)
But the actual data is  sunrise is 03:13 UT and sunset in 16:14 UT
What may be the case? Thanks!
I suppose you mean the 'extraordinary results' are the correct ones, I can reproduce them both on my PC and in jdoodle. I also suppose the 'actual data' are the wrong results you get on Android.
You are maybe using something from Java (Calendar class or whatever) that could explain the differences, not entering at the end exactly the same parameters. The program has 0 dependencies intentionally so you should find exactly the same results on both platforms. Another idea in case you have some kind of loop or thread in your Android example code (likely since you did a clock) is to recreate the SunMoonCalculator instance all the time (or call setUTDate), since some methods could leave a wrong UTC date inside the instance. In my opinion, something must be wrong when you call the calculator from your Android code, although I admit I never tried this code on Android. There are above some messages from Bryan with a similar problem, at the end it was the Calendar class and he finally found exactly the same results.
Thanks for reply, but the problem is the opposite one.
I meant that 'extraordinary results' are the not the correct ones (which I got from the algo). And the 'actual data' are results which I got from other services (timeanddate.com and time.is).
Yes, algo works every time the same way on any platform (android, jdoodle, etc.), don't argue with that, but the results are not correct.
I don't use any Calendar class or whatever, just hardcode time and location values (lat = 54.68653, lon = 40.64941, timezone = +3, year = 2019, month = 4, hour = 63, min = 28, sec = 52) and got Sun Rise at 2019/04/01 02:50:46 UT. But the value represented by any other servises  sunrise is at 03:13 UT. It's around 25 minutes of difference.
More over, the c++ version is showing no difference from other services (timeanddate.com and time.is).
I'm confused.
So the point is you think the web services and c++ port are ok and the output from the calculator is not (I think is the opposite!). The code in this class is not wrong, I can guarrantee that. Checked with my complete JPARSEC library the difference is 2s at most.
I think the only answer is that the location you use in the web service / c++ port is not the same as the one you use in the java code. Please check carefully the input you provided above, specially in the web service/c++ code since we agree in the results of my java code. What you obtain from the web service is wrong, or more likely is calculated for a different observer location / date. You can try a different location with a well known lon/lat, and also a different c++ port available in a previous post above.
Here are the results using my JPARSEC library:
Sun Right ascension: 00h 40m 54.0270s Sun Declination: 04° 23' 51.458” Sun Distance: 0.999117 AU Sun Rise: 20190401 02:50:47 UT Sun Transit: 20190401 09:21:24 UT Sun Transit elevation: 39° 49' 32.3” Sun Set: 20190401 15:53:13 UT Sun Azimuth: 88° 57' 38.589” Sun Elevation: 4° 49' 28.781”
And here the results using this calculator, which are similar:
Sun Az: 88.964134° El: 4.8308864° Dist: 0.999117 AU RA: 10.225616° DEC: 4.396833° Ill: 100.0% ang.R: 0.26679978° Rise: 2019/04/01 02:50:46 UT Set: 2019/04/01 15:53:11 UT Transit: 2019/04/01 09:21:22 UT (elev. 39.824867°)
Input data:
I have a question, while everything appears to work correctly, I notice that the siderealtime does not match on here: http://www.neoprogrammics.com/moon/DE405_Lunar_Ephemeris_Calculator.php
However everything else does (within error)
Is there a difference between “apparent sidereal time” and “sidereal time”?, or is this an error?
Hi Sascha,
Right after entering the page you mention the default coordinates and UT time I can see are:
The results for the sidereal time on that page are:
Now on this program, the sidereal time is not directly accesible, you need to map it to another external variable or similar. One point here is that Sun/Moon position are recomputed to obtain the accurate rise/set times, and this will modify the sidereal time. In addition, the value is in radians without the 02PI reduction. Assuming you create the static variables globalGMST and globalLST to map the corresponding values (gmst and lst variables in doCalc method), at the end of the program you should repeat the position of a body to reset the sidereal time, something like:
My results are:
The difference between these values and those in the page you provide are 0.3 degrees or 1.2 minutes of time. Despite it seems little, the difference is too much, there is an error somewhere. If you go to the Horizons service of JPL at https://ssd.jpl.nasa.gov/horizons.cgi and take the time to input everything with the same values, you can check the AST provided there is 340.48913 degrees. Using my JPARSEC library, the value I obtain is 340.4899 degrees. There are different algorithms for the sidereal time and I know Horizons use an old one. Anyway, the error in the sidereal time provided in this calculator is in the 0.005 degrees level.
This difference of 1.2 minutes is suspicious, because the difference TTUT1 is 1.15 minutes currently. It seems the page you mention uses TT to compute the sidereal time instead of UT, so I think there is a bug on that page that indeed also affects the azimuth and elevation of the Moon reported there. Those values agree to better than 0.01 degrees level here and in Horizons, but on that page they are off by 0.2 degrees. In case that page is broadly used by people it would be good to report this bug there.
Regards
with this input:
SunMoonCalculator output:
How to force it to output details for only the input date (20190726), I mean:
At the end of the doCalc method you can read:
The rise and set times are calculated for the inmediate past and future, and one of those values is selected (the closest in time instead the one corresponding to today). This decision is justified because of the iteration required for accurate results in moving objects. You can change this behavior and select the one corresponding to today, although I think the convergency of the iteration will not be completely guarranteed, specially if there is no such result for today.
The easiest and probably better way is to repeat computations for the next and previous days (with the same time) and to check manually if there is a result for a given calendar day in any of the computations, with this strategy you can also account for possible time zone corrections to obtain a result, if any, for a given calenday day in local time and not UT.
To be more specific, it should work if you replace the if sentences with
You will have to compute riseToday1 and setToday1 also and check if jdToday is equal to riseToday1 or riseToday2, same for set. If not for both values you will have the no moon rise/set event (set rise/set variables to 1 for that). But remember this will only be valid for UT time scale.
How to calculate rise and set times taking in consideration the observers's elevation (altitude)? I think it is more accurate.
Being height the elevation in meters, compute the depresion of the horizon (<plus> is the plus symbol)
and in the twilight switch statement (doCalc method) include it, for instance for the standard 34' refraction at horizon:
I added this code just after the switch :
will this work even when tmp is positive?
Yes, but tmp is never positive unless you want to get when the object is at a given positive elevation above the effect of refraction.
Exactly! I edited the switch as follow:
need to allow the user to get geocentric details just by passing a parameter to:
so we get geocentric elevation and azimuth and other details, without setting radiusAU to 0 manually (I mean it will be done by the method if geocentric = true).
need also to implement methods to calculate crescent width and visibility using known creterion:
Sorry but I don't have time to help with that, if you need it then it is your homework.
How to get moon Geocentric Conjunction and Topocentric Conjunction times?
While there is only one Geocentric Conjunction per lunar month, will an observer have many Topocentric Conjunctions?
Once again, it is your homework, just use a searchlike iterative algorithm.
You are saying last logic seems working can you give the example please i am trying since last 3 days
nakshtras = (int) (normalizeRadians(longitude * DEG_TO_RAD) * 27 / TWO_PI);
its showing 28 but actual value of nakshatras if Aswini which is 1
To get a chance of an answer from Vikrant you should reply to him at the top. Anyway, the expression you post is a number between 0 and 26, you will only get a result >=28 if you add a constant. Of course, due to that constant in case the output is outside 127 you should reduce it to that range with +27 or 27, so just do nakshtras=nakshtras27 and you will get 1.
Hi! Thanks for this amazing resource, it calculates all the data I need and seems very accurate.
I've noticed that when the moon barely rises above the horizon, the rise and set events can be missed. This seems to happen because in the obtainAccurateRiseSetTransit method, the final value of step after iteration is slightly more than (1.0 / SECONDS_PER_DAY)  although it did converge.
For an example, see 4th June 2020 in Nordkapp, Norway  https://www.timeanddate.com/moon/@778460
int year = 2020, month = 6, day = 4, h = 0, m = 0, s = 0; double obsLon = 25.7813 * DEG_TO_RAD, obsLat = 71.1696 * DEG_TO_RAD;
Even though the apparent elevation reaches 0.19 degrees no rise or set is calculated. Sometimes, changing the date I started with by a few hours can result in a small enough step value that missed events are calculated, so I have ended up calculating rise and set for several times through the day to try to find all events.
Should I increase the iteration count, or is there a better way to resolve this?
You are right, there is a problem inherent to the iteration process used. For high latitudes the Moon or Sun could rise almost parallel to the horizon, and this requires more iterations since the exact second (within the accuracy limitations of the program) is hard to obtain. If you change the input time you could get the output if you are lucky, but the iteration is not guarranteed to reach one second accuracy. I have tested that 7 iterations are needed in the example you mentioned, and the program was using 5 for the Moon, and 3 for the Sun.
To try to solve this I have changed the code to 15 iterations for both bodies, that should be enough for all cases, but again I haven't tested this in depth. Since this would slow down the computation for most locations, I have also added a convergency check in obtainAccurateRiseSetTransit to leave the loop when enough accuracy is reached before all those iterations.
Thanks a lot for reporting this problem.
Thanks for your quick reply! Breaking when one second accuracy is reached is a good solution, and I may also reduce the accuracy requirement a bit to avoid excessive iterations. Until now my app has used 1 minute accuracy from a previous algorithm.
Hi Tomás, this is really good  it's the first library I've tried which has decent accuracy.
I'd like to port this to JavaScript for a project of mine  is there an opensource licence for the code you have shared here?
If not, would it be ok to release my port under the MIT licence? More info on licensing can be found at https://choosealicense.com/.
Thanks Josh. MIT license is ok, it is also my intention for this calculator, so my complete library (JPARSEC) continues as GPL. I would appreciate a link to your javascript code (in this discussion or sending me an email), I'm thinking about adding a little section in this blog entry about other ports.
2y I ported it to C for the Arduino/ESP8266. I gave credit in the code like this: https://github.com/ThingPulse/esp8266weatherstation/blob/master/src/SunMoonCalc.cpp#L23 My code is also MIT.
Today, June 10 2020, I have revised the SunMoonCalculator class and added another class to compute the positions of the planets. I did only basic tests but everything seems to work fine. I think is a nice addition that could be useful in other projects. Since some of you are subscribed to comments I think is a good idea to share/spread it here.
I also would like to say thanks to all of you for using the code in your own projects.
Hi!, a small contribution to this Project in Python language. Credits included.
GitHub respository: https://github.com/bokepasa/SunMoonCalculator
Thanks Tomás.
Not so small contribution,
¡Gracias Leandro!
Hi.
Love your code. I was attempting to implement the moon calcs from Meeus chapter 47 when I came across your project. It saved me so much writing/debugging/swearing. Your moon calcs differ from Meeus and your precision is truncated but as they say 'The test of the pudding is in the tasting' and your efforts product good numbers so who am I to argue.
While porting your code I came across a small typo/cutpaste error in the extended precision moon calculations. The following term was duplicated.
l += E * 1.22E3 * Math.sin(4*phasesanomalyanomaly)……..
l += E * 1.22E3 * Math.sin(4*phasesanomalyanomaly)……..
I do not think the lunar orbit will destabilize due to this error but.. better safe than sorry. …………..
When I am finished one other feature I will need it the ability to supply an objects RA/Dec and return the localized El/Az.. i.e where is Messier M31.
Any suggestions??\
Thanks for your efforts.
Ups! Yes, a cut/paste error because I did the code in two times and left that line two times… Thanks for the typo, maybe the lunar position will be now even better than 0.005 degrees for current dates.
Respect computing the position of stars or deep sky objects from catalog coordinates, you have to reduce them in a similar way respect the planets. Assuming you have J2000 equatorial coordinates (RA/DEC referred to J2000 equinox, which is the typical situation), you have to:
1. Compute ecliptic J2000 rectangular coordinates from RA/DEC using the obliquity at J2000.
2. Compute mean ecliptic coordinates of date, by calling method precessFromJ2000 listed in the section for computing the position of the planets. p = (x, y, z)
3. Convert p back again to spherical and reduce using doCalc. Use a high value for r and whatever for angRadius. Here I add the approximate correction for aberration before calling doCalc. The instance of the SunMoonCalculator class to be used to call doCalc should have as input date the final desired date for the computations, the same for step 2.
NOTE: haven't tested myself, but that's the process.
The latest and most complete code is available at a new repository for the Sun/Moon calculator: https://bitbucket.org/talonsoalbi/sunmooncalculator/commits/. The repository provides the history of all code changes with diff between versions, so it can be useful to update the code in your projects, especially for ports to other languages. For instance, the diff file between the July 3, 2020 and the August 7, 2020 versions would be this link.
Hi Tomás, thanks for providing the code! I have a question with regard to the TTminusUT calculation. I do not understand why TTminusUT (when going back in the past) decreases until about year 1800 and then starts increasing up to 19000sec for the year 600. I would expect a constantly decrease because earth's rotation speeds up when going back in time.
Hi Christoph, I'm not an expert on that so without reading the experts I can only suggest that the answer is probably the little amount of time in which the formula can be applied. After millions of years I would expect the same, but in very short time scales (few millenia is very short in this case) the variation is irregular in any direction and can only be tracked with historical records. I suppose a very longterm variation has been modeled/computed, but the additional shortterm variations on it are unpredictable. The same happens with the exact direction of the rotation axis in the range of ten meters, I think has to do with changes in the properties or mass distribution of the Earth mantle and nucleus.
Another similar case is the longterm climatic change versus shortterm climatology. Today we have the greatest snowstorm in recent history here in Madrid, but the trend after some decades towards a global warming is clear.
Thanks for the quick reply. Yes, climate can be crazy, the news today showed people skiing in Madrid City… I think I found the answer of my question. It has nothing to do with shortterm variations. TTminusUT or deltaT is always positive and the reason why it is increasing also when going back in time is the following: If in the past the moment was for example 12:00 noon TT, the time in UT was 11:00 (because earth was rotating faster at that time). DeltaT was positiv. An UT watch was faster than a TT watch. This difference however decreased with time as earth rotation and a such the UT time slowed down, until around 1800, where both times where synchronous (deltaT =0). With further time the UT watch more and more slows down compared to the TT watch. So sometimes in future there is a moment where it is 12:00 TT and 11:00UT. Again deltaT is positiv. If we consider only the long term change in earth rotation, the TTminusUT should follow a simple const*T^2 polynomial. Your polynomial is much more complicated as it considers also short term effects!
Thanks Christoph, I had that quite forgotten. The idea is that the deceleration of the rotation w is measured in rad s^2, similar to the linear acceleration of a body in m s^2. The radians can be transformed to seconds of time knowing the length of the day (86400s for 2PI) and the s^2 to centuries^2. To obtain deltaT in s you need to compute w * t * t + referenceConstant, and the approximate result is deltaT = 20 + 32 * t * t, with t = (year  1820) / 100 because the length of the day in 1820 was close to 86400s and was taken arbitrary as reference following some work of Newcomb (explained in section 4 at https://royalsocietypublishing.org/doi/10.1098/rspa.2016.0404). The 32 in the approximate deceleration without shortterm effects is caused by the Moon of course, and can change with time but is constant for only a few millenia, producing a simple parabola in deltaT. For nonhistorical times it can be modeled knowing the position of the Moon and has been done in the JPL ephemerides. The shortterm effects are caused by the mantle of the Earth (mass displacement, heat difussion), the ocean and climatology effects, and the magnetic field.
Thank you for the explanation!
One followup question on deltaT: I played around with your code to calculate sun eclipses in the very past millenniums. Of course the code is somewhat limited for instances which are far away in the past, however I could improve the accuracy by implementing a correction to deltaT (using deltaT = 20 + 32*t*t with t = (year  1820)/100 ). This correction is done similar to the explanation found here: https://eclipse.gsfc.nasa.gov/SEcat5/secular.html Now it seems that the moon algorithm provided above with less terms is interestingly predicting the eclipse time better than the more precise moon position calculation. The value 32 in the deltaT equation above was lower in both cases with about 30 for the moon routine with less terms and about 26.5 for the more precise routine. For the extreme example of a sun eclipse in Oct 04 1996 (https://eclipse.gsfc.nasa.gov/SEsearch/SEsearchmap.php?Ecl=19961004) I had to reduce deltaT = 46560 sec by 7500 sec for the more precise routine and by 2500 sec for the shorter routine. It seems that the updated moon algorithm is less accurate in predicting the secular moon acceleration as the correction in deltaT was larger. Can this be the case? Is there a way to improve the secular precision?
Sorry for the previous answer, it was completely wrong, I answered too fast without testing both methods and making again a similar error I did in my library years ago.
First, you are right, the method with less terms is more accurate in the ancient past. This is due to the expansions I used, from other sources, for the mean anomaly and mean longitude of the Moon and the Sun (secular terms). In the method with less terms I use expansions up to t^4, coming from the book Calendrical Calculations, which are similar to the ones published by Simon et al. 1994. But Peter DuffetSmith, in his MOON program, only uses terms up to t^3, that cannot be used in the ancient past without losing accuracy. In my library I didn't implement the MOON program, I just used it many years ago in Basic. In case I changed here his parameters with others because I tested the others were more adequate for ancient times, I cannot remember. Maybe I did that or I improved the code accidentally, to later go back to his original implementation that resulted worse for ancient dates. The original MOON program is only better for current dates. But as you noticed, to get good accuracy for ancient times the secular arguments computed at the beginning needs to be improved with even more terms.
The most complete expansions I'm aware of are those by S. L. Moshier. He fitted the mean elements to t^7 to fit the ELP2000 lunar ephemerides to DE404 and extend it to ancient dates. The rest of sources, including SOFA and others, do not use his work and only reach t^4. The accuracy when using all those terms up to t^7 are a lot better without adding many lines or computations, so I will update the code soon with both the TTUT and the expansion by Moshier.
Respect the 32, I would not change it to 'fit' a given eclipse. Testing the methods against an accurate lunar ephemerides shows that the difference can be sometimes very little, so you would get different values around 32 for different eclipses. The next version of the code should provide results much closer to the values published by Spenak for that eclipse in year 1996.
Thanks a lot for pointing me to this, my intention wasn't to make a good code for several millenia into the past, but I didn't realize how important was to add powers beyond t^4 to the secular arguments and how much could this improve the results in ancient times.
Thanks for your effort. I have seen that you have already updated the getMoon code accordingly. I tested it and it works! In addition to your code I included also the normalizeRadians to the secular terms, which was necessary in my case. From the figures above it seems that the precision in lon/lat is nearly independent of the selected millennium, is that correct? In the new code the correction for TTminusUT is applied also for calculation the sun position. Is this correct? However since the correction is very small it would not significantly affect the sun's position.
The new TTminusUT correction is for years before 600 and after 2200, and connects quite well with the polynomial fits. The little additional correction outside 19552005 is also documented in the Spenak web page. They are also correct for the Sun and planets, although the algorithm for the planets is probably very bad outside the valid time span of the method.
Respect the Moon, as you see the accuracy do not degrade easily with time respect DE404. This is thanks to the expansions by Moshier (see his work http://adsabs.harvard.edu/pdf/1992A%26A...262..613M and the cmoon and plan404 programs at http://moshier.net/). He first extended to 11000 years around the 2000 the ELP2000 theory (based on DE200 but using periodic terms) by fitting the coefficients to the DE200, finding new polynomial expansions for the secular terms, but later he revised that and did the same with DE404, extending also the polynomials to t^7, something never done before or after, up to my knowledge. The expansions were used for analytical theories and now numerical integration is used instead, especially for the ancient past. It is not clear if he did that fitting also to 11000 years into the past and future, but it seems so since there is no other good reason to include powers of 7. The moon program he did provides extremely accurate lunar positions only between 1370 B.C. and 3000 A.D., but, surprisingly, the expansions up to t^7 are extremely useful for reasonably good positions, even with few terms coming from other simple theories. Of course, later integrations consider more asteroids, better models, or a more accurate tidal acceleration, that will affect the position of the Moon for several millenia in the past, so a comparison with DE430 and a long term expansion fitted to DE430 would improve the ephemerides even more. A few tests with the Horizons server show discrepancies for the Moon around 150” respect DE430 for year 2000 (affecting only 10s in eclipses or to rise/set times), and about half of that for 1000. But note the previous discrepancy was 10x more. Accuracy means discrepancy respect reality, but the best we can do is to compare with the most recent JPL numerical integration.
In the previous version of the code or the DE200 integration, the mean values for the position of the Sun/Moon or secular terms start to deviate clearly from later works for calculations more than 500 years far from year 2000. When using the new secular terms it seems the shortterm periodic terms, like those from the complete MOON program, still contribute to improve the accuracy even several millenia far from 2000. But in this process, more terms can arise (Moshier also added more terms to ELP2000). The two terms I added have periods of 270 and 19 years, they could be related to objects from the Kuiper belt and the Earth's figure.
The lunar eclipse computations by Spenak use the ELP2000 theory and newer expansions for the mean elements provided by the authors of that theory (https://www.aanda.org/articles/aa/pdf/2002/20/aa2201.pdf). I assume they will be accurate up to year 2000, but they only reach t^4. It could be interesting to compare both expansions, for eclipses enough far those by Moshier are probably better. Anyway, the uncertainty in TTminusUT will affect more.
So the new code can be used probably up to 10 millenia into the past and future, but discrepancy with recent or future JPL integrations will increase with time, although not to very high levels.
Interesting work you have provided. That was a good idea to do tests with the Horizons server. When I did this I could identify several bugs in my code (floating point issues)! However, when everything worked, the discrepancy where surprisingly slightly different to your plots shown above. For the sun there was only a purely linear discrepancy which I could correct by one additional term
l += 3.6306/3600.0 * t;
in the getSun routine. For the date 1 000 000 JDTT for example this adds about 144” to lon! With this term the precision was similar to your plot. However I did not observe the decrease in precision towards year 4000.
For the modified getMoon I used the two additional periodic terms you have suggested. But in addition, a quadratic term
l += (0.07 * t * t + 0.3675 * t) / 3600.0;
was needed to improve accuracy in my case (to +/15” between 2000 to 4000). I don't know why this term is needed only in my case. I used the same modifications with all t^7 terms as you did. Could the moon's secular acceleration be the reason in the way, that the acceleration was different to the Horizons calculation?
I did not compare with DE430, the current ephemerides of Horizons, because I thought I had to download the entire dataset, although maybe the Horizons can directly provide thousands of points between 2000 and 4000. In addition, I had the expansions for the mean elements for DE404 so it was natural to compare with DE404 and maybe add some periodic terms to improve the ephemerides. For the Moon I used t but I could have used phase instead. For the Sun I didn't pay much attention to the problem around year >3000, since the periodic terms added where good enough for the rest of the interval. The point of using periodic terms is that any correction would not get bigger and bigger with time, and if that is necessary then the problem are the mean elements.
I suppose the models in DE404 and DE430 are similar, with minor modifications to the masses of the bodies, the initial conditions of the integration for positions and velocities, and the secular acceleration of the Moon. But these differences get bigger with time and add quadratic or other terms to the mean elements, which is what you are fitting. I suppose you found a dominant effect in the mean longitude of the Moon, but you should also check the ecliptic latitude because some effect may appear also in the longitude of the ascending node. Another possible difference is the precession. Since I had control over the Moshier algorithm I used the same precession algorithm in both calculations, for DE430 probably Horizons allows to get the ecliptic coordinates of the Moon for J2000 so you can apply the precession yourself. In case everything is ok and with only that change you get a similar fit for the three position parameters of the ephemerides to the average of thousands of points from DE430, welldone, since DE430 will be better for ancient times, but I would suggest not going too much beyond the fitting window you used, because any extrapolation could lead to considerable errors because of the longterm behavior of the mean elements.
Ok that explain why the discrepancies are slightly different in your case. Please ignore what I wrote about the linear correction for the sun. The reason was insufficient floating point accuracy in my code (sun_element array not declared ad double and insufficient dec places for Pi). After fixing the accuracy was in the arcsecs range.
Good to know, 3.6” per century for the longitude of the Sun was strange, that's why I mentioned the possible difference in the precesion. The Sun should be very similar in both integrations. There is still something strange that may deserve a close comparison with Horizons or a fix, because for the solar longitude this program diverges from DE200 with a parabola with a minimum at year 700. Something is wrong with the method published by Bretagnon in 1986, so the fix I suggest above is quite bad. I may try to use the mean solar longitude by Moshier with a previous algorithm to see if it gets any better, but I agree a comparison with DE430 is the way to go.