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

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

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.02 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 10s or better for the Moon.

To obtain the difference TT-UT1 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. The effect of this correction in the output values is not very relevant, since the main source of uncertainty is the reduced number of terms in the calculation of the ecliptic coordinates of the Moon. Anyway, I have included enough terms to reach an accuracy in the arcminute range over several centuries, 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. It would be easy to do that with few lines of code, but the point is that this correction is not going to provide a better accuracy since the approximations used within this code produce a greater level of errors in the output positions, typically about a hundred km for the Moon and more for the other bodies.

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 final coordinates (azimuth and elevation, although RA and 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.

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 UT. The method calcSunAndMoon is used to compute everything at once comfortably, you maybe only need to call the Sun or the Moon (doCalc(getSun()) for the Sun) but an 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 several 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.

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. 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 Jean-Louis Simon, Willman-Bell, 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 original.

SunMoonCalculator.java
/**
 * A very simple yet accurate Sun/Moon calculator without using JPARSEC library.
 * @author T. Alonso Albi - OAN (Spain), email t.alonso@oan.es
 * @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 2000-01-01. */
	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. */
	protected double jd_UT = 0;
	protected double t = 0;
	private double obsLon = 0, obsLat = 0;
	private double TTminusUT = 0;
	private TWILIGHT twilight = TWILIGHT.HORIZON_34arcmin;
	private double nutLon = 0, nutObl = 0;
	protected double meanObliquity = 0;
 
	private 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);
 
		TTminusUT = 0;
		if (year > -600 && year < 2200) {
			double x = year + (month - 1 + day / 30.0) / 12.0;
			double x2 = x * x, x3 = x2 * x, x4 = x3 * x;
			if (year < 1600) {
				TTminusUT = 10535.328003 - 9.9952386275 * x + 0.00306730763 * x2 - 7.7634069836E-6 * x3 + 3.1331045394E-9 * x4 + 
						8.2255308544E-12 * x2 * x3 - 7.4861647156E-15 * x4 * x2 + 1.936246155E-18 * x4 * x3 - 8.4892249378E-23 * x4 * x4;
			} else {
				TTminusUT = -1027175.34776 + 2523.2566254 * x - 1.8856868491 * x2 + 5.8692462279E-5 * x3 + 3.3379295816E-7 * x4 + 
						1.7758961671E-10 * x2 * x3 - 2.7889902806E-13 * x2 * x4 + 1.0224295822E-16 * x3 * x4 - 1.2528102371E-20 * x4 * x4;
			}
		}
		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.2e-6 * T0) + 9.3104e-2) * T0) + 8640184.812866) * T0) + 24110.54841;
		double msday = 1.0 + (((((-1.86e-5 * 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 Jean-Louis 
	// Simon, Willman-Bell, 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() {
		// MOON PARAMETERS (Formulae from "Calendrical Calculations")
		double phase = normalizeRadians((297.8502042 + 445267.1115168 * t - 0.00163 * t * t + t * t * t / 538841 - t * t * t * t / 65194000) * DEG_TO_RAD);
 
		// Anomalistic phase
		double anomaly = (134.9634114 + 477198.8676313 * t + .008997 * t * t + t * t * t / 69699 - t * t * t * t / 14712000);
		anomaly = anomaly * DEG_TO_RAD;
 
		// Degrees from ascending node
		double node = (93.2720993 + 483202.0175273 * t - 0.0034029 * t * t - t * t * t / 3526000 + t * t * t * t / 863310000);
		node = node * DEG_TO_RAD;
 
		double E = 1.0 - (.002495 + 7.52E-06 * (t + 1.0)) * (t + 1.0);
 
		// Solar anomaly
		double sanomaly = (357.5291 + 35999.0503 * t - .0001559 * t * t - 4.8E-07 * t * t * t) * DEG_TO_RAD;
 
		// 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
		double l = (218.31664563 + 481267.8811958 * t - .00146639 * t * t + t * t * t / 540135.03 - t * t * t * t / 65193770.4);
		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.695E-3 * Math.sin(2 * anomaly - sanomaly) + 2.602E-3 * Math.sin(anomaly - 2*(node+phase));
		l += E * 2.396E-3 * Math.sin(2*(phase - anomaly) - sanomaly) - 2.349E-3 * Math.sin(anomaly+phase);
		l += E * E * 2.249E-3 * Math.sin(2*(phase-sanomaly)) - E * 2.125E-3 * Math.sin(2*anomaly+sanomaly);
		l += -E * E * 2.079E-3 * Math.sin(2*sanomaly) + E * E * 2.059E-3 * Math.sin(2*(phase-sanomaly)-anomaly);
		l += -1.773E-3 * Math.sin(anomaly+2*(phase-node)) - 1.595E-3 * Math.sin(2*(node+phase));
		l += E * 1.22E-3 * Math.sin(4*phase-sanomaly-anomaly) - 1.11E-3 * 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.73E-4 * Math.cos(3 * anomaly) + 1.67E-4 * Math.cos(4*phase-anomaly);
 
		// 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.472E-3 * Math.sin(2 * phase + node - sanomaly - anomaly);
		l += E * 2.222E-3 * Math.sin(2 * phase + node - sanomaly);
		l += E * 2.072E-3 * 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 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
		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://icts-yebes.oan.es/reports/doc/IT-OAN-2003-2.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;
 
		if (jd < 2299160.0 && jd >= 2299150.0)
			throw new Exception("invalid julian day " + jd + ". This date does not exist.");
 
		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 {
		if (jd < 2299160.0 && jd >= 2299150.0)
			throw new Exception("invalid julian day " + jd + ". This date does not exist.");
 
		// 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();
		}
	}
}

Extending the program for drawing purposes

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 (p-par) 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};
}

Demostration program for Java desktop

The previous code to compute the position of the Moon and the Sun and the orientation angles of the Moon has been recently 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.


 

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.

Demostration program for iOS (IPhone)

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.

Ports available of the Sun/Moon calculator to other languages

Besides the original Java version I created there are some ports to other languages I would like to emphasize. The Swift and C++ ports are from old versions of the calculator while those for Kotlin and Python are recent.

Swift (iOS):

https://github.com/kanchudeep/SunMoonCalculator

C++ (for Arduino):

https://github.com/ThingPulse/esp8266-weather-station/blob/master/src/SunMoonCalc.cpp

Kotlin (includes planetary positions using a different code from the one presented below):

https://github.com/davemorrissey/sundroid/tree/master/app/src/main/kotlin/uk/co/sundroid/util/astro/smc

Python:

https://github.com/bokepasa/SunMoonCalculator

Computing some events using iterative algorithms

There are many different kind of computations that can be solved using iterative and search-like algorithms. For Android-like 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 thread-safe, 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 += 1-age;
		}
		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;
}

