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.

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 but not returned). 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 force radiusAU to 0 to get the geocentric ones, although rise/set/transit values in that case will make no sense.

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. Input and output times are in UT.

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 can easily be extended to include the ephemerides of other planets. You only need some code to obtain the geocentric ecliptic longitude, latitude, and distance to a given planet referred to the equinox of date and use it instead of the getSun() and getMoon() methods. For this task I would recommend the classes VsopData and Vsop contained in the AstroLib package, although you should first ask for explicit permission from M. Huss. The method obtainAccurateRiseSetTransit would require also some update, replacing the boolean sun with an identifier for the planet.

Note: the program was extended in July 24, 2018, 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. In this work there are expansions to compute the positions of the planets that can be used instead of the method suggested above. The code has been changed to improve its quality and readability, so it is quite different from the original. The position of the Moon can probably be improved a little, I may revise this in the future.

SunMoonCalculator.java
/**
 * A very simple Sun/Moon calculator without using JPARSEC library.
 * @author T. Alonso Albi - OAN (Spain), email t.alonso@oan.es
 * @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;
 
	/** Radians to hours. */
	public static final double RAD_TO_HOUR = 180.0 / (15.0 * Math.PI);
 
	/** Radians to days. */
	public static final double RAD_TO_DAY = RAD_TO_HOUR / 24.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;
 
	/** The inverse of two times Pi. */
	public static final double TWO_PI_INVERSE = 1.0 / (2.0 * Math.PI);
 
	/** Four times Pi. */
	public static final double FOUR_PI = 4.0 * Math.PI;
 
	/** Pi divided by two. */
	public static final double PI_OVER_TWO = Math.PI / 2.0;
 
	/** Length of a sidereal day in days according to IERS Conventions. */
	public static final double SIDEREAL_DAY_LENGTH = 1.00273781191135448;
 
	/** 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;
 
	/** Our default epoch. <BR>
	 * 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.
		 */
		TWILIGHT_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.
		 */
		TWILIGHT_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.
		 */
		TWILIGHT_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;
		}
	}
 
	/** Input values. */
	private double jd_UT = 0, t = 0, obsLon = 0, obsLat = 0, TTminusUT = 0;
	private TWILIGHT twilight = TWILIGHT.HORIZON_34arcmin;
 
	/**
	 * 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) {
			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;
		}
 
		/** 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;
	}
 
	/** 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.328003326353 - 9.995238627481024 * x + 0.003067307630020489 * x2 - 7.76340698361363E-6 * x3 + 3.1331045394223196E-9 * x4 + 
						8.225530854405553E-12 * x2 * x3 - 7.486164715632051E-15 * x4 * x2 + 1.9362461549678834E-18 * x4 * x3 - 8.489224937827653E-23 * x4 * x4;
			} else {
				TTminusUT = -1027175.3477559977 + 2523.256625418965 * x - 1.885686849058459 * x2 + 5.869246227888417E-5 * x3 + 3.3379295816475025E-7 * x4 + 
						1.7758961671447929E-10 * x2 * x3 - 2.7889902806153024E-13 * x2 * x4 + 1.0224295822336825E-16 * x3 * x4 - 1.2528102370680435E-20 * x4 * x4;
			}
		}
		this.obsLon = obsLon;
		this.obsLat = obsLat;
		setUTDate(jd);
	}
 
	private 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;
	}
 
	/**
	 * 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;
	}
 
	private void setUTDate(double jd) {
		this.jd_UT = jd;
		this.t = (jd + TTminusUT / SECONDS_PER_DAY  - J2000) / JULIAN_DAYS_PER_CENTURY;
	}
 
	/** Calculates everything for the Sun and the Moon. */
	public void calcSunAndMoon() {
		double jd = this.jd_UT;
 
		// First the Sun
		sun = doCalc(getSun(), false);
 
		int niter = 3; // 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(), false).transitElevation;
		}
 
		// Now Moon
		setUTDate(jd);
		moon = doCalc(getMoon(), false);
		double ma = moonAge;
 
		niter = 5; // 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);
			getSun();
			moon.transitElevation = doCalc(getMoon(), false).transitElevation;
		}
		setUTDate(jd);
		moonAge = ma;
 
		// Compute illumination phase percentage for the Moon (do not use for other bodies!)
		double dlon = moon.rightAscension - sun.rightAscension;
		double elong = Math.acos(Math.sin(sun.declination) * Math.sin(moon.declination) + 
				Math.cos(sun.declination) * Math.cos(moon.declination) * Math.cos(dlon));
		moon.illuminationPhase = 100 * (1.0 - Math.cos(elong)) * 0.5;
	}
 
	// Formulae here is a simplifcation of the expansion from 
	// "Planetary Programs and Tables" by Pierre Bretagnon and
	// Jean-Louis Simon, Willman-Bell, 1986. This source also 
	// have expansions for ephemerides of planets
	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 }
	};
 
	private double[] getSun() {
		double L = 0.0, R = 0.0;
		double t2 = t * 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 lon = normalizeRadians(4.9353929 + normalizeRadians(62833.196168 * t2) + L / 10000000.0) * RAD_TO_DEG;
		double sdistance = 1.0001026 + R / 10000000.0;
 
		// Now, let calculate 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;
		double d = -.00569 - .0047785 * Math.sin(M1) - .0003667 * Math.sin(M2);
 
		double slongitude = lon + d; // apparent longitude (error<0.001 deg)
		double slatitude = 0; // Sun's ecliptic latitude is always negligible
 
		return new double[] {slongitude, slatitude, sdistance, Math.atan(696000 / (AU * sdistance))};
	}
 
	private 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;
 
		// Let's add nutation here also
		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;
		double d = - .0047785 * Math.sin(M1) - .0003667 * Math.sin(M2);
		longitude += d;
 
		// Get accurate Moon age
		double Psin = 29.530588853;
		moonAge = normalizeRadians(longitude * DEG_TO_RAD - sun.eclipticLongitude) * 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;
 
		return new double[] {longitude, latitude, distance * EARTH_RADIUS / AU, Math.atan(1737.4 / (distance * EARTH_RADIUS))};
	}
 
	private Ephemeris doCalc(double[] pos, boolean geocentric) {
		// Ecliptic to equatorial coordinates
		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;
		double angle = (23.4392911111111 + tmp) * DEG_TO_RAD; // obliquity
 
		// Add nutation in obliquity
		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;
		double d = .002558 * Math.cos(M1) - .00015339 * Math.cos(M2);
		angle += d * DEG_TO_RAD;
 
		pos[0] *= DEG_TO_RAD;
		pos[1] *= DEG_TO_RAD;
		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]);
		tmp = y * Math.cos(angle) - z * Math.sin(angle);
		z = y * Math.sin(angle) + z * Math.cos(angle);
		y = tmp;
 
		if (geocentric) return new Ephemeris(0, 0, -1, -1, -1, -1, normalizeRadians(Math.atan2(y, x)), 
				Math.atan2(z / Math.sqrt(x * x + y * y), 1.0), 	Math.sqrt(x * x + y * y + z * z), pos[0], pos[1]);
 
		// 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;
		double lst = gmst + obsLon;
 
		// Obtain topocentric rectangular coordinates
		double radiusAU = EARTH_RADIUS / AU;
		double correction[] = new double[] {
				radiusAU * Math.cos(obsLat) * Math.cos(lst),
				radiusAU * Math.cos(obsLat) * Math.sin(lst),
				radiusAU * Math.sin(obsLat)};
		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.sqrt(xtopo * xtopo + ytopo * 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 sinlat = Math.sin(obsLat); 
		double coslat = Math.cos(obsLat); 
		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
 
		// Get apparent elevation
		if (alt > -3 * DEG_TO_RAD) {
			double r = 0.016667 * DEG_TO_RAD * Math.abs(Math.tan(PI_OVER_TWO - (alt * RAD_TO_DEG +  7.31 / (alt * RAD_TO_DEG + 4.4)) * DEG_TO_RAD));
			double refr = r * ( 0.28 * 1010 / (10 + 273.0)); // Assuming pressure of 1010 mb and T = 10 C
			alt = Math.min(alt + refr, PI_OVER_TWO); // This is not accurate, but acceptable
		}
 
		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 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;
		}
 
		// Compute cosine of hour angle
		tmp = (Math.sin(tmp) - Math.sin(obsLat) * Math.sin(dec)) / (Math.cos(obsLat) * Math.cos(dec));
		double celestialHoursToEarthTime = RAD_TO_DAY / SIDEREAL_DAY_LENGTH;
 
		// 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(Math.sin(dec) * Math.sin(obsLat) + Math.cos(dec) * Math.cos(obsLat));
		if (transit_alt > -3 * DEG_TO_RAD) {
			double r = 0.016667 * DEG_TO_RAD * Math.abs(Math.tan(PI_OVER_TWO - (transit_alt * RAD_TO_DEG +  7.31 / (transit_alt * RAD_TO_DEG + 4.4)) * DEG_TO_RAD));
			double refr = r * ( 0.28 * 1010 / (10 + 273.0)); // Assuming pressure of 1010 mb and T = 10 C
			transit_alt = Math.min(transit_alt + refr, PI_OVER_TWO); // This is not accurate, but acceptable
		}
 
		// 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]);
		return out;
	}
 
	/**
	 * 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.
	 */
	public 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[] = SunMoonCalculator.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.
	 */
	public static double normalizeRadians(double r)
	{
		if (r < 0 && r >= -TWO_PI) return r + TWO_PI;
		if (r >= TWO_PI && r < FOUR_PI) return r - TWO_PI;
		if (r >= 0 && r < TWO_PI) return r;
 
		r -= TWO_PI * Math.floor(r * TWO_PI_INVERSE);
		if (r < 0.) r += TWO_PI;
 
		return r;
	}
 
	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(), false);
			} else {
				getSun();
				out = doCalc(getMoon(), false);
			}
 
			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) return -1; // did not converge => without rise/set/transit in this date
		return riseSetJD;
	}
 
	/**
	 * 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;
		while (true) {
			double sunLon = getSun()[0]; // Compute Sun before always !
			double age = normalizeRadians((getMoon()[0] - sunLon) * 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;
			}
			setUTDate(jd_UT);
		}
		double out = jd_UT;
		setUTDate(oldJD);
		return out;
	}
 
	/**
	 * Returns the dates of the Spring and Autumn equinoxes.
	 * Note the algorithm is not speed-optimized.
	 * @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];
 
		int year = getDate(jd_UT)[0];
		double precTime = 30 / 86400.0; // 30s
 
		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;
			while (true) {
				double dec = doCalc(getSun(), true).declination;
				if (Math.abs(dec) < min || min == -1) {
					min = Math.abs(dec);
					minT = jd;
				}
				if (Math.abs(dec) > min && min >= 0) {
					out[i] = minT;
					break;
				}
				jd += precTime;
				setUTDate(jd);
			}
		}
		setUTDate(jdOld);
		doCalc(getSun(), false);
		return out;
	}
 
	/**
	 * Returns the dates of the Summer and Winter solstices.
	 * Note the algorithm is not speed-optimized.
	 * @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];
 
		int year = getDate(jd_UT)[0];
		double precTime = 30 / 86400.0; // 30s
 
		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;
			while (true) {
				double dec = doCalc(getSun(), true).declination;
				if (Math.abs(dec) > max || max == -1) {
					max = Math.abs(dec);
					maxT = jd;
				}
				if (Math.abs(dec) < max && max >= 0) {
					out[i] = maxT;
					break;
				}
				jd += precTime;
				setUTDate(jd);
			}
		}
		setUTDate(jdOld);
		doCalc(getSun(), false);
		return out;
	}
 
	/**
	 * Main test program.
	 * @param args Not used
	 */
	public static void main(String[] args) {
		System.out.println("SunMoonCalculator test run");
 
		try {
			int year = 2018, month = 8, day = 1, h = 12, m = 0, s = 0; // in UT !!!
			double obsLon = -73.56 * DEG_TO_RAD, obsLat = 40.40 * DEG_TO_RAD;
			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(" Rise:    "+SunMoonCalculator.getDateAsString(smc.sun.rise));
			System.out.println(" Set:     "+SunMoonCalculator.getDateAsString(smc.sun.set));
			System.out.println(" Transit: "+SunMoonCalculator.getDateAsString(smc.sun.transit)+" (max. elev. "+(float) (smc.sun.transitElevation * RAD_TO_DEG)+degSymbol+")");
			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(" 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)+" (max. elev. "+(float) (smc.moon.transitElevation * RAD_TO_DEG)+degSymbol+")");
 
			smc.setTwilight(TWILIGHT.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 in over 1800 - 2200:
			// - Sun: 0.001 deg or 3 arcsec. <1s in rise/set/transit times. 1 min in Equinoxes/Solstices
			// - Mon: 0.02 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 from 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. The following algorithm uses some methods from the previous program to compute the necessary angles. First two values 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 can be directly computed from the moon age, knowing that the cycle is 29.5306 days and the lunar disk will be fully illuminated at cycle/2.

Of course 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. Some code lines could be removed by integrating the method in the previous code, as well as the initial calls to computations.

public double[] getMoonDiskOrientationAngles()
{
	sun = doCalc(getSun(), false);
	double moonPos[] = getMoon();
	moon = doCalc(moonPos, false);
	double moonLon = moonPos[0], moonLat = moonPos[1], 
        	moonRA = moon.rightAscension, moonDEC = moon.declination;
	double sunRA = sun.rightAscension, sunDEC = sun.declination;
 
	// 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;
	// Obliquity of ecliptic (approx, better formulae up)
	double eps = 23.43929 * DEG_TO_RAD;
 
	// Obtain optical librations lp and bp
	double W = moonLon - omega;
	double sinA = Math.sin(W) * Math.cos(moonLat) * Math.cos(I) - Math.sin(moonLat) * Math.sin(I);
	double cosA = Math.cos(W) * Math.cos(moonLat);
	double A = Math.atan2(sinA, cosA);
	double lp = normalizeRadians(A - F);
	double sinbp = - Math.sin(W) * Math.cos(moonLat) * Math.sin(I) - Math.sin(moonLat) * Math.cos(I);
	double bp = Math.asin(sinbp);
 
	// Obtain position angle of axis p
	double x = Math.sin(I) * Math.sin(omega);
	double y = Math.sin(I) * Math.cos(omega) * Math.cos(eps) - Math.cos(I) * Math.sin(eps);
	double w = Math.atan2(x, y);
	double sinp = Math.sqrt(x*x + y*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) * Math
			.sin(moonDEC) * Math.cos(moonRA - sunRA) - Math
			.sin(sunDEC) * Math.cos(moonDEC)));
 
	// Paralactic angle par (first 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;
	double lst = gmst + obsLon;
 
	y = Math.sin(lst - moonRA);
	x = Math.tan(obsLat) * Math.cos(moonDEC) - Math.sin(moonDEC) * 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.

At the end of the discussion there is also a link to a port to C++ for an Arduino-based project.

NOTE: this program may still use the old version of the Sun/Moon calculator.

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

Enter your comment
 
 
blog/sun_moon_position.txt · created: 2013/06/21 15:21 (Last modified 2018/09/06 17:40) 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