package lunisolar; import java.util.*; /* * Copyright 1996 by John D. Ramsdell * * This program is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License * as published by the Free Software Foundation; either version 2 * of the License, or (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */ /** Computes the phase of the moon and the sun. * @version March 1996 * @author John D. Ramsdell */ public final class Lunisolar { private Lunisolar() { } private final static double PI = Math.PI; private final static double HALF_PI = 0.5 * PI; private final static double TWO_PI = 2.0 * PI; private final static double RADIANS_PER_DEGREE = PI / 180.0; private final static long MINUTES_PER_DAY = 60 * 24; private final static long SECONDS_PER_DAY = 60 * MINUTES_PER_DAY; private final static long MILLISECONDS_PER_DAY = 1000 * SECONDS_PER_DAY; private final static long DAYS_PER_JULEAN_CENTURY = 36525; /* Time is most often represented as a double precision number */ /* in units of days. Angles are in radians. */ /* J2000 is the date January 1, 2000; 12:00:00 GMT */ /* This date is really called J2000.0. */ private static long compute_J2000() { GregorianCalendar gc = new GregorianCalendar(TimeZone.getTimeZone("GMT")); gc.set(Calendar.YEAR, 2000); gc.set(Calendar.MONTH, 0); gc.set(Calendar.DATE, 1); gc.set(Calendar.HOUR_OF_DAY, 12); gc.set(Calendar.MINUTE, 0); gc.set(Calendar.SECOND, 0); return gc.getTime().getTime(); } private static final long J2000 = compute_J2000(); private static final long fifty_years = DAYS_PER_JULEAN_CENTURY * MILLISECONDS_PER_DAY / 2; /* The formulas a valid between 1950 and 2050. */ private static final long J1950 = J2000 - fifty_years; private static final long J2050 = J2000 + fifty_years; /* Returns the current time, in units of days, after J2000.0. */ private static double days_after_J2000(long msecs) { return ((double)(msecs - J2000)) / MILLISECONDS_PER_DAY; } private static long the_time(double days) { return (long)Math.floor((double)MILLISECONDS_PER_DAY * days) + J2000; } /* Returns the angle between -PI < angle <= PI. */ private static double normalize_angle(double angle) { if (angle > PI) do angle -= TWO_PI; while (angle > PI); else while (angle <= -PI) angle += TWO_PI; return angle; } /*******************************************************************/ /* Astronomical almanac */ /* * All formulas are from: * The Astronomical Almanac for the Year 1984, * US Naval Observatory and Royal Greenwich Observatory, * US Government Printing Office, Washington DC, 1984. */ /* Angular position of the sun to a */ /* precision of 0.01 degrees. (Page C24). */ private final static double SUN0 = RADIANS_PER_DEGREE * 280.460; private final static double SUN1 = RADIANS_PER_DEGREE * 0.9856474; private final static double SUN2 = RADIANS_PER_DEGREE * 357.528; private final static double SUN3 = RADIANS_PER_DEGREE * 0.9856003; private final static double SUN4 = RADIANS_PER_DEGREE * 1.915; private final static double SUN5 = RADIANS_PER_DEGREE * 0.020; private static double sun_position(double days, double phase_offset) { double mean_longitude_of_sun = normalize_angle(SUN0 + SUN1 * days); double mean_anomaly = normalize_angle(SUN2 + SUN3 * days); double ecliptic_longitude = normalize_angle(mean_longitude_of_sun + SUN4 * Math.sin(mean_anomaly) + SUN5 * Math.sin(2.0 * mean_anomaly) - phase_offset); return ecliptic_longitude; } /* Angular velocity of the sun. Derivative of sun_position. */ private static double sun_velocity(double days) { double mean_anomaly = normalize_angle(SUN2 + SUN3 * days); return SUN1 + SUN4 * SUN3 * Math.cos(mean_anomaly) + SUN5 * 2.0 * SUN3 * Math.cos(2.0 * mean_anomaly); } /* Angular position of the moon to a */ /* precision of 0.3 degrees. (Page D46). */ private final static double RADIAN_CENTURY = RADIANS_PER_DEGREE / DAYS_PER_JULEAN_CENTURY; private final static double MOON0 = RADIANS_PER_DEGREE * 218.32; private final static double MOON1 = RADIAN_CENTURY * 481267.883; private final static double MOON2A = RADIANS_PER_DEGREE * 6.29; private final static double MOON2B = RADIANS_PER_DEGREE * 134.9; private final static double MOON2C = RADIAN_CENTURY * 477198.85; private final static double MOON3A = RADIANS_PER_DEGREE * -1.27; private final static double MOON3B = RADIANS_PER_DEGREE * 259.2; private final static double MOON3C = RADIAN_CENTURY * -413335.38; private final static double MOON4A = RADIANS_PER_DEGREE * 0.66; private final static double MOON4B = RADIANS_PER_DEGREE * 235.7; private final static double MOON4C = RADIAN_CENTURY * 890534.23; private final static double MOON5A = RADIANS_PER_DEGREE * 0.21; private final static double MOON5B = RADIANS_PER_DEGREE * 269.9; private final static double MOON5C = RADIAN_CENTURY * 954397.70; private final static double MOON6A = RADIANS_PER_DEGREE * -0.19; private final static double MOON6B = RADIANS_PER_DEGREE * 357.5; private final static double MOON6C = RADIAN_CENTURY * 035999.05; private final static double MOON7A = RADIANS_PER_DEGREE * -0.11; private final static double MOON7B = RADIANS_PER_DEGREE * 186.6; private final static double MOON7C = RADIAN_CENTURY * 966404.05; private static double moon_position(double days) { return normalize_angle(MOON0 + MOON1 * days + MOON2A * Math.sin(MOON2B + MOON2C * days) + MOON3A * Math.sin(MOON3B + MOON3C * days) + MOON4A * Math.sin(MOON4B + MOON4C * days) + MOON5A * Math.sin(MOON5B + MOON5C * days) + MOON6A * Math.sin(MOON6B + MOON6C * days) + MOON7A * Math.sin(MOON7B + MOON7C * days)); } /** * Computes the phase of the moon. * The formulas are valid only between 1950 and 2050. * @param msecs the time in milliseconds * @return the phase of the moon as a double in radians. * @exception Exception when date is not between 1950 and 2050 */ public static double moon_phase(long msecs) throws Exception { if (msecs < J1950 || msecs > J2050) throw new Exception("bad time value in moon_phrase"); double days = days_after_J2000(msecs); return -sun_position(days, moon_position(days)); } /****************************************************************/ private static double day_of_phase(double days, double phase) { for (int i = 0; i < 20; i++) // Newton's method days -= sun_position(days, phase) / sun_velocity(days); double last; do { last = days; // Newton's method checking result. days -= sun_position(days, phase) / sun_velocity(days); } while (Math.abs(days - last) >= (double)SECONDS_PER_DAY); return days; } /** * Finds a date giving the phase of the sun near the starting time. * @param msecs the start time in milliseconds * @param phase the phase of the sun as a double in radians. * @return the time in milliseconds * @exception Exception when date is not between 1950 and 2050 */ public static long time_of_phase(long msecs, double phase) throws Exception { if (msecs < J1950 || msecs > J2050) throw new Exception("bad time value in time_of_phase"); else return the_time(day_of_phase(days_after_J2000(msecs), phase)); } /****************************************************************/ /* Constructs an English sentence giving the current phase of the moon. */ private final static double PHASE_LIMIT = MOON1; /** * Constructs a phrase that describes the phase of the moon * @param msecs the time in milliseconds * @return the phase of the moon as a double in radians. * @exception Exception when date is not between 1950 and 2050 */ public static String moon_phrase(long msecs) throws Exception { double phase = moon_phase(msecs); /* Visable fraction. */ long percent = Math.round(50.0 * (1.0 - Math.cos(phase))); if (Math.abs (phase) < PHASE_LIMIT) return "new"; else if (Math.abs(normalize_angle (phase + PI)) < PHASE_LIMIT) return "full"; else if (Math.abs(phase - HALF_PI) < PHASE_LIMIT) return "first quarter (" + percent + "% of full)"; else if (Math.abs(phase + HALF_PI) < PHASE_LIMIT) return "last quarter (" + percent + "% of full)"; else if (phase > HALF_PI) return "waxing and gibbous (" + percent + "% of full)"; else if (phase > 0.0) return "a waxing crescent (" + percent + "% of full)"; else if (phase > -HALF_PI) return "a waning crescent (" + percent + "% of full)"; else return "waning and gibbous (" + percent + "% of full)"; } public static void main(String args[]) { try { System.out.println("Today's moon is " + moon_phrase(System.currentTimeMillis()) + "."); } catch (Exception e) { System.out.println(e); } } }