More accuracy in the position of the Moon ?

The implementation presented above for the position of the Moon can be improved by adding all the terms from the MOON program by Peter Duffet-Smith. In addition, in that case it is convenient to implement the theory in the same way he did instead of obtaining the mean elements from a different source (something I did because it required less lines of code). Here we have the problem of accuracy versus speed, if you prefer accuracy you will 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.

protected double[] getMoon() {
	// Implementation following P. Duffet's MOON program
	double td = t + 1, td2 = t * t;
	double qd = td * JULIAN_DAYS_PER_CENTURY * 360.0;
 
	double M1 = qd / 2.732158213E1, M2 = qd / 3.652596407E2, M3 = qd / 2.755455094E1;
	double M4 = qd / 2.953058868E1, M5 = qd / 2.721222039E1, M6 = qd / 6.798363307E3;
 
	double l = normalizeRadians((2.70434164E2 + M1 - (1.133E-3 - 1.9E-6 * td) * td2) * DEG_TO_RAD);
	double sanomaly = normalizeRadians((3.58475833E2 + M2 - (1.5E-4 + 3.3E-6 * td) * td2) * DEG_TO_RAD);
	double anomaly = normalizeRadians((2.96104608E2 + M3 + (9.192E-3 + 1.44E-5 * td) * td2) * DEG_TO_RAD);
	double phase = normalizeRadians((3.50737486E2 + M4 - (1.436E-3 - 1.9E-6 * td) * td2) * DEG_TO_RAD);
	double node = normalizeRadians((11.250889 + M5 - (3.211E-3 + 3E-7 * td) * td2) * DEG_TO_RAD);
	double NA = normalizeRadians((2.59183275E2 - M6 + (2.078E-3 + 2.2E-6 * td) * td2) * DEG_TO_RAD);
	double A = DEG_TO_RAD * (51.2 + 20.2 * td);
	double S1 = Math.sin(A), S2 = Math.sin(NA);
	double B = 346.56 + (132.87 - 9.1731E-3 * td) * td;
	double S3 = 3.964E-3 * Math.sin(DEG_TO_RAD * B);
	double C = NA + DEG_TO_RAD * (275.05 - 2.3 * td);
	double S4 = Math.sin(C);
	l = l * RAD_TO_DEG + (2.33E-4 * S1 + S3 + 1.964E-3 * S2);
	sanomaly = sanomaly - (1.778E-3 * S1) * DEG_TO_RAD;
	anomaly = anomaly + (8.17E-4 * S1 + S3 + 2.541E-3 * S2) * DEG_TO_RAD;
	node = node + (S3 - 2.4691E-2 * S2 - 4.328E-3 * S4) * DEG_TO_RAD;
	phase = phase + (2.011E-3 * S1 + S3 + 1.964E-3 * S2) * DEG_TO_RAD;
	double E = 1 - (2.495E-3 + 7.52E-6 * td) * td, E2 = E * E;
 
	// Now longitude, with the three main correcting terms of evection,
	// variation, and equation of year, plus other terms (error<0.01 deg)
	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.695E-3 * Math.sin(2 * anomaly - sanomaly) + 2.602E-3 * Math.sin(anomaly - 2*(node+phase));
	l += E * 2.396E-3 * Math.sin(2*(phase - anomaly) - sanomaly) - 2.349E-3 * Math.sin(anomaly+phase);
	l += E * E * 2.249E-3 * Math.sin(2*(phase-sanomaly)) - E * 2.125E-3 * Math.sin(2*anomaly+sanomaly);
	l += -E * E * 2.079E-3 * Math.sin(2*sanomaly) + E * E * 2.059E-3 * Math.sin(2*(phase-sanomaly)-anomaly);
	l += -1.773E-3 * Math.sin(anomaly+2*(phase-node)) - 1.595E-3 * Math.sin(2*(node+phase));
	l += E * 1.22E-3 * Math.sin(4*phase-sanomaly-anomaly) - 1.11E-3 * Math.sin(2*(anomaly+node));
	l += E * 1.22E-3 * Math.sin(4*phase-sanomaly-anomaly) - 1.11E-3 * Math.sin(2*(anomaly+node));
	l += 8.92E-4 * Math.sin(anomaly - 3 * phase) - E * 8.11E-4 * Math.sin(sanomaly + anomaly + 2 * phase);
	l += E * 7.61E-4 * Math.sin(4 * phase - sanomaly - 2 * anomaly);
	l += E2 * 7.04E-4 * Math.sin(anomaly - 2 * (sanomaly + phase));
	l += E * 6.93E-4 * Math.sin(sanomaly - 2 * (anomaly - phase));
	l += E * 5.98E-4 * Math.sin(2 * (phase - node) - sanomaly);
	l += 5.5E-4 * Math.sin(anomaly + 4 * phase) + 5.38E-4 * Math.sin(4 * anomaly);
	l += E * 5.21E-4 * Math.sin(4 * phase - sanomaly) + 4.86E-4 * Math.sin(2 * anomaly - phase);
	l += E2 * 7.17E-4 * 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.73E-4 * Math.cos(3 * anomaly) + 1.67E-4 * Math.cos(4*phase-anomaly);
	p += -E * 1.11E-4 * Math.cos(sanomaly) + 1.03E-4 * Math.cos(4 * phase - 2 * anomaly);
	p += -8.4E-5 * Math.cos(2 * anomaly - 2 * phase) - E * 8.3E-5 * Math.cos(2 * phase + sanomaly);
	p += 7.9E-5 * Math.cos(2 * phase + 2 * anomaly) + 7.2E-5 * Math.cos(4 * phase);
	p += E * 6.4E-5 * Math.cos(2 * phase - sanomaly + anomaly)
			- E * 6.3E-5 * Math.cos(2 * phase + sanomaly - anomaly);
	p += E * 4.1E-5 * Math.cos(sanomaly + phase) + E * 3.5E-5 * Math.cos(2 * anomaly - sanomaly);
	p += -3.3E-5 * Math.cos(3 * anomaly - 2 * phase) - 3E-5 * Math.cos(anomaly + phase);
	p += -2.9E-5 * Math.cos(2 * (node - phase)) - E * 2.9E-5 * Math.cos(2 * anomaly + sanomaly);
	p += E2 * 2.6E-5 * Math.cos(2 * (phase - sanomaly)) - 2.3E-5 * Math.cos(2 * (node - phase) + anomaly);
	p += E * 1.9E-5 * 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.472E-3 * Math.sin(2 * phase + node - sanomaly - anomaly);
	l += E * 2.222E-3 * Math.sin(2 * phase + node - sanomaly);
	l += E * 2.072E-3 * Math.sin(2 * phase - node - sanomaly - anomaly);
	l += E * 1.877E-3 * Math.sin(node - sanomaly + anomaly) + 1.828E-3 * Math.sin(4 * phase - node - anomaly);
	l += -E * 1.803E-3 * Math.sin(node + sanomaly) - 1.75E-3 * Math.sin(3 * node);
	l += E * 1.57E-3 * Math.sin(anomaly - sanomaly - node) - 1.487E-3 * Math.sin(node + phase);
	l += -E * 1.481E-3 * Math.sin(node + sanomaly + anomaly) + E * 1.417E-3 * Math.sin(node - sanomaly - anomaly);
	l += E * 1.35E-3 * Math.sin(node - sanomaly) + 1.33E-3 * Math.sin(node - phase);
	l += 1.106E-3 * Math.sin(node + 3 * anomaly) + 1.02E-3 * Math.sin(4 * phase - node);
	l += 8.33E-4 * Math.sin(node + 4 * phase - anomaly) + 7.81E-4 * Math.sin(anomaly - 3 * node);
	l += 6.7E-4 * Math.sin(node + 4 * phase - 2 * anomaly) + 6.06E-4 * Math.sin(2 * phase - 3 * node);
	l += 5.97E-4 * Math.sin(2 * (phase + anomaly) - node);
	l += E * 4.92E-4 * Math.sin(2 * phase + anomaly - sanomaly - node)
			+ 4.5E-4 * Math.sin(2 * (anomaly - phase) - node);
	l += 4.39E-4 * Math.sin(3 * anomaly - node) + 4.23E-4 * Math.sin(node + 2 * (phase + anomaly));
	l += 4.22E-4 * Math.sin(2 * phase - node - 3 * anomaly)
			- E * 3.67E-4 * Math.sin(sanomaly + node + 2 * phase - anomaly);
	l += -E * 3.53E-4 * Math.sin(sanomaly + node + 2 * phase) + 3.31E-4 * Math.sin(node + 4 * phase);
	l += E * 3.17E-4 * Math.sin(2 * phase + node - sanomaly + anomaly);
	l += E2 * 3.06E-4 * Math.sin(2 * (phase - sanomaly) - node) - 2.83E-4 * Math.sin(anomaly + 3 * node);
	double W1 = 4.664E-4 * Math.cos(NA);
	double W2 = 7.54E-5 * 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))};
}

