* * Permission to use, copy, modify, and redistribute this software and its * documentation for personal, non-commercial use is hereby granted provided that * this copyright notice and appropriate documentation appears in all copies. This * software may not be distributed for fee or as part of commercial, "shareware," * and/or not-for-profit endevors including, but not limited to, CD-ROM collections, * online databases, and subscription services without specific license.
*
* @author Martin Minow
* @version 1.0
* Set tabs every 4 characters.
*/
package org.minow.applets.sunsphere;
import java.util.*;
/**
* SunSphereRiseSet computes rise, set, and twilight of the sun and moon.
*/
public class SunSphereRiseSet
{
/**
* ABOVE_HORIZON and BELOW_HORIZON are returned for sun and moon
* calculations where the astronomical object does not cross the
* horizon.
*/
public static final double ABOVE_HORIZON = Double.POSITIVE_INFINITY;
public static final double BELOW_HORIZON = Double.NEGATIVE_INFINITY;
/*
* These selectors are passed to sunRiseSet
*/
public static final double RISING = +1.0; /* RISING and SETTING */
public static final double SETTING = -1.0; /* configure ephemeris */
public static final double DegRad = (Math.PI / 180.0);
public static final double javaEpochMJD = 40587.0;
public static final double msecPerDay = 86400000.0;
/*
* Return the MJD for the given date. Note that this is always in
* the local timezone.
* @param javaDate The desired date.
* @return the Modified Julian Date for this date.
* Test: 1957.10.04:19:28:00 -> 2 436 116.31 (Launch of Sputnik 1)
* 1985.02.17:06:00:00 -> 2 446 113.75 (Duffett-Smith sample)
* 1980.07.23:12.00.00 -> 2 444 444.00 (Practical Astronomy)
*/
public static double getMJD(
Date javaDate
)
{
long instant = javaDate.getTime() + Utility.getTimeZoneOffset(javaDate);
return (((double) instant) / ((double) msecPerDay) + javaEpochMJD);
}
/**
* Compute the time of sunrise and sunset for this date. This
* is a low-precision iterative algorithm following the "systematic" approach
* outlined in the Explanatory Supplement to the Astronomical Almanac.
* It is loses accuracy above 60¡ latitude, but is suitable for image
* filtering.
*
* The times are returned in the specified timezone.
* @parm MJD The MJD at midnight at the observer's location.
* @param zone The timezone offset in hours, East postive.
* @param latitude The observer's latitude, North positive
* @param longitude The observer's longitude, East positive
* @param horizon The adopted true altitude of the horizon
* in degrees. Use one of the values defined
* in Astro.java, as noted below.
* @param riseSet +1.0 for rise, -1.0 for set.
* @return The time that the sun crosses the horizon.
*
ABOVE_HORIZON Returned if the sun never rises
*
BELOW_HORIZON Returned if the sun never sets
*
* The following values maybe used for the horizon parameter:
*
-0¡50' SUNRISE
*
-6¡00' CIVIL_TWILIGHT
*
-12¡00' NAUTICAL_TWILIGHT
*
-18¡00' ASTRONOMICAL_TWILIGHT
*
* Here are some test values. The "correct" values were taken from
* the United States Naval
* Observatory Web page.
* 1997.01.01, latitude 0.0, longitude 0.0, zone 0.0
* SunRiseSet Civil Twil. Naut Twil.
* 06:00 18:07 05:37 18:30 05:11 18:56
*
* 1997.02.21, latitude 37.8N, longitude 122.4W, zone -8.0
* SunRiseSet Civil Twil.
* 06:52 17:56 06:25 18:22
*
*/
public static double sunRiseSet(
double MJD, /* MJD at midnight */
double zone, /* Timezone offset in hours */
double latitude, /* latitude, North positive */
double longitude, /* longitude, East positive */
double horizon, /* horizon (for twilight) */
double riseSet /* +1.0 == rising, -1.0 = setting */
)
{
double sinLatitude = Math.sin(latitude * DegRad);
double cosLatitude = Math.cos(latitude * DegRad);
double sinHorizon = Math.sin(horizon * DegRad);
double UT = 0.0;
int aboveCount = 0;
int belowCount = 0;
/*
* Ensure that there is a rise/set at this latitude. This will
* be the case if startBelow is below the horizon and endBelow
* above the horizon or vice-versa
*/
double UT0 = 12.0;
for (int i = 0; i < 20 && aboveCount < 3 && belowCount < 3; i++) {
/*
* Compute the solar ephemeris
* T Number of centuries from J2000.
* (Original used JD, we use MJD).
*/
double T = (MJD + (UT / 24.0) - 51544.5) / 36525.0;
/*
* Solar arguments
* L Mean longitude corrected for aberration (degrees)
* G Mean anomaly (radians)
* Lambda Ecliptic longitude (radians)
* Epsilon The obliquity of the ecliptic (radians)
*/
double L = MOD(280.460 + 36000.770 * T, 360.0);
double G = MOD(357.528 + 35999.050 * T, 360.0);
double gRad = G * DegRad;
double sinG = Math.sin(gRad);
double sin2G = Math.sin(gRad * 2.0);
double gFactor = 1.915 * sinG + 0.020 * sin2G;
double Lambda = L + gFactor;
double lambdaRad = Lambda * DegRad;
double Epsilon = 23.4393 - 0.01300 * T;
double epsilonRad = Epsilon * DegRad;
/*
* Compute the ephemeris
* E Equation of time
* GHA Greenwich Hour Angle (degrees)
* sinDeclin sin(Declination)
* cosDeclin cos(Declination)
*/
double E = -gFactor
+ 2.466 * Math.sin(lambdaRad * 2.0)
- 0.053 * Math.sin(lambdaRad * 4.0);
double GHA = ((15.0 * UT) - 180.0 + E);
double sinDeclin = Math.sin(epsilonRad) * Math.sin(lambdaRad);
double rDeclination = Math.asin(sinDeclin);
double cosDeclin = Math.cos(rDeclination);
/*
* The hour angle at time UT0 is HA, where cosHA is
*/
double HA;
double cosHA;
try {
cosHA = (sinHorizon - (sinLatitude * sinDeclin))
/ (cosLatitude * cosDeclin);
}
catch (ArithmeticException e) { /* At +/- 90 degrees latitude */
cosHA = (sinDeclin >= 0) ? 2.0 : -2.0; /* Hack */
}
if (cosHA > +1.0) {
++belowCount;
HA = 0;
}
else if (cosHA < -1.0) {
++aboveCount;
HA = 180.0;
}
else {
aboveCount = 0;
belowCount = 0;
HA = Math.acos(cosHA) / DegRad;
}
UT = UT0 - (GHA + longitude + (riseSet * HA)) / 15.0;
UT = MOD(UT, 24.0);
if (Math.abs(UT - UT0) < 0.0008)
break;
UT0 = UT;
}
if (belowCount > 0)
UT = BELOW_HORIZON;
else if (aboveCount > 0)
UT = ABOVE_HORIZON;
else {
UT = MOD(UT + zone, 24.0);
}
return (UT);
}
/*
* Modulus function that always returns a positive value. For example,
* AstroMath.mod(-3, 24) == 21
*/
public static double MOD(
double numerator,
double denomenator
)
{
double result = Math.IEEEremainder(numerator, denomenator);
if (result < 0.0)
result += denomenator;
return (result);
}
}