The previous code is accurate enough to start thinking about implementing a better 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 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.03 to 0.005 degrees (just a few arcseconds) in current dates, which means around 1-2 seconds in the rise, set, and transit times. Moon distance will be accurate to 30 km or better. Even around year 1000 A.D. the accuracy will be still around 0.1 degrees, well below the apparent radius of the Moon.

Extending the program to compute planetary positions

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 1000-3000 A.D., but will give better results between 1800-2100. 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 Earth-Moon 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 light-time 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.

PlanetCalculator.java
/**
 * 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       semi-major 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,         3e-10,     0.0 },  /* Mars    */
	       new double[] {  5.2026032092,     19132e-10, -39e-10 },  /* Jupiter */
	       new double[] {  9.5549091915, -0.0000213896, 444e-10 },  /* Saturn  */
	       new double[] { 19.2184460618,     -3716e-10, 979e-10 },  /* Uranus  */
	       new double[] { 30.1103868694,    -16635e-10, 686e-10 }   /* 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,    -28349e-10 },
			new double[] { 0.0067719164, -0.0004776521,     98127e-10 },
			new double[] { 0.0167086342, -0.0004203654, -0.0000126734 },
			new double[] { 0.0934006477,  0.0009048438,    -80641e-10 },
			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 semi-major 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 light-time 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., Chapront-Touze, M., Francou, G., and 
	 * Laskar, J., Astronomy and Astrophysics 282, 663 (1994).
	 *  
	 * Expected maximum error over the interval 1800-2100 A.D.: 
	 * better than 110" for Jupiter, Saturn, and Uranus. 40" for Mars, <=20" for the others.
	 * Over 1000-3000 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)) * 1e-7;
			dl += (cl[np][k] * Math.cos(argl) + sl[np][k] * Math.sin(argl)) * 1e-7;
		}
		arga = kp[np][8] * dmu;
		da += t * (ca[np][8] * Math.cos(arga) + sa[np][8] * Math.sin(arga)) * 1e-7;
		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)) * 1e-7;
		}
		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) > 1e-12) {
			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 1000-3000, but 1800-2100 is recommended
			// Expected accuracy (maximum errors) 1800-2100:
			// 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();
		}
	}
}

Overall the code presented here 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 millenia, 0.005 for the Moon in current dates (0.1 degrees 1000 years far from year 2000, probably around 1-2 degrees 2000 years far), and for the planets 0.05 degrees or better up to 1000 year fars 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.

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 typical accuracies 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 millenia. The method is less suitable for low powerful devices but extremely more accurate than this post. Have a look at the final section of that post. In case of using the methods described in this post the accuracy is limited to 1”.

Discussion

vikrant bhave, 2013/12/12 11:39

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

Tomás Alonso Albi, 2013/12/13 11:00

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.

vikrant bhave, 2013/12/20 15:26

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.

Tomás Alonso Albi, 2013/12/20 16:59

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))

vikrant bhave, 2013/12/23 13:16

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.

Bryan, 2016/04/24 05:06

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.

Tomás Alonso Albi, 2016/04/25 18:19

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.

Bryan, 2016/04/26 00:30

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 moon-rise that day I will display something like “none” or “no moon-rise 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).

Tomás Alonso Albi, 2016/04/26 17:46

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.

Bryan, 2016/04/27 01:37

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)?

Bryan, 2016/04/27 01:43

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?

Tomás Alonso Albi, 2016/04/27 18:04

- 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.

Bryan, 2016/04/28 03:33

– 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.

Bryan, 2016/04/28 03:45

Just FYI I am allowing the user to use a date picker to choose dates in the future and past.

Tomás Alonso Albi, 2016/04/28 15:48

← 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.

Bryan, 2016/04/29 05:03

–> 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.

Tomás Alonso Albi, 2016/04/29 15:24

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.

Bryan, 2016/05/03 03:09

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 {

  static final double RAD_TO_DEG = 180.0 / Math.PI;
  /** Degrees to radians. */
  static final double DEG_TO_RAD = 1.0 / RAD_TO_DEG;
  //  calendar and date formatter objects
  static Calendar calendar = Calendar.getInstance();
  static long JD = calendar.getTimeInMillis();
  static int year = calendar.get(Calendar.YEAR);
  static int calMonth = calendar.get(Calendar.MONTH); // first month in calendar is 0 so add 1
  static int month = calMonth + 1;
  static int day = calendar.get(Calendar.DAY_OF_MONTH); // day number (first day of month is 1, second is 2, etc.).
  static int h = calendar.get(Calendar.HOUR); // hours
  static int m = calendar.get(Calendar.MINUTE); // minutes
  static int s = calendar.get(Calendar.SECOND); // seconds
  static int timeZone = (calendar.get(Calendar.ZONE_OFFSET) / (1000 * 60 * 60) % 24); // time zone offset
  static int daylightSavings = (calendar.get(Calendar.DST_OFFSET) / (1000 * 60 * 60) % 24);
  //int am_pm = calendar.get(Calendar.AM_PM); // am / pm indicator
  static double obsLon = -94.56 * DEG_TO_RAD;  // observer longitude
  static double obsLat = 40.88 * DEG_TO_RAD; // observer latitude

… the rest of your setup code if I left any thing out

  public void calculate() { // < this is just your Main method trimmed down to give only      the data I want to display
      //  get references to programaticlly manipulated TextViews
      sunRiseText = (TextView) findViewById(R.id.sunRiseText);
      

… more of the same as above

      try {
          SunMoonCalculator smc = new SunMoonCalculator(year, month, day, h, m, s, obsLon, obsLat);
          smc.setTwilight(SunMoonCalculator.TWILIGHT.TWILIGHT_CIVIL);
          smc.calcSunAndMoon();

… other unrelated code that uses moonAge to change the picture displayed to the correct

    moon phase picture, and display the approprite text (full moon, first crescent, etc...)
          //sunrise   < my attempt at using your 3 tries solution, it does not work for me
          int date1[];
          int date2[];
          int date3[];
          int sunriseTime[] = {0, 0, 0};
          date1 = SunMoonCalculator.getDate(smc.sunRise);
          date2 = SunMoonCalculator.getDate(smc.sunRise + 24);
          date3 = SunMoonCalculator.getDate(smc.sunRise - 24);
          if (date1[2] == day) {
              sunriseTime[0]= date1[3];
              sunriseTime[1]= date1[4];
              sunriseTime[2]= date1[5];
          }
          if (date2[2] == day) {
              sunriseTime[0]= date2[3];
              sunriseTime[1]= date2[4];
              sunriseTime[2]= date2[5];
          }
          if (date3[2] == day) {
              sunriseTime[0]= date3[3];
              sunriseTime[1]= date3[4];
              sunriseTime[2]= date3[5];
          } // < something I tried before that did not work
          //int sunriseTime[] = SunMoonCalculator.getDate(smc.sunRise + timeZone);
          // the rest of the 3 tries method
          int sunrisehour = sunriseTime[0];
          int sunriseminutes = sunriseTime[1];
          int sunriseseconds = sunriseTime[2];

some other things I have tried sunset

          int sunsetTime[] = SunMoonCalculator.getDate(smc.sunSet);
          int sunsethour = sunsetTime[3];
          int sunsetminutes = sunsetTime[4];
          int sunsetseconds = sunsetTime[5];
          //suntransit
          int suntransitTime[] = SunMoonCalculator.getDate(smc.sunTransit + timeZone);
          int suntransithour = suntransitTime[3];
          int suntransitminutes = suntransitTime[4];
          int suntransitseconds = suntransitTime[5];
          //moonrise
          int moonriseTime[] = SunMoonCalculator.getDate(smc.moonRise + timeZone);
          int moonrisehour = moonriseTime[3]; // + timeZone;
          int moonriseminutes = moonriseTime[4];
          int moonriseseconds = moonriseTime[5];
          // moonset
          int moonsetTime[] = SunMoonCalculator.getDate(smc.moonSet + timeZone / 24);
          int moonsethour = moonsetTime[3];
          int moonsetminutes = moonsetTime[4];
          int moonsetseconds = moonsetTime[5];
Bryan, 2016/05/03 03:15

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.

Bryan, 2016/05/03 03:41

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.

Tomás Alonso Albi, 2016/05/03 12:37

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[]) {

try {
	// Check time zone is in hours
	int timeZone = -5 + 1; // (+1 since it is summer)
	// Check these values are in radians, they are for Washington DC
	double obsLon = -77.03 * DEG_TO_RAD, obsLat = 38.9 * DEG_TO_RAD;
	// Set to false to return results for a different day (if local hour is close to 0 or 24)
	boolean forceResultsForToday = true;
	// Test date in UT, h (UT) = 20 (LT) - TZ
	int year = 2016, month = 5, day = 3, h = 20 - timeZone, m = 0, s = 0;
	
	// For today, h = 20 LT and not 12 to return here results for day 4, not 3
	SunMoonCalculator smc = new SunMoonCalculator(year, month, day, h, m, s, obsLon, obsLat);
	smc.calcSunAndMoon();
	// Save results in LT
	double sunRise = smc.sunRise + timeZone / 24.0;
	double sunSet = smc.sunSet + timeZone / 24.0;
	double sunTransit = smc.sunTransit + timeZone / 24.0;
	double moonRise = smc.moonRise + timeZone / 24.0;
	double moonSet = smc.moonSet + timeZone / 24.0;
	double moonTransit = smc.moonTransit + timeZone / 24.0;
	if (forceResultsForToday) {
		// Next day
		SunMoonCalculator smcTomo = new SunMoonCalculator(year, month, day+1, h, m, s, obsLon, obsLat);
		smcTomo.calcSunAndMoon();
		// Previous day
		SunMoonCalculator smcYest = new SunMoonCalculator(year, month, day-1, h, m, s, obsLon, obsLat);
		smcYest.calcSunAndMoon();
		// Select desired result for today
		if (getDate(sunRise)[2] != day) sunRise = smcTomo.sunRise + timeZone / 24.0;
		if (getDate(sunRise)[2] != day) sunRise = smcYest.sunRise + timeZone / 24.0;
		if (getDate(sunSet)[2] != day) sunSet = smcTomo.sunSet + timeZone / 24.0;
		if (getDate(sunSet)[2] != day) sunSet = smcYest.sunSet + timeZone / 24.0;
		if (getDate(sunTransit)[2] != day) sunTransit = smcTomo.sunTransit + timeZone / 24.0;
		if (getDate(sunTransit)[2] != day) sunTransit = smcYest.sunTransit + timeZone / 24.0;
		if (getDate(moonRise)[2] != day) moonRise = smcTomo.moonRise + timeZone / 24.0;
		if (getDate(moonRise)[2] != day) moonRise = smcYest.moonRise + timeZone / 24.0;
		if (getDate(moonSet)[2] != day) moonSet = smcTomo.moonSet + timeZone / 24.0;
		if (getDate(moonSet)[2] != day) moonSet = smcYest.moonSet + timeZone / 24.0;
		if (getDate(moonTransit)[2] != day) moonTransit = smcTomo.moonTransit + timeZone / 24.0;
		if (getDate(moonTransit)[2] != day) moonTransit = smcYest.moonTransit + timeZone / 24.0;
	}
	
	// Output removing 'UT'
	String o = "";
	System.out.println("Sun");
	o = SunMoonCalculator.getDateAsString(sunRise);
	if (o.endsWith(" UT")) o = o.substring(0, o.length()-3);
	System.out.println(" Rise:    "+o);
	o = SunMoonCalculator.getDateAsString(sunTransit);
	if (o.endsWith(" UT")) o = o.substring(0, o.length()-3);
	System.out.println(" Transit: "+o);
	o = SunMoonCalculator.getDateAsString(sunSet);
	if (o.endsWith(" UT")) o = o.substring(0, o.length()-3);
	System.out.println(" Set:     "+o);
	System.out.println("Moon");
	o = SunMoonCalculator.getDateAsString(moonRise);
	if (o.endsWith(" UT")) o = o.substring(0, o.length()-3);
	System.out.println(" Rise:    "+o);
	o = SunMoonCalculator.getDateAsString(moonTransit);
	if (o.endsWith(" UT")) o = o.substring(0, o.length()-3);
	System.out.println(" Transit: "+o);
	o = SunMoonCalculator.getDateAsString(moonSet);
	if (o.endsWith(" UT")) o = o.substring(0, o.length()-3);
	System.out.println(" Set:     "+o);
} catch (Exception exc) {
	exc.printStackTrace();
}

}

Tomás Alonso Albi, 2016/05/03 20:20

At the end of the forceResultsForToday block I would add this to also solve the no Moon event problem once a month

if (getDate(moonRise)[2] != day) moonRise = -1;
if (getDate(moonSet)[2] != day) moonSet = -1;
if (getDate(moonTransit)[2] != day) moonTransit = -1;
Bryan, 2016/05/04 01:30

Ok thanks, I will give it a try and let you know.

Bryan, 2016/05/05 04:27

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.

Bryan, 2016/05/05 05:16

I test against: http://aa.usno.navy.mil/data/docs/RS_OneDay.php - the U.S. Naval Observatory

Bryan, 2016/05/05 06:44

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

        rise          transit          set

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)

Tomás Alonso Albi, 2016/05/06 11:43

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: 1-2s with Sun, around 10s with Moon.

You mention you modified something to simplify the program. Check your program, you must be doing something wrong.

Bryan, 2016/05/06 21:52

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:

          double ageOfMoon = moonAge;
          double moonCycle = 29.530588853;
          double moonPhaseLength = (moonCycle / 8);
          double newMoonStart = moonCycle - (moonPhaseLength / 2);
          double newMoonEnd = moonPhaseLength / 2;
          double waxingCrescentStart = newMoonEnd;
          double waxingCrescentEnd = waxingCrescentStart + moonPhaseLength;
          double firstQuarterStart = waxingCrescentEnd;
          double firstQuarterEnd = firstQuarterStart + moonPhaseLength;
          double waxingGibbousStart = firstQuarterEnd;
          double waxingGibbousEnd = waxingGibbousStart + moonPhaseLength;
          double fullMoonStart = waxingGibbousEnd;
          double fullMoonEnd = fullMoonStart + moonPhaseLength;
          double waningGibbousStart = fullMoonEnd;
          double waningGibbousEnd = waningGibbousStart + moonPhaseLength;
          double lastQuarterStart = waningGibbousEnd;
          double lastQuarterend = lastQuarterStart + moonPhaseLength;
          double waningCrescentStart = lastQuarterend;
          double waningCrescentend =newMoonStart;
          String moonPhaseString = "None"; // for trouble shooting, if we get "none" it never calculated
          ImageView imageView = new ImageView(this);
          if ((ageOfMoon >= newMoonStart) && (ageOfMoon <= newMoonEnd)) {
              moonPhaseString = "New Moon";
              ImageView img= (ImageView) findViewById(R.id.phaseImgView);
              img.setImageResource(R.drawable.moon_new);
          }
          if ((ageOfMoon >= waxingCrescentStart) && (ageOfMoon <= waxingCrescentEnd)) {
              moonPhaseString = "Waxing Crescent";
              ImageView img= (ImageView) findViewById(R.id.phaseImgView);
              img.setImageResource(R.drawable.moon_waxcres); // imageView.setImageResource(@drawable/moon_waxcres);
          }
          if ((ageOfMoon >= firstQuarterStart) && (ageOfMoon <= waxingCrescentEnd)) {
              moonPhaseString = "First Quarter";
              ImageView img= (ImageView) findViewById(R.id.phaseImgView);
              img.setImageResource(R.drawable.moon_firstqtr);
          }
          if ((ageOfMoon >= waxingGibbousStart) && (ageOfMoon <= waxingCrescentEnd)) {
              moonPhaseString = "Waxing Gibbous";
              ImageView img= (ImageView) findViewById(R.id.phaseImgView);
              img.setImageResource(R.drawable.moon_waxgib);
          }
          if ((ageOfMoon >= fullMoonStart) && (ageOfMoon <= fullMoonEnd)) {
              moonPhaseString = "Full Moon";
              ImageView img= (ImageView) findViewById(R.id.phaseImgView);
              img.setImageResource(R.drawable.moon_full);
          }
          if (ageOfMoon >= waningGibbousStart && ageOfMoon <= waningGibbousEnd) {
              moonPhaseString = "Waning Gibbous";
              ImageView img= (ImageView) findViewById(R.id.phaseImgView);
              img.setImageResource(R.drawable.moon_wangib);
          }
          if (ageOfMoon >= lastQuarterStart && ageOfMoon <= lastQuarterend) {
              moonPhaseString = "Last Quarter";
              ImageView img= (ImageView) findViewById(R.id.phaseImgView);
              img.setImageResource(R.drawable.moon_lastqtr);
          }
          if (ageOfMoon >= waningCrescentStart && ageOfMoon <= waningCrescentend) {
              moonPhaseString = "Waning Crescent";
              ImageView img= (ImageView) findViewById(R.id.phaseImgView);
              img.setImageResource(R.drawable.moon_wancres);
          }
          //else {
          //    moonPhaseString = "Error"; // for trouble shooting, if we get "error" it calculated but was not caught by one of the if / if elses's
          //    ImageView img = (ImageView) findViewById(R.id.phaseImgView);
          //    img.setImageResource(R.drawable.error);
          //}

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:

  /**
   * Returns a time as a string.
   * @param jd The Juliand day.
   * @return The String.
   * @throws Exception If the date does not exists.
   */
  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 date[3]+":"+date[4]+":"+date[5];
  }

I think it is the problem with the incorrect month, but I have not had a chance to test it yet. Thank you.

Bryan, 2016/05/06 23:28

The sunrise and moon-rise 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 moon-rise 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.

Bryan, 2016/05/07 00:08

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, moon-rise 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.

Bryan, 2016/05/10 06:15

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 moon-rise 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).

Tomás Alonso Albi, 2016/05/10 13:16

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!

Bryan, 2016/05/11 02:24

Except for the moon-age 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 moon-age 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]);”

Tomás Alonso Albi, 2016/05/11 13:34

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.

Bryan, 2016/05/11 20:42

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.

Bryan, 2016/05/11 21:59

If I understand correctly, since the moon-age is just a count of days since the last new moon, no timezone correction is necessary, is this correct?.

Bryan, 2016/05/12 04:40

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 re-post 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.

Bryan, 2016/05/12 05:37

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.

Tomás Alonso Albi, 2016/05/12 11:34

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.

Bryan, 2016/05/12 18:08

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?

Bryan, 2016/05/12 19:44

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.

Bryan, 2016/05/12 20:09

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.

Bryan, 2016/05/12 20:49

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.

Bryan, 2016/05/23 18:54

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.

Tim Robinson, 2017/09/03 00:32

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.

Tomás Alonso Albi, 2017/09/05 10:29

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

Tim Robinson , 2017/09/05 19:29

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.

Tim Robinson , 2017/09/05 19:34

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.

Bruce Jin, 2017/11/05 16:18

I know you can get age of the moon. Is it possible to calculate the phase and the Illumination (12%, 98% etc.) as well?

Tomás Alonso Albi, 2017/11/07 16:21

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):

phase = (1 - Math.cos(elong)) * 0.5;
ill = phase * 100;

Angular distance can be approximated from RA/DEC of both objects (you need to compute positions of both objects):

dlon = RA2 - RA1;
elong = Math.acos(Math.sin(DEC1) * Math.sin(DEC2) + Math.cos(DEC1) * Math.cos(DEC2) * Math.cos(dlon));

although you could use azimuth/elevation, its better to use RA/DEC because elevation could be affected by refraction.

Bruce Jin, 2017/11/19 01:38

Thanks!

Where can I get the values of RA1,RA2,DEC1,DEC2? Thanks!

Tomás Alonso Albi, 2017/11/20 17:38

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:

// Sun section
sunRA = out[6];
sunDEC = out[7];
 
// ...
// Moon section
 
moonRA = out[6];
moonDEC = out[7];

as you can see, RA and DEC are the elements 6 and 7 in the output array 'out', which I decided to ignore.

Bruce Jin, 2017/12/04 01:54

Thanks! did your code. Looks like the result matches internet!

Bruce Jin, 2018/01/14 17:31

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.

Tomás Alonso Albi, 2018/01/15 17:10

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.

Bruce Jin, 2018/01/15 22:48

Sorry. I was using local PC time which caused incorrect results. Thanks!

Ranjeet C, 2018/05/13 12:01

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?

Tomás Alonso Albi, 2018/05/17 11:03

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:

String phaseS[] = new String[] {"New Moon", "Crescent quarter", "Full Moon", "Descent quarter"};
double phaseD[] = new double[] {0, 0.25, 0.5, 0.75};
int phase = 1; // 0-3
double accuracy = 10 / (30 * SECONDS_PER_DAY); // 10s
double refPhase = phaseD[phase], jd_UT = smc.jd_UT;
while (true) {
	double age = normalizeRadians((smc.getMoon()[0] - smc.getSun()[0]) * DEG_TO_RAD) / 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 += 1-age;
	}
	smc.setUTDate(jd_UT);
}
System.out.println(phaseS[phase]+": "+getDateAsString(jd_UT));

where phase is the phase to compute (0-3).

Ranjeet C, 2018/05/20 11:24

Thanks to Tomás Alonso Albi, I got the idea now.

Ranjeet C, 2018/05/30 09:57

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?

Tomás Alonso Albi, 2018/05/30 17:59

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.

Ranjeet C, 2018/06/03 06:36

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.

Tomás Alonso Albi, 2018/06/04 18:01

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.

Marcel Stör, 2018/07/07 14:19

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/esp8266-weather-station.

Proper credit was given in three places: the README, the release notes, and of course the code

Tomás Alonso Albi, 2018/07/07 18:45

Glad you find the class useful, and thanks for the detailed credits

Tim Robinson, 2018/07/24 03:21

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

Tomás Alonso Albi, 2018/07/25 14:49

I have updated the program, with some reorganization of the code, to support those kind of calculations and to add also the lunar phases.

Marcel Stör, 2018/07/24 07:05

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.

Tomás Alonso Albi, 2018/07/25 14:53

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.

Tim Robinson, 2018/07/25 16:54

Thank you, much appreciated. Regards, Tim

Alfan, 2018/10/11 09:05

It's great. Thank you very much. It's very useful for creating Android project and save a lot of development time.

Sergey Shavarnaev, 2018/10/21 11:52

Hello. Im try to use moon.eclipticLongitude But get some unnormal value. I wait degree 0-360, 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

Tomás Alonso Albi, 2018/10/22 18:06

I didn't apply the reduction 0-2PI on that value. Just call the included function normalizeRadians on that value before converting radians to degrees.

ÇINAR, 2019/01/09 15:06

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

Tomás Alonso Albi, 2019/01/10 12:54

Sorry this is Java code, you will have to port it to javascript. You can partially save work using http://www.jsweet.org/.

Dev, 2019/01/18 08:34

I want to get longitude of sun and moon. How can i get ?

Tomás Alonso Albi, 2019/01/18 10:14

Can't understand what you ask for. Celestial longitude = right ascension, sun-centered 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 …

Ziva, 2019/04/01 14:37

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 = 6-3, 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!

Tomás Alonso Albi, 2019/04/01 15:08

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 re-create 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.

Ziva, 2019/04/02 15:42

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 = 6-3, 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.

Tomás Alonso Albi, 2019/04/02 16:33

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: 2019-04-01 02:50:47 UT Sun Transit: 2019-04-01 09:21:24 UT Sun Transit elevation: 39° 49' 32.3” Sun Set: 2019-04-01 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:

int year = 2019, month = 4, day = 1, h = 6-3, m = 28, s = 52; // in UT !!!
double obsLon = 40.64941 * Constant.DEG_TO_RAD, obsLat = 54.68653 * Constant.DEG_TO_RAD; // lon is negative to the west
Sascha Elble, 2019/06/17 00:07

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?

Tomás Alonso Albi, 2019/06/17 13:00

Hi Sascha,

Right after entering the page you mention the default coordinates and UT time I can see are:

int year = 2019, month = 6, day = 17, h = 10, m = 7, s = 45;
double obsLon = -76.8230111 * DEG_TO_RAD, obsLat = 42.6764593 * DEG_TO_RAD;

The results for the sidereal time on that page are:

APPARENT SIDEREAL TIME

Greenwich  03h 50m 26.29s = 3.8406367 h = 57.6095505 deg
Local      22h 43m 08.77s = 22.7191026 h = 340.7865394 deg

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 0-2PI 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:

smc.doCalc(smc.getSun(), false);
globalGMST = normalizeRadians(globalGMST);
globalLST = normalizeRadians(globalLST);
System.out.println("GMST: "+(float) (globalGMST * RAD_TO_DEG)+degSymbol);
System.out.println("AST:  "+(float) (globalLST * RAD_TO_DEG)+degSymbol);

My results are:

GMST: 57.31713°
AST:  340.4941°

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 TT-UT1 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

Jalil, 2019/07/24 19:59

with this input:

  lat = 36.0, lon = 3.0
  year = 2019, month = 7, day = 26
  hour = 0, minute = 0, second = 0

SunMoonCalculator output:

  moonrise : 2019-07-27 00:26:16
  moonset  : 2019-07-26 13:30:43

How to force it to output details for only the input date (2019-07-26), I mean:

  moonrise : no moon rise
  moonset  : 2019-07-26 13:30:43
Tomás Alonso Albi, 2019/07/25 10:38

At the end of the doCalc method you can read:

// 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 and 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 and Math.abs(set_time2) < Math.abs(set_time1)) set_time = set_time2;
rise = jd_UT + rise_time;
set = jd_UT + set_time;

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.

Tomás Alonso Albi, 2019/07/25 11:22

To be more specific, it should work if you replace the if sentences with

if (jdToday == riseToday2) rise_time = rise_time2;
 
...
 
if (jdToday == setToday2) set_time = set_time2;

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.

Jalil, 2019/08/06 14:40

How to calculate rise and set times taking in consideration the observers's elevation (altitude)? I think it is more accurate.

Tomás Alonso Albi, 2019/08/06 16:51

Being height the elevation in meters, compute the depresion of the horizon (<plus> is the plus symbol)

double eqRad = 6378.14, poRad = 6356.75; // IERS2003 Ellipsoid
double ratio =  (poRad / eqRad) * 0.5;
double rho = eqRad * (1.0 - ratio <plus> ratio * Math.cos(2.0 * obsLat));
double depresion = Math.acos(Math.sqrt(rho / (rho <plus> height / 1000.0)));

and in the twilight switch statement (doCalc method) include it, for instance for the standard 34' refraction at horizon:

tmp = -(34.0 / 60.0) * DEG_TO_RAD - pos[3] - depresion;
Jalil, 2019/08/06 22:22

I added this code just after the switch :

      if (obsElev != 0) {
          double eqRad = 6378.14, poRad = 6356.75; // IERS2003 Ellipsoid
          double ratio = (poRad / eqRad) * 0.5;
          double rho = eqRad * (1.0 - ratio + ratio * Math.cos(2.0 * obsLat));
          double depression = Math.acos(Math.sqrt(rho / (rho + obsElev / 1000.0)));
          tmp -= depression;
      }

will this work even when tmp is positive?

Tomás Alonso Albi, 2019/08/07 10:26

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.

Jalil, 2019/08/07 10:42

Exactly! I edited the switch as follow:

      switch (givenElevation) {
          case HORIZON_34arcmin_LIMB_UPPER:
              tmp = -(34.0 / 60.0) * DEG_TO_RAD - pos[3];
              break;
          case CUSTOM:
              tmp = customElevationDegrees * DEG_TO_RAD;
              break;
          case HORIZON_34arcmin_LIMB_LOWER:
              tmp = -(34.0 / 60.0) * DEG_TO_RAD + pos[3];
              break;
          case TWILIGHT_CIVIL:
              tmp = -6 * DEG_TO_RAD;
              break;
          case TWILIGHT_NAUTICAL:
              tmp = -12 * DEG_TO_RAD;
              break;
          case TWILIGHT_ASTRONOMICAL:
              tmp = -18 * DEG_TO_RAD;
              break;
      }
Jalil, 2019/08/07 23:53

need to allow the user to get geocentric details just by passing a parameter to:

  calcSunAndMoon(bool geocentric) 

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:

  Odeh, NACSA, Yallop, SAAO...
Tomás Alonso Albi, 2019/08/08 00:31

Sorry but I don't have time to help with that, if you need it then it is your homework.

Jalil, 2019/09/17 08:35

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?

Tomás Alonso Albi, 2019/09/17 10:57

Once again, it is your homework, just use a search-like iterative algorithm.

Applone, 2019/10/12 06:33

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

Tomás Alonso Albi, 2019/10/13 00:13

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 1-27 you should reduce it to that range with +27 or -27, so just do nakshtras=nakshtras-27 and you will get 1.

Dave, 2020/06/03 00:00

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?

Tomás Alonso Albi, 2020/06/03 11:51

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.

Dave, 2020/06/03 12:35

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.

Josh, 2020/06/06 23:08

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 open-source 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/.

Tomás Alonso Albi, 2020/06/07 11:21

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.

Marcel, 2020/06/07 11:30

2y I ported it to C for the Arduino/ESP8266. I gave credit in the code like this: https://github.com/ThingPulse/esp8266-weather-station/blob/master/src/SunMoonCalc.cpp#L23 My code is also MIT.

Tomás Alonso Albi, 2020/06/10 12:59

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.

Leandro, 2020/06/28 11:41

Hi!, a small contribution to this Project in Python language. Credits included.

GitHub respository: https://github.com/bokepasa/SunMoonCalculator

Thanks Tomás.

Tomás Alonso Albi, 2020/06/28 11:48

Not so small contribution,

¡Gracias Leandro!

Enter your comment
   __ __   __      __   ___    ____
  / // /  / /  __ / /  / _ |  / __/
 / _  /  / /__/ // /  / __ | / _/  
/_//_/  /____/\___/  /_/ |_|/_/
 
 
blog/sun_moon_position.txt · created: 2013/06/21 15:21 (Last modified 2020/08/07 18:36) by Tomás Alonso Albi
 
Recent changes RSS feed Creative Commons License Donate Powered by PHP Valid XHTML 1.0 Valid CSS Driven by DokuWiki