// This is a library of calculations for calculating astronomical positions.
// Most algorithms are based on Jean Meeus book _Astronomical Algorithms_
// Sign conventions
// It is *highly* recommended that you use these constants instead of
// assuming a sign convention. These are the conventions used in Meeus and
// many other astronomical calculations, which often differ than sign
// conventions used in geoinformatics.
North := +1
South := -1
West := +1
East := -1
// Turn a latitude measurement into the appropriate letter or word. (e.g.
// "N" or "S"
// Pass in the latitude and an optional flag to produce the letter
// or word. Default is letter only.
latitudeName[lat, short=true] :=
{
if lat >= 0 degrees
short ? "N" : "North"
else
short ? "S" : "South"
}
// Turn a longitude measurement into the appropriate letter or word. (e.g.
// "N" or "S"
// Pass in the longitude and an optional flag to produce the letter
// or word. Default is letter only.
longitudeName[long, short=true] :=
{
if long >= 0 degrees
short ? "W" : "West"
else
short ? "E" : "East"
}
// Returns the value of omega from Chapter 25 of Jean Meeus.
meeusOmega[date] :=
{
T = meeusT[date]
return (125.04 - 1934.136 T) degree
}
// Returns the apparent longitude of the sun (lambda) (low accuracy)
// from Meeus, Chapter 25
sunApparentLongitude[date] :=
{
T = meeusT[date]
omega = meeusOmega[date]
// println["omega: " + (omega mod (360 degree) -> "degrees")];
return sunTrueLongitude[date] - 0.00569 degree - 0.00478 degree sin[omega]
}
// Returns geometric mean longitude L0 of the sun, referred to the mean equinox
// of the date. Equation 25.2 in Jean Meeus.
sunGeometricMeanLongitude[date] :=
{
T = meeusT[date]
return (280.46646 + 36000.76983 T + 0.0003032 T^2) degree
}
// Returns true longitude (sunSymbol) of the sun for a specified date.
// From Chapter 25 in Jean Meeus.
sunTrueLongitude[date] :=
{
return sunGeometricMeanLongitude[date] + sunCenter[date]
}
// Returns the mean anomaly M of the sun for a given date.
// Equation 25.3 in Jean Meeus.
sunMeanAnomaly[date] :=
{
T = meeusT[date]
return ((357.52911 + 35999.05029 T - 0.0001537 T^2) degree) mod (360 degree)
}
// Returns the true anomaly nu of the sun for a given date.
// From Chapter 25 of Meeus.
sunTrueAnomaly[date] :=
{
T = meeusT[date]
return sunMeanAnomaly[date] + sunCenter[date]
}
// Equation for center of the Sun C from Chapter 25 of Meeus.
sunCenter[date] :=
{
T = meeusT[date]
M = sunMeanAnomaly[date]
return ((1.914602 - 0.004817 T - 0.000014 T^2) degree) sin[M] +
((0.019993 - 0.000101 T) degree) sin[2 M] +
(0.000289 degree) sin[3 M]
}
// Returns the eccentricity of the earth's orbit for a given date.
// Equation 25.4 in Meeus.
earthEccentricity[date] :=
{
T = meeusT[date]
return 0.016708634 - 0.000042037 T - 0.0000001267 T^2
}
// Distance between center of Sun and Earth
sunDistance[date] :=
{
eccentricity = earthEccentricity[date]
trueAnomaly = sunTrueAnomaly[date]
return 1.000001018 au ( 1 - eccentricity^2 ) / (1 + eccentricity cos[trueAnomaly])
}
// Calculates the apparent sun radius angle, corrected for distance, but not
// corrected for refraction. It gives a true arcsin geometric interpretation.
sunRadiusAngle[date] := arcsin[sunradius / (sunradius + sunDistance[date])]
// This returns the number of Julian centuries since JDE 2451545 as a
// dimensionless number.
// Equation 25.1 in Jean Meeus
meeusT[date] := (JDE[date] - 2451545 days) / juliancentury
// Mean obliquity of the ecliptic epsilon0 from Jean Meeus, eq. 22.3. This
// is a high accuracy version of the function that is valid for a period
// of 10000 years before or after J2000.0
// It is NOT corrected for nutation.
// In Meeus, this symbol is referred to as epsilon_0
meanObliquityOfEcliptic[date] :=
{
T = meeusT[date]
U = T/100
return 23 degree + 26 arcmin + 21.448 arcsec +
(-4680.93 U - 1.55 U^2 + 1999.25 U^3 - 51.38 U^4 - 249.67 U^5 +
-39.05 U^6 + 7.12 U^7 + 27.87 U^8 + 5.79 U^9 + 2.45 U^10) arcsec
}
// True obliquity of the ecliptic corrected for nutation. This is called
// epsilon in Meeus.
trueObliquityOfEcliptic[date] :=
{
epsilon0 = meanObliquityOfEcliptic[date]
[deltapsi, deltaepsilon] = highAccuracyNutation[date]
return epsilon0 + deltaepsilon
}
// Mean obliquity of the ecliptic epsilon0 from Jean Meeus, eq. 22.3. This
// is a high accuracy version of the function that is only valid for a period
// of 10000 years before or after J2000.0
// In Meeus, this symbol is referred to as epsilon_0
narrowValidityMeanObliquityOfEcliptic[date] :=
{
T = meeusT[date]
return 23 degree + 26 arcmin + 21.448 arcsec +
(-46.8150 T - 0.00059 T^2 + 0.001813 T^3) arcsec
}
// Returns the apparent right ascension and declination of the sun for a
// given date. It implements the corrected versions of equation 25.6
// and 25.7 in Meeus.
// Returns [ra, decl]
sunApparentRADecl[date] :=
{
epsilon = trueObliquityOfEcliptic[date]
epsilonprime = epsilon + 0.00256 degree cos[meeusOmega[date]]
sunLong = sunApparentLongitude[date]
ssl = sin[sunLong]
ra = arctan[cos[epsilonprime] ssl, cos[sunLong]]
decl = arcsin[sin[epsilonprime] ssl]
return [ra, decl]
}
// Calculate mean Sidereal time theta_0 at Greenwich for any instant as an
// angle.
// Equation 12.4 in Meeus.
meanGreenwichSiderealAngle[date] :=
{
T = meeusT[date]
return ((280.46061837 + 360.98564736629 (JD[date]/day - 2451545) +
0.000387933 T^2 - T^3/38710000) degree) mod circle
}
// Calculate apparent Sidereal time as an angle at Greenwich for any instant.
// From Chapter 12 of Meeus.
apparentGreenwichSiderealAngle[date] :=
{
deltapsi = highAccuracyNutation[date]@0
trueObliquity = trueObliquityOfEcliptic[date]
correction = deltapsi cos[trueObliquity]
return (meanGreenwichSiderealAngle[date] + correction) mod circle
}
// Calculate apparent local sidereal time as an angle.
apparentLocalSiderealAngle[date, longitude] :=
{
angle = apparentGreenwichSiderealAngle[date]
return angle - longitude
}
// Calculate the hour angle H of a body with right ascension ra.
// See chapter 13 in Meeus.
hourAngle[date, longitude, ra] :=
{
return apparentLocalSiderealAngle[date, longitude] - ra
}
// Returns the parallactic angle q for a given body from a point [lat, long]
// on the earth of a body at right ascension ra
// and declination decl.
// See equation 14.1 in Meeus.
parallacticAngle[date, lat, long, ra, decl] :=
{
H = hourAngle[date, long, ra]
return arctan[sin[H], tan[lat] cos[decl] - sin[decl] cos[H]]
}
// Calculate the approximate date of the equinoxes and solstices
// (to medium accuracy.) The value should be correct within a minute;
// see Table 27.D in Meeus.
// This is from chapter 27 of Meeus.
// This function is called by the equinox and solstice methods below.
equinoxDate[JDE0] :=
{
T = (JDE0 - 2451545.0) / 36525
W = (35999.373 T - 2.47) degrees
deltaLambda = 1 + 0.0334 cos[W] + 0.0007 cos[2 W]
// These values are from Table 27.C
S = 485 cos[(324.96 + 1934.136 T) degrees] +
203 cos[(337.23 + 32964.467 T) degrees] +
199 cos[(342.08 + 20.186 T) degrees] +
182 cos[( 27.85 + 445267.112 T) degrees] +
156 cos[( 73.14 + 45036.886 T) degrees] +
136 cos[(171.52 + 22518.443 T) degrees] +
77 cos[(222.54 + 65928.934 T) degrees] +
74 cos[(296.72 + 3034.906 T) degrees] +
70 cos[(243.58 + 9037.513 T) degrees] +
58 cos[(119.81 + 33718.147 T) degrees] +
52 cos[(297.17 + 150.678 T) degrees] +
50 cos[( 21.02 + 2281.226 T) degrees] +
45 cos[(247.54 + 29929.562 T) degrees] +
44 cos[(325.15 + 31555.956 T) degrees] +
29 cos[( 60.93 + 4443.417 T) degrees] +
18 cos[(155.12 + 67555.328 T) degrees] +
17 cos[(288.79 + 4562.452 T) degrees] +
16 cos[(198.04 + 62894.029 T) degrees] +
14 cos[(199.76 + 31436.921 T) degrees] +
12 cos[( 95.39 + 14577.848 T) degrees] +
12 cos[(287.11 + 31931.756 T) degrees] +
12 cos[(320.81 + 34777.259 T) degrees] +
9 cos[(227.73 + 1222.114 T) degrees] +
8 cos[( 15.45 + 16859.074 T) degrees]
ss = JDE0 + 0.00001 S/deltaLambda
return JDE[ss]
}
// Spring Equinox for the years +1000 to +3000, Table 27.B in Meeus.
// The year passed in must be an integer!
springEquinox[year] :=
{
Y = (year - 2000)/1000
JDE0 = 2451623.80984 + 365242.37404 Y + 0.05169 Y^2 - 0.00411 Y^3 -
0.00057 Y^4
return equinoxDate[JDE0]
}
// Summer Solstice for the years +1000 to +3000, Table 27.B in Meeus.
// The year passed in must be an integer!
summerSolstice[year] :=
{
Y = (year - 2000)/1000
JDE0 = 2451716.56767 + 365241.62603 Y + 0.00325 Y^2 + 0.00888 Y^3 -
0.00030 Y^4
return equinoxDate[JDE0]
}
// Autumn Equinox for the years +1000 to +3000, Table 27.B in Meeus.
// The year passed in must be an integer!
autumnEquinox[year] :=
{
Y = (year - 2000)/1000
JDE0 = 2451810.21715 + 365242.01767 Y - 0.11575 Y^2 + 0.00337 Y^3 +
0.00078 Y^4
return equinoxDate[JDE0]
}
// Winter Solstice for the years +1000 to +3000, Table 27.B in Meeus.
// The year passed in must be an integer!
winterSolstice[year] :=
{
Y = (year - 2000)/1000
JDE0 = 2451900.05952 + 365242.74049 Y - 0.06223 Y^2 - 0.00823 Y^3 +
0.00032 Y^4
return equinoxDate[JDE0]
}
// Calculate Nutation in longitude (delta psi) and
// obliquity (delta epsilon) to low accuracy.
// Returns [deltapsi, deltaepsilon]
// From chapter 22 in Meeus.
lowAccuracyNutation[date] :=
{
T = meeusT[date]
omega = (125.04452 - 1934.136261 T) degree
// Mean longitude of the sun
L = (280.4665 + 36000.7698 T) degree
// Mean longitude of the moon
Lmoon = (218.3165 + 481267.8813 T) degree
deltapsi = (-17.20 sin[omega] - 1.32 sin[2 L] - 0.23 sin[2 Lmoon] +
0.21 sin[2 omega] ) arcsec
deltaepsilon = (9.20 cos[omega] + 0.57 cos[2 L] + 0.10 cos[2 Lmoon] -
0.09 cos[2 omega] ) arcsec
return [deltapsi, deltaepsilon]
}
// Calculate Nutation in longitude (delta psi) and
// obliquity (delta epsilon) to high
// Returns [deltapsi, deltaepsilon]
// From chapter 22 in Meeus.
highAccuracyNutation[date] :=
{
T = meeusT[date]
// println["T: $T"]
// Mean elongation of the moon from the sun:
D = (297.85036 + 445267.111480 T - 0.0019142 T^2 + T^3/189474) degree
// println["D: " + (D->"degree")]
// Mean anomaly of the Sun:
M = (357.52772 + 35999.050340 T - 0.0001603 T^2 - T^3/300000) degree
// println["M: " + (M->"degree")]
// Mean anomaly of the Moon:
MP = (134.96298 + 477198.867398 T + 0.0086972 T^2 + T^3/56250) degree
// println["MP: " + (MP->"degree")]
// Moon's argument of Latitude:
F = (93.27191 + 483202.017538 T - 0.0036825 T^2 + T^3/327270) degree
// println["F: " + (F->"degree")]
// Longitude of the ascending node of the moon's mean orbit on the ecliptic
// measured from the mean equinox of the date.
Omega = (125.04452 - 1934.136261 T + 0.0020708 T^2 + T^3/450000) degree
// println["Omega: " + (Omega->"degree")]
// These lines generated by iau1980.frink and pasted in here.
deltapsi = 0.0001 arcsec ( + (-171996 + -174.2 T) sin[ + Omega] + (-13187 + -1.6 T) sin[ -2 D + 2 F + 2 Omega] + (-2274 + -0.2 T) sin[ + 2 F + 2 Omega] + (2062 + 0.2 T) sin[ + 2 Omega] + (1426 + -3.4 T) sin[ + M] + (712 + 0.1 T) sin[ + MP] + (-517 + 1.2 T) sin[ -2 D + M + 2 F + 2 Omega] + (-386 + -0.4 T) sin[ + 2 F + Omega] + (-301) sin[ + MP + 2 F + 2 Omega] + (217 + -0.5 T) sin[ -2 D -1 M + 2 F + 2 Omega] + (-158) sin[ -2 D + MP] + (129 + 0.1 T) sin[ -2 D + 2 F + Omega] + (123) sin[ -1 MP + 2 F + 2 Omega] + (63) sin[ + 2 D] + (63 + 0.1 T) sin[ + MP + Omega] + (-59) sin[ + 2 D -1 MP + 2 F + 2 Omega] + (-58 + -0.1 T) sin[ -1 MP + Omega] + (-51) sin[ + MP + 2 F + Omega] + (48) sin[ -2 D + 2 MP] + (46) sin[ -2 MP + 2 F + Omega] + (-38) sin[ + 2 D + 2 F + 2 Omega] + (-31) sin[ + 2 MP + 2 F + 2 Omega] + (29) sin[ + 2 MP] + (29) sin[ -2 D + MP + 2 F + 2 Omega] + (26) sin[ + 2 F] + (-22) sin[ -2 D + 2 F] + (21) sin[ -1 MP + 2 F + Omega] + (17 + -0.1 T) sin[ + 2 M] + (16) sin[ + 2 D -1 MP + Omega] + (-16 + 0.1 T) sin[ -2 D + 2 M + 2 F + 2 Omega] + (-15) sin[ + M + Omega] + (-13) sin[ -2 D + MP + Omega] + (-12) sin[ -1 M + Omega] + (11) sin[ + 2 MP -2 F] + (-10) sin[ + 2 D -1 MP + 2 F + Omega] + (-8) sin[ + 2 D + MP + 2 F + 2 Omega] + (7) sin[ + M + 2 F + 2 Omega] + (-7) sin[ -2 D + M + MP] + (-7) sin[ -1 M + 2 F + 2 Omega] + (-8) sin[ + 2 D + 2 F + Omega] + (6) sin[ + 2 D + MP] + (6) sin[ -2 D + 2 MP + 2 F + 2 Omega] + (6) sin[ -2 D + MP + 2 F + Omega] + (-6) sin[ + 2 D -2 MP + Omega] + (-6) sin[ + 2 D + Omega] + (5) sin[ -1 M + MP] + (-5) sin[ -2 D -1 M + 2 F + Omega] + (-5) sin[ -2 D + Omega] + (-5) sin[ + 2 MP + 2 F + Omega] + (4) sin[ -2 D + 2 MP + Omega] + (4) sin[ -2 D + M + 2 F + Omega] + (4) sin[ + MP -2 F] + (-4) sin[ -1 D + MP] + (-4) sin[ -2 D + M] + (-4) sin[ + D] + (3) sin[ + MP + 2 F] + (-3) sin[ -2 MP + 2 F + 2 Omega] + (-3) sin[ -1 D -1 M + MP] + (-3) sin[ + M + MP] + (-3) sin[ -1 M + MP + 2 F + 2 Omega] + (-3) sin[ + 2 D -1 M -1 MP + 2 F + 2 Omega] + (-3) sin[ + 3 MP + 2 F + 2 Omega] + (-3) sin[ + 2 D -1 M + 2 F + 2 Omega])
deltaepsilon = 0.0001 arcsec ( + (92025 + 8.9 T) cos[ + Omega] + (5736 + -3.1 T) cos[ -2 D + 2 F + 2 Omega] + (977 + -0.5 T) cos[ + 2 F + 2 Omega] + (-895 + 0.5 T) cos[ + 2 Omega] + (54 + -0.1 T) cos[ + M] + (-7) cos[ + MP] + (224 + -0.6 T) cos[ -2 D + M + 2 F + 2 Omega] + (200) cos[ + 2 F + Omega] + (129 + -0.1 T) cos[ + MP + 2 F + 2 Omega] + (-95 + 0.3 T) cos[ -2 D -1 M + 2 F + 2 Omega] + (-70) cos[ -2 D + 2 F + Omega] + (-53) cos[ -1 MP + 2 F + 2 Omega] + (-33) cos[ + MP + Omega] + (26) cos[ + 2 D -1 MP + 2 F + 2 Omega] + (32) cos[ -1 MP + Omega] + (27) cos[ + MP + 2 F + Omega] + (-24) cos[ -2 MP + 2 F + Omega] + (16) cos[ + 2 D + 2 F + 2 Omega] + (13) cos[ + 2 MP + 2 F + 2 Omega] + (-12) cos[ -2 D + MP + 2 F + 2 Omega] + (-10) cos[ -1 MP + 2 F + Omega] + (-8) cos[ + 2 D -1 MP + Omega] + (7) cos[ -2 D + 2 M + 2 F + 2 Omega] + (9) cos[ + M + Omega] + (7) cos[ -2 D + MP + Omega] + (6) cos[ -1 M + Omega] + (5) cos[ + 2 D -1 MP + 2 F + Omega] + (3) cos[ + 2 D + MP + 2 F + 2 Omega] + (-3) cos[ + M + 2 F + 2 Omega] + (3) cos[ -1 M + 2 F + 2 Omega] + (3) cos[ + 2 D + 2 F + Omega] + (-3) cos[ -2 D + 2 MP + 2 F + 2 Omega] + (-3) cos[ -2 D + MP + 2 F + Omega] + (3) cos[ + 2 D -2 MP + Omega] + (3) cos[ + 2 D + Omega] + (3) cos[ -2 D -1 M + 2 F + Omega] + (3) cos[ -2 D + Omega] + (3) cos[ + 2 MP + 2 F + Omega])
return[deltapsi, deltaepsilon]
}
// Calculate the right ascension and declination of an object given ecliptical
// coordinates. This is an implementation of eq. 13.3 and 13.4 in Meeus.
//
// arguments:
// [lambda, beta, epsilon]
//
// returns:
// [ra, decl]
eclipticalToRADecl[lambda, beta, epsilon] :=
{
alpha = arctan[sin[lambda] cos[epsilon] - tan[beta] sin[epsilon], cos[lambda]]
delta = arcsin[sin[beta] cos[epsilon] + cos[beta] sin[epsilon] sin[lambda]]
return [alpha, delta]
}
// Calculate the azimuth and altitude of an object.
// Arguments:
// [date, objectRA, objectDecl, lat, long]
//
// date: time of observation
// objectRA: Right Ascension of object (can be angle or time)
// objectDecl: Declination of object
// lat: Latitude of observer (North is positive)
// long: Longitude of observer (West is positive)
//
// Returns:
// [azimuth, altitude] as angles. Azimuth is calculated as angle west from
// south. To get compass bearing, use
// (azimuth + 180 degrees) mod circle
//
// Implements equations 13.5 and 13.6 from Meeus.
// NOTE: This is not corrected for refraction nor parallax! These are
// geocentric coordinates, that is, relative to the center of the earth.
// Due to parallax, the altitude figure may be significantly in error for
// nearby bodies like the moon. To get the corrected altitudes, use the
// refractedAzimuthAltitude function below.
airlessAzimuthAltitude[date, objectRA, objectDecl, lat, long] :=
{
// We want all of this in angles... correct if necessary.
if (objectRA conforms hour)
objectRA = objectRA (circle/day)
H = hourAngle[date, long, objectRA]
azimuth = arctan[sin[H], cos[H] sin[lat] - tan[objectDecl] cos[lat]]
altitude = arcsin[sin[lat] sin[objectDecl] + cos[lat] cos[objectDecl] cos[H]]
return [azimuth, altitude]
}
// Calculate the apparent azimuth and altitude of an object with atmospheric
// refraction and parallax corrections.
// Arguments:
// [date, objectRA, objectDecl, lat, long, distance, temp, pressure]
//
// date: time of observation
// objectRA: Right Ascension of object (can be angle or time)
// objectDecl: Declination of object
// lat: Latitude of observer (North is positive)
// long: Longitude of observer (West is positive)
// distance: Distance to object (undef means don't correct for parallax)
// temp: Temperature of atmosphere
// pressure: Atmospheric pressure
//
// Returns:
// [azimuth, altitude] as angles. Azimuth is calculated as angle west from
// south. To get compass bearing, use
// (azimuth + 180 degrees) mod circle
//
// Implements equations 13.5 and 13.6 from Meeus.
//
refractedAzimuthAltitude[date, objectRA, objectDecl, lat, long, distance=undef, temp=283K, pressure=1010 millibars] :=
{
[azimuth, altitude] = airlessAzimuthAltitude[date, objectRA, objectDecl, lat, long]
// Correct for parallax is distance is defined.
if distance != undef
altitude = altitude - parallaxAngleAlt[distance, altitude]
// println["Parallax angle is " + (parallaxAngleAlt[distance, altitude] -> "degrees")]
corrAltitude = altitude + refractionAngle[altitude, temp, pressure]
return[azimuth, corrAltitude]
}
// Calculate angle of refraction given the altitude in airless coordinates.
// This is equation 16.4 in Meeus, although the equation in the book
// is explained ridiculously poorly, and I had to dig through old
// Sky & Telescope program listings to figure out how to apply the equation.
// This is the angle to be *added* to the airless altitude to get the
// refracted altitude.
refractionAngle[altitude, temp = 283 K, pressure = 1010 millibars] :=
{
// h5 is altitude in degrees
h5 = altitude/degree
v5 = (h5 + 10.3/(h5+5.11)) degree
pressCorr = pressure/(1010 millibars) (283 K)/temp
R = pressCorr 1.02 arcmin / tan[v5]
return R
}
// Returns the apparent azimuth and altitude for the sun, adjusted for
// nutation, but not adjusted for refraction.
//
// Returns:
// [azimuth, altitude] Azimuth is calculated as angle west from
// south. To get compass bearing, use
// (azimuth + 180 degrees) mod circle
airlessSunAzimuthAltitude[date, lat, long] :=
{
[sunRA, sunDecl] = sunApparentRADecl[date]
return airlessAzimuthAltitude[date, sunRA, sunDecl, lat, long]
}
// Returns the apparent azimuth and altitude for the sun, adjusted for
// nutation, atmospheric refraction, and parallax.
//
// Returns:
// [azimuth, altitude] as angles. Azimuth is calculated as angle west from
// south. To get compass bearing, use
// (azimuth + 180 degrees) mod circle
refractedSunAzimuthAltitude[date, lat, long, temp = 283K, pressure = 1010 millibars] :=
{
[sunRA, sunDecl] = sunApparentRADecl[date]
return refractedAzimuthAltitude[date, sunRA, sunDecl, lat, long, sunDistance[date], temp, pressure]
}
// Secant method solver to find when the sun is at a given altitude to
// an observer at the specified latitude and longitude on the Earth.
// It takes a first guess at date1 and adds an hour to that to find
// its two starting points. Thus, to find sunrise, give it a guess
// of sometime around 6 AM local time and a desired altitude of 0 degrees,
// or about -16 arcmin to account for the time that the sun's upper limb
// is just crosses the horizon. To find sunset, give it a guess
// around 6 PM local time.
//
// This currently doesn't handle the case where the sun doesn't cross through
// the specified azimuth on that date and may give unpredictable results in
// this case or in the case when the initial guess is far off.
sunSecantAltitude[date1, lat, long, desiredAltitude, temperature = 283 K,
pressure=1010 millibars ] :=
{
date2 = date1 + 1 hour
[azimuth1, altitude1] = refractedSunAzimuthAltitude[date1, lat, long, temperature, pressure]
[azimuth2, altitude2] = refractedSunAzimuthAltitude[date2, lat, long, temperature, pressure]
altitude1 = altitude1 - desiredAltitude
altitude2 = altitude2 - desiredAltitude
// println["Altitude 1: $altitude1"]
// println["Altitude 2: $altitude2"]
while (true)
{
correction = (altitude1 * (date1 - date2)) / (altitude1-altitude2)
date = date1 - correction
// println[date -> ### HH:mm:ss.SSS ###]
if abs[correction] < 1 second
return date
date2 = date1
date1 = date
altitude2 = altitude1
[azimuth1, altitude1] = refractedSunAzimuthAltitude[date, lat, long, temperature, pressure]
altitude1 = altitude1 - desiredAltitude
}
}
// Secant method solver to find when the sun is at a given altitude to
// an observer at the specified latitude and longitude on the Earth.
// This version does *not* correct for refraction. To correct for refraction,
// use the function sunSecantAltitude[] instead.
// It takes a first guess at date1 and adds an hour to that to find
// its two starting points. Thus, to find sunrise, give it a guess
// of sometime around 6 AM local time and a desired altitude of 0 degrees,
// or about -16 arcmin to account for the time that the sun's upper limb
// is just crosses the horizon. To find sunset, give it a guess
// around 6 PM local time.
//
// This currently doesn't handle the case where the sun doesn't cross through
// the specified azimuth on that date and may give unpredictable results in
// this case or in the case when the initial guess is far off.
sunSecantAirlessAltitude[date1, lat, long, desiredAltitude ] :=
{
date2 = date1 + 1 hour
[azimuth1, altitude1] = airlessSunAzimuthAltitude[date1, lat, long]
[azimuth2, altitude2] = airlessSunAzimuthAltitude[date2, lat, long]
altitude1 = altitude1 - desiredAltitude
altitude2 = altitude2 - desiredAltitude
// println["Altitude 1: $altitude1"]
// println["Altitude 2: $altitude2"]
while (true)
{
correction = (altitude1 * (date1 - date2)) / (altitude1-altitude2)
date = date1 - correction
// println[date -> ### HH:mm:ss.SSS ###]
if abs[correction] < 1 second
return date
date2 = date1
date1 = date
altitude2 = altitude1
[azimuth1, altitude1] = airlessSunAzimuthAltitude[date, lat, long]
altitude1 = altitude1 - desiredAltitude
}
}
// Secant method solver to find when the sun is at a given azimuth to
// an observer at the specified latitude and longitude on the Earth.
// It takes a first guess at date1 and adds an hour to that to find
// its two starting points. Thus, to find sunrise, give it a guess
// of sometime around 6 AM local time and a desired altitude of 0 degrees,
// or about -16 arcmin to account for the time that the sun's upper limb
// is just crosses the horizon. To find sunset, give it a guess
// around 6 PM local time.
//
// This function corrects for refraction.
//
// This currently doesn't handle the case where the sun doesn't pass through
// the desired azimuth on that date and may give unpredictable results in
// this case or the case where the initial guess is far off.
//
// desiredAzimuth should be specified in the coordinate system specified by
// Meeus, which is angle west of south. To convert from normal true compass
// bearings, use (trueAzimuth + 180 degrees) mod circle
sunSecantAzimuth[date1, lat, long, desiredAzimuth, temperature = 283 K,
pressure=1010 millibars ] :=
{
date2 = date1 + 5 minutes
[azimuth1, altitude1] = refractedSunAzimuthAltitude[date1, lat, long, temperature, pressure]
[azimuth2, altitude2] = refractedSunAzimuthAltitude[date2, lat, long, temperature, pressure]
azimuth1 = azimuth1 - desiredAzimuth
azimuth2 = azimuth2 - desiredAzimuth
while (true)
{
// Don't overcorrect... if we're correcting more than half a phase
// from the initial guess, then go back a phase.
if azimuth1 > 180 degrees
azimuth1 = azimuth1 - circle
if azimuth1 < -180 degrees
azimuth1 = azimuth1 + circle
if azimuth2 > 180 degrees
azimuth2 = azimuth2 - circle
if azimuth2 < -180 degrees
azimuth2 = azimuth2 + circle
// did we wrap around the circle?
if abs[azimuth1 - azimuth2] > 180 degrees
if (azimuth1 < azimuth2)
azimuth2 = azimuth2 - circle
else
azimuth1 = azimuth1 - circle
// println["Azimuth 1: $azimuth1"]
// println["Azimuth 2: $azimuth2"]
correction = (azimuth1 * (date1 - date2)) / (azimuth1-azimuth2)
// println["Correction: $correction"]
date = date1 - correction
// println[date]
if abs[correction] < 1 second
return date
date2 = date1
date1 = date
azimuth2 = azimuth1
[azimuth1, altitude1] = refractedSunAzimuthAltitude[date, lat, long, temperature, pressure]
azimuth1 = azimuth1 - desiredAzimuth
}
}
// Secant method solver to find when the moon is at a given altitude to
// an observer at the specified latitude and longitude on the Earth.
// It takes a first guess at date1 and adds an hour to that to find
// its two starting points.
//
// This currently doesn't handle the case where the moon doesn't cross through
// the specified altitude on that date, or if initial guess is very wrong,
// and may give unpredictable results.
moonSecantAltitude[date1, lat, long, desiredAltitude, temperature = 283 K,
pressure=1010 millibars ] :=
{
date2 = date1 + 1 hour
[azimuth1, altitude1] = refractedMoonAzimuthAltitude[date1, lat, long, temperature, pressure]
[azimuth2, altitude2] = refractedMoonAzimuthAltitude[date2, lat, long, temperature, pressure]
altitude1 = altitude1 - desiredAltitude
altitude2 = altitude2 - desiredAltitude
// println["Altitude 1: $altitude1"]
// println["Altitude 2: $altitude2"]
while (true)
{
correction = (altitude1 * (date1 - date2)) / (altitude1-altitude2)
if (correction > 24 hours)
correction = 24 hours
if (correction < -24 hours)
correction = -24 hours
date = date1 - correction
// println[date -> ### HH:mm:ss.SSS ###]
// println[date]
if abs[correction] < 1 second
return date
date2 = date1
date1 = date
altitude2 = altitude1
[azimuth1, altitude1] = refractedMoonAzimuthAltitude[date, lat, long, temperature, pressure]
altitude1 = altitude1 - desiredAltitude
}
}
// Calculate the approximate mean time of various phases of the moon.
// This will give the approximate time that the specified phase
// will occur near the given date. This is only accurate to within
// some number of hours.
//
// The phaseAdjust can be only one of 4 values:
// 0 for New Moon
// 1/4 for First Quarter
// 1/2 for Full Moon
// 3/4 for Last Quarter
//
// All other values give meaningless results!
//
// This guess serves to supply a good starting guess to moonPhase
// which is much higher accuracy.
moonApproximatePhase[date, phaseAdjust] :=
{
// This was found experimentally by Alan Eliasen to give a better initial
// estimate that centers around the closest phase of the moon, with a
// slight bias of about 4-5 days toward selecting the "next" event.
d1 = date - phaseAdjust lunarmonth
// Get approximate k, Meeus eq. 49.2
k = round[(JDE[d1] - JDE[#2000-01-01#])/year * 12.3685] + phaseAdjust
T = k / 1236.85
approxDate = JDE[2451550.09766 + 29.530588861 k + 0.00015437 T^2 -
0.000000150 T^3 + 0.00000000073 T^4]
return approxDate
}
// This finds the excess of the apparent geocentric latitude of the moon
// over the apparent geocentric longitude of the sun. It is used to find
// the phases of the moon to medium accuracy. The accuracy of this is limited
// to the somewhat low accuracy calculation of sunApparentLongitude.
// This could be improved by using the high-accuracy VSOP87 theory
// implemented in planets.frink and Meeus p. 166
// See Meeus chapter 49.
moonExcessLongitude[date] := (moonApparentLongitude[date] - sunApparentLongitude[date]) mod circle
// Find high accuracy phases of the moon.
// The phaseAdjust can be only one of 4 values:
// 0 for New Moon
// 1/4 for First Quarter
// 1/2 for Full Moon
// 3/4 for Last Quarter
//
// All other values give meaningless results!
//
// For simplicity, you should probably use one of the convenience functions:
// newMoon[date], fullMoon[date], etc.
//
// This function uses a secant method successive approximation to find the
// moment of the desired phase to high accuracy using the full-precision
// moon calculations. The convergence limit is set at one second or less.
moonPhase[date, phaseAdjust] :=
{
desiredAngle = phaseAdjust * circle
date1 = moonApproximatePhase[date, phaseAdjust]
// println[date1]
date2 = date1 + 5 min
excess1 = moonExcessLongitude[date1]
excess2 = moonExcessLongitude[date2]
// println["$excess1\t$excess2"]
// println["Phase is " + excess1 / circle]
// Don't overcorrect... if we're correcting more than half a phase forward
// from the initial guess, then go back a phase.
// This generally corrects the problem around the new moon where the angle
// jumps from close to 360 degrees back to 0 degrees, and may prevent other
// as-yet-undetected weird jumps.
if (excess1 - desiredAngle > 180 degrees)
excess1 = excess1 - circle
if (excess2 - desiredAngle > 180 degrees)
excess2 = excess2 - circle
// did we wrap around the circle?
if abs[excess1 - excess2] > 180 degrees
if (excess1 < excess2)
excess2 = excess2 - circle
else
excess1 = excess1 - circle
// println["$excess1\t$excess2"]
excess1 = excess1 - desiredAngle
excess2 = excess2 - desiredAngle
// println["$excess1\t$excess2"]
while (true)
{
correction = (excess1 * (date1 - date2)) / (excess1-excess2)
date = date1 - correction
// println[date]
// println["$excess1\t$excess2"]
if abs[correction] < 1 second
return date
date2 = date1
date1 = date
excess2 = excess1
excess1 = moonExcessLongitude[date1]
excess1 = excess1 - desiredAngle
// did we wrap around the circle?
if abs[excess1 - excess2] > 180 degrees
if (excess1 < excess2)
excess2 = excess2 - circle
else
excess1 = excess1 - circle
}
return date
}
// Finds the date of the New Moon close to the specifed date.
newMoon[date] := moonPhase[date, 0]
// Finds the date of the First Quarter moon close to the specifed date.
firstQuarterMoon[date] := moonPhase[date, 1/4]
// Finds the date of the Full Moon close to the specifed date.
fullMoon[date] := moonPhase[date, 1/2]
// Finds the date of the Last Quarter moon close to the specifed date.
lastQuarterMoon[date] := moonPhase[date, 3/4]
// Secant method solver to find when the moon is at a given azimuth to
// an observer at the specified latitude and longitude on the Earth.
//
// This currently doesn't handle the case where the moon doesn't pass through
// a certain azimuth on that date and may give unpredictable results in that
// case or the case where the initial guess is far off.
//
// desiredAzimuth should be specified in the coordinate system specified by
// Meeus, which is angle west of south. To convert from normal true compass
// bearings, use (trueAzimuth + 180 degrees) mod circle
moonSecantAzimuth[date1, lat, long, desiredAzimuth, temperature = 283 K,
pressure=1010 millibars ] :=
{
date2 = date1 + 5 minutes
[azimuth1, altitude1] = refractedMoonAzimuthAltitude[date1, lat, long, temperature, pressure]
[azimuth2, altitude2] = refractedMoonAzimuthAltitude[date2, lat, long, temperature, pressure]
azimuth1 = azimuth1 - desiredAzimuth
azimuth2 = azimuth2 - desiredAzimuth
// println["Azimuth 1: " + (azimuth1->"degrees")]
// println["Azimuth 2: " + (azimuth2->"degrees")]
while (true)
{
// Don't overcorrect... if we're correcting more than half a phase
// from the initial guess, then go back a phase.
if azimuth1 > 180 degrees
azimuth1 = azimuth1 - circle
if azimuth1 < -180 degrees
azimuth1 = azimuth1 + circle
if azimuth2 > 180 degrees
azimuth2 = azimuth2 - circle
if azimuth2 < -180 degrees
azimuth2 = azimuth2 + circle
// did we wrap around the circle?
if abs[azimuth1 - azimuth2] > 180 degrees
if (azimuth1 < azimuth2)
azimuth2 = azimuth2 - circle
else
azimuth1 = azimuth1 - circle
// println["Diff is " + (azimuth1-azimuth2 -> "degrees")]
correction = (azimuth1 * (date1 - date2)) / (azimuth1-azimuth2)
// println["Correction: $correction"]
date = date1 - correction
// println[date]
if abs[correction] < 1 second
return date
date2 = date1
date1 = date
azimuth2 = azimuth1
[azimuth1, altitude1] = refractedMoonAzimuthAltitude[date, lat, long, temperature, pressure]
azimuth1 = azimuth1 - desiredAzimuth
}
}
// Calculate the time when the upper limb of the sun just crosses the
// horizon.
sunrise[date, lat, long, temp = 283 K, pressure=1010 millibars ] := sunSecantAltitude[approxMidnight[date, long] + 5.5 hours,
lat, long, -16 arcmin, temp, pressure]
// Calculate the time when the upper limb of the sun just crosses the
// horizon.
sunset[date, lat, long, temp = 283 K, pressure=1010 millibars ] := sunSecantAltitude[approxMidnight[date, long] + 17.5 hours,
lat, long, -16 arcmin, temp, pressure]
civilTwilightBegin[date, lat, long] := sunSecantAirlessAltitude[approxMidnight[date, long] + 5 hour, lat, long, -6 degree]
civilTwilightEnd[date, lat, long] := sunSecantAirlessAltitude[approxMidnight[date, long] + 18 hour, lat, long, -6 degree]
nauticalTwilightBegin[date, lat, long] := sunSecantAirlessAltitude[approxMidnight[date, long] + 4.7 hour, lat, long, -12 degree]
nauticalTwilightEnd[date, lat, long] := sunSecantAirlessAltitude[approxMidnight[date, long] + 18.3 hour, lat, long, -12 degree]
astronomicalTwilightBegin[date, lat, long] := sunSecantAirlessAltitude[approxMidnight[date, long] + 4.3 hour, lat, long, -18 degree]
astronomicalTwilightEnd[date, lat, long] := sunSecantAirlessAltitude[approxMidnight[date, long] + 18.7 hour, lat, long, -18 degree]
// Get a time object representing the approximate midnight that begins the
// the specified date at the specified longitude. It is highly recommended
// that you use a date significantly after midnight to eliminate issues with
// the local meridian and daylight savings time.
approxMidnight[date, longitude] :=
{
d0 = #2003-01-01 UTC# + (longitude / (circle/day))
r = d0 + floor[(date - d0) / day] day
// println[r]
return r
}
// Besselian Date, as defined by http://maia.usno.navy.mil/
BesselianDate[date] := 2000 + (MJD[date]/day - 51544.03) / 365.2422
//
// Moon position calculations
//
// Moon mean longitude L', referred to the mean equinox of date, including the
// constant term of the effect of the light-time. This is equation 47.1
// in Meeus.
moonMeanLongitude[date] :=
{
T = meeusT[date]
return ((218.3164477 + 481267.88123421 T - 0.0015786 T^2 + T^3/538_841 -
T^4/65_194_000) degree) mod circle
}
// Moon mean elongation as an angle D. This is equation 47.2 in Meeus.
moonMeanElongation[date] :=
{
T = meeusT[date]
return ((297.8501921 + 445267.1114034 T - 0.0018819 T^2 + T^3/545_868 -
T^4/113_065_000) degree) mod circle
}
// Calculate the mean anomaly of the sun, M. This is equation 47.3 in Meeus.
// Note that these numbers are very slightly different than the sun anomaly
// function given by equation 25.3 and called sunMeanAnomaly above.
moonCalcSunMeanAnomaly[date] :=
{
T = meeusT[date]
return ((357.5291092 + 35999.0502909 T - 0.0001536 T^2 +
T^3/24_490_000) degree) mod circle
}
// Moon mean anomaly M'. This is equation 47.4 in Meeus.
moonMeanAnomaly[date] :=
{
T = meeusT[date]
return ((134.9633964 + 477198.8675055 T + 0.0087414 T^2 + T^3/69_699 -
T^4 / 14_712_000) degree) mod circle
}
// Moon's argument of latitude F, the mean distance of the Moon from its
// ascending node. This is equation 47.5 in Meeus.
moonArgumentOfLatitude[date] :=
{
T = meeusT[date]
return ((93.2720950 + 483202.0175233 T - 0.0036539 T^2 - T^3 / 3_526_000 +
T^4 / 863_310_000) degree) mod circle
}
// Moon correction term A1 (listed after equation 47.5 in Meeus)
// This is due to the action of Venus.
moonA1[date] :=
{
T = meeusT[date]
return ((119.75 + 131.849 T) degree) mod circle
}
// Moon correction term A2 (listed after equation 47.5 in Meeus)
// This is due to the action of Jupiter.
moonA2[date] :=
{
T = meeusT[date]
return ((53.09 + 479264.290 T) degree) mod circle
}
// Moon correction term A3 (listed after equation 47.5 in Meeus)
moonA3[date] :=
{
T = meeusT[date]
return ((313.45 + 481266.484 T) degree) mod circle
}
// Earth eccentricity E for Moon calculations. This is equation 47.6 in
// Meeus. This is not the eccentricity given by eq. 25.4
moonCalcEarthEccentricity[date] :=
{
T = meeusT[date]
return 1 - 0.002516 T - 0.0000074 T^2
}
// Moon sum of L from Meeus, chapter 47. This is with the extra additive
// terms listed on p. 342.
moonSumL[date] :=
{
D = moonMeanElongation[date]
M = moonCalcSunMeanAnomaly[date]
MP = moonMeanAnomaly[date]
F = moonArgumentOfLatitude[date]
E = moonCalcEarthEccentricity[date]
LP = moonMeanLongitude[date]
sumL = ( + (6288774) sin[ + MP] + (1274027) sin[ + 2 D - MP] + (658314) sin[ + 2 D] + (213618) sin[ + 2 MP] + (-185116 E) sin[ + M] + (-114332) sin[ + 2 F] + (58793) sin[ + 2 D -2 MP] + (57066 E) sin[ + 2 D - M - MP] + (53322) sin[ + 2 D + MP] + (45758 E) sin[ + 2 D - M] + (-40923 E) sin[ + M - MP] + (-34720) sin[ + D] + (-30383 E) sin[ + M + MP] + (15327) sin[ + 2 D -2 F] + (-12528) sin[ + MP + 2 F] + (10980) sin[ + MP -2 F] + (10675) sin[ + 4 D - MP] + (10034) sin[ + 3 MP] + (8548) sin[ + 4 D -2 MP] + (-7888 E) sin[ + 2 D + M - MP] + (-6766 E) sin[ + 2 D + M] + (-5163) sin[ + D - MP] + (4987 E) sin[ + D + M] + (4036 E) sin[ + 2 D - M + MP] + (3994) sin[ + 2 D + 2 MP] + (3861) sin[ + 4 D] + (3665) sin[ + 2 D -3 MP] + (-2689 E) sin[ + M -2 MP] + (-2602) sin[ + 2 D - MP + 2 F] + (2390 E) sin[ + 2 D - M -2 MP] + (-2348) sin[ + D + MP] + (2236 E^2) sin[ + 2 D -2 M] + (-2120 E) sin[ + M + 2 MP] + (-2069 E^2) sin[ + 2 M] + (2048 E^2) sin[ + 2 D -2 M - MP] + (-1773) sin[ + 2 D + MP -2 F] + (-1595) sin[ + 2 D + 2 F] + (1215 E) sin[ + 4 D - M - MP] + (-1110) sin[ + 2 MP + 2 F] + (-892) sin[ + 3 D - MP] + (-810 E) sin[ + 2 D + M + MP] + (759 E) sin[ + 4 D - M -2 MP] + (-713 E^2) sin[ + 2 M - MP] + (-700 E^2) sin[ + 2 D + 2 M - MP] + (691 E) sin[ + 2 D + M -2 MP] + (596 E) sin[ + 2 D - M -2 F] + (549) sin[ + 4 D + MP] + (537) sin[ + 4 MP] + (520 E) sin[ + 4 D - M] + (-487) sin[ + D -2 MP] + (-399 E) sin[ + 2 D + M -2 F] + (-381) sin[ + 2 MP -2 F] + (351 E) sin[ + D + M + MP] + (-340) sin[ + 3 D -2 MP] + (330) sin[ + 4 D -3 MP] + (327 E) sin[ + 2 D - M + 2 MP] + (-323 E^2) sin[ + 2 M + MP] + (299 E) sin[ + D + M - MP] + (294) sin[ + 2 D + 3 MP])
sumL = sumL + 3958 sin[moonA1[date]] +
1962 sin[LP - F] +
318 sin[moonA2[date]]
return sumL
}
// Moon longitude calculation from section 47 of Meeus. This is not corrected
// for Earth's nutation. For that, use moonApparentLongitude[date]
moonLongitude[date] :=
{
return moonMeanLongitude[date] + (moonSumL[date]/1000000 degree)
}
// Returns the moon's apparent longitude, corrected for Earth's
// nutation. It is not corrected for refraction.
moonApparentLongitude[date] :=
{
[deltaPsi, deltaEpsilon] = highAccuracyNutation[date]
return moonLongitude[date] + deltaPsi
}
// Moon distance calculation Delta from section 47 of Meeus.
moonDistance[date] :=
{
D = moonMeanElongation[date]
M = moonCalcSunMeanAnomaly[date]
MP = moonMeanAnomaly[date]
F = moonArgumentOfLatitude[date]
E = moonCalcEarthEccentricity[date]
LP = moonMeanLongitude[date]
sumR = ( + (-20905355) cos[ + MP] + (-3699111) cos[ + 2 D - MP] + (-2955968) cos[ + 2 D] + (-569925) cos[ + 2 MP] + (48888 E) cos[ + M] + (-3149) cos[ + 2 F] + (246158) cos[ + 2 D -2 MP] + (-152138 E) cos[ + 2 D - M - MP] + (-170733) cos[ + 2 D + MP] + (-204586 E) cos[ + 2 D - M] + (-129620 E) cos[ + M - MP] + (108743) cos[ + D] + (104755 E) cos[ + M + MP] + (10321) cos[ + 2 D -2 F] + (79661) cos[ + MP -2 F] + (-34782) cos[ + 4 D - MP] + (-23210) cos[ + 3 MP] + (-21636) cos[ + 4 D -2 MP] + (24208 E) cos[ + 2 D + M - MP] + (30824 E) cos[ + 2 D + M] + (-8379) cos[ + D - MP] + (-16675 E) cos[ + D + M] + (-12831 E) cos[ + 2 D - M + MP] + (-10445) cos[ + 2 D + 2 MP] + (-11650) cos[ + 4 D] + (14403) cos[ + 2 D -3 MP] + (-7003 E) cos[ + M -2 MP] + (10056 E) cos[ + 2 D - M -2 MP] + (6322) cos[ + D + MP] + (-9884 E^2) cos[ + 2 D -2 M] + (5751 E) cos[ + M + 2 MP] + (-4950 E^2) cos[ + 2 D -2 M - MP] + (4130) cos[ + 2 D + MP -2 F] + (-3958 E) cos[ + 4 D - M - MP] + (3258) cos[ + 3 D - MP] + (2616 E) cos[ + 2 D + M + MP] + (-1897 E) cos[ + 4 D - M -2 MP] + (-2117 E^2) cos[ + 2 M - MP] + (2354 E^2) cos[ + 2 D + 2 M - MP] + (-1423) cos[ + 4 D + MP] + (-1117) cos[ + 4 MP] + (-1571 E) cos[ + 4 D - M] + (-1739) cos[ + D -2 MP] + (-4421) cos[ + 2 MP -2 F] + (1165 E^2) cos[ + 2 M + MP] + (8752) cos[ + 2 D - MP -2 F])
return (385000560 + sumR) meter
}
// Moon latitude calculation beta from section 47 of Meeus.
moonLatitude[date] :=
{
D = moonMeanElongation[date]
M = moonCalcSunMeanAnomaly[date]
MP = moonMeanAnomaly[date]
F = moonArgumentOfLatitude[date]
E = moonCalcEarthEccentricity[date]
LP = moonMeanLongitude[date]
A1 = moonA1[date]
sumB = ( + (5128122) sin[ + F] + (280602) sin[ + MP + F] + (277693) sin[ + MP - F] + (173237) sin[ + 2 D - F] + (55413) sin[ + 2 D - MP + F] + (46271) sin[ + 2 D - MP - F] + (32573) sin[ + 2 D + F] + (17198) sin[ + 2 MP + F] + (9266) sin[ + 2 D + MP - F] + (8822) sin[ + 2 MP - F] + (8216 E) sin[ + 2 D - M - F] + (4324) sin[ + 2 D -2 MP - F] + (4200) sin[ + 2 D + MP + F] + (-3359 E) sin[ + 2 D + M - F] + (2463 E) sin[ + 2 D - M - MP + F] + (2211 E) sin[ + 2 D - M + F] + (2065 E) sin[ + 2 D - M - MP - F] + (-1870 E) sin[ + M - MP - F] + (1828) sin[ + 4 D - MP - F] + (-1794 E) sin[ + M + F] + (-1749) sin[ + 3 F] + (-1565 E) sin[ + M - MP + F] + (-1491) sin[ + D + F] + (-1475 E) sin[ + M + MP + F] + (-1410 E) sin[ + M + MP - F] + (-1344 E) sin[ + M - F] + (-1335) sin[ + D - F] + (1107) sin[ + 3 MP + F] + (1021) sin[ + 4 D - F] + (833) sin[ + 4 D - MP + F] + (777) sin[ + MP -3 F] + (671) sin[ + 4 D -2 MP + F] + (607) sin[ + 2 D -3 F] + (596) sin[ + 2 D + 2 MP - F] + (491 E) sin[ + 2 D - M + MP - F] + (-451) sin[ + 2 D -2 MP + F] + (439) sin[ + 3 MP - F] + (422) sin[ + 2 D + 2 MP + F] + (421) sin[ + 2 D -3 MP - F] + (-366 E) sin[ + 2 D + M - MP + F] + (-351 E) sin[ + 2 D + M + F] + (331) sin[ + 4 D + F] + (315 E) sin[ + 2 D - M + MP + F] + (302 E^2) sin[ + 2 D -2 M - F] + (-283) sin[ + MP + 3 F] + (-229 E) sin[ + 2 D + M + MP - F] + (223 E) sin[ + D + M - F] + (223 E) sin[ + D + M + F] + (-220 E) sin[ + M -2 MP - F] + (-220 E) sin[ + 2 D + M - MP - F] + (-185) sin[ + D + MP + F] + (181 E) sin[ + 2 D - M -2 MP - F] + (-177 E) sin[ + M + 2 MP + F] + (176) sin[ + 4 D -2 MP - F] + (166 E) sin[ + 4 D - M - MP - F] + (-164) sin[ + D + MP - F] + (132) sin[ + 4 D + MP - F] + (-119) sin[ + D - MP - F] + (115 E) sin[ + 4 D - M - F] + (107 E^2) sin[ + 2 D -2 M + F])
sumB = sumB - 2235 sin[LP] + 382 sin[moonA3[date]] + 175 sin[A1 - F] +
175 sin[A1 + F] + 127 sin[LP - MP] - 115 sin[LP + MP]
return (sumB / 1000000) degree
}
// Calculates the apparent Right ascension and declination of the moon
// at a given date.
// returns [ra, decl]
moonApparentRADecl[date] :=
{
epsilon = trueObliquityOfEcliptic[date]
moonLong = moonApparentLongitude[date]
moonLat = moonLatitude[date]
se = sin[epsilon]
ce = cos[epsilon]
sm = sin[moonLong]
ra = arctan[sm ce - tan[moonLat] se,
cos[moonLong]]
decl = arcsin[sin[moonLat] ce + cos[moonLat] sin[epsilon] sm]
return[ra,decl]
}
// Returns the apparent azimuth and altitude for the moon, adjusted for
// nutation, but not adjusted for refraction.
//
// Returns:
// [azimuth, altitude] Azimuth is calculated as angle west from
// south. To get compass bearing, use
// (azimuth + 180 degrees) mod circle
airlessMoonAzimuthAltitude[date, lat, long] :=
{
[moonRA, moonDecl] = moonApparentRADecl[date]
return airlessAzimuthAltitude[date, moonRA, moonDecl, lat, long]
}
// Returns the apparent azimuth and altitude for the moon, adjusted for
// nutation, atmospheric refraction, and parallax.
//
// Returns:
// [azimuth, altitude] as angles. Azimuth is calculated as angle west from
// south. To get compass bearing, use
// (azimuth + 180 degrees) mod circle
refractedMoonAzimuthAltitude[date, lat, long, temp = 283K, pressure=1010 millibars] :=
{
[moonRA, moonDecl] = moonApparentRADecl[date]
return refractedAzimuthAltitude[date, moonRA, moonDecl, lat, long, moonDistance[date], temp, pressure]
}
// Calculates the apparent moon radius angle, corrected for distance, but not
// corrected for refraction. It gives a true arcsin geometric interpretation.
moonRadiusAngle[date] := arcsin[moonradius / (moonradius + moonDistance[date])]
// Calculates the geocentric elongation of the moon from the sun
// by Meeus equation 48.2 (second equation)
moonGeocentricElongationFromSun[date] :=
{
arccos[cos[moonLatitude[date]] * cos[moonApparentLongitude[date] - sunTrueLongitude[date]]]
}
// Calculates the geocentric elongation of the moon from the sun
// by Meeus equation 48.2 (first equation)
moonGeocentricElongationFromSun2[date] :=
{
[alpha0, delta0] = sunApparentRADecl[date]
[alpha, delta] = moonApparentRADecl[date]
return arccos[sin[delta0] sin[delta] +
cos[delta0] cos[delta] cos[alpha0 - alpha]]
}
// Phase angle of the moon as seen from earth, by equation 48.3 of Meeus
// This is constrained to be in the range 0-180 degrees as expected by the
// illuminated fraction calculations.
moonPhaseAngle[date] :=
{
R = sunDistance[date]
delta = moonDistance[date]
phi = moonGeocentricElongationFromSun[date]
return arctan[R sin[phi] / (delta - R cos[phi])] mod (180 degrees)
}
// Illuminated fraction of the moon as seen from earth. This is found by
// equation 48.1 of Meeus.
moonIlluminatedFraction[date] := (1 + cos[moonPhaseAngle[date]]) / 2
// The position angle of the moon's bright limb, reckoned eastward from the
// North Point of the moon. (That is, relative to the earth's North/South
// axis projected onto the moon's disc, *not* relative to the moon's axis.)
// This is equation 48.5 in Meeus.
moonPositionAngle[date] :=
{
[alpha0, delta0] = sunApparentRADecl[date]
[alpha, delta] = moonApparentRADecl[date]
return arctan[cos[delta0] sin[alpha0 - alpha],
sin[delta0] cos[delta] - cos[delta0] sin[delta] cos[alpha0-alpha]]
}
// The position angle of the moon's bright limb, reckoned "eastward" (?)
// from the zenith. Note that this gives the direction *counterclockwise*
// from the zenith as seen from an observer standing on earth. The "eastward"
// terminology is confusing.
// This is equation 48.5 in Meeus.
moonPositionAngleRelativeToZenith[date, lat, long] :=
{
chi = moonPositionAngle[date]
[ra, decl] = moonApparentRADecl[date]
q = parallacticAngle[date, lat, long, ra, decl]
return chi - q
}
// Returns a GeneralPath that represents the position of the moon.
// arguments:
// cx, cy: centerpoint of moon
// radius: radius of moon
// angle: an angle, as obtained by one of the moonPositionAngle...
// functions above.
// illum: the illuminated fraction of the moon as obtained from the
// moonIlluminatedFraction function above.
// filled: A boolean value indicating if the polygon is filled or not.
drawMoonPolygon[cx, cy, radius, angle, illumFraction, filled] :=
{
// The shorter semidiameter of the inner ellipse.
re = radius * (2 illumFraction - 1)
ca = cos[angle]
sa = sin[angle]
if (filled)
gp = new filledGeneralPath
else
gp = new GeneralPath
// Draw outer limb
gp.moveTo[radius ca + cx, -radius sa + cy]
gp.circularArc[cx, cy, 180 degrees]
// println["illumFraction is $illumFraction"]
// println["re is $re"]
// println["angle is " + (angle->degrees)]
// Draw inner ellipse
for theta = -90 degrees to 90 degrees step 2 degrees
{
x = re cos[theta]
y = radius sin[theta]
yp = x ca - y sa + cy
xp = x sa + y ca + cx
gp.lineTo[xp,yp]
}
gp.close[]
return gp
}
// Draws a moon polygon relative to to the zenith.
// Params:
// date: The date to draw for
// lat, long: The latitude and longitude of the viewing location
// cx, cy: center coordinates to draw to
// radius: radius to draw moon
// filled: boolean flag describing whether to draw moon.
drawMoonPolygonRelativeToZenith[date, lat, long, cx, cy, radius, filled] :=
{
return drawMoonPolygon[cx, cy, radius, moonPositionAngleRelativeToZenith[date, lat, long], moonIlluminatedFraction[date], filled]
}
// Parallax angle corrections (in azimuth and altitude) of an object as seen
// from earth. This is needed because most of the other calculations here
// are geocentric, i.e., taken relative to the center of the earth. The
// parallax becomes especially significant when taken to objects like the
// moon which are close to earth, and when looking at them close to the
// horizon. Sigh. I found this out the hard way.
// RETURNS:
// altParallax
parallaxAngleAlt[distance, origAltitude] :=
{
// TODO: Add a term for altitude above the reference geoid
// Horizontal parallax, called pi in Meeus
// This is the equation *before* equation 40.1 in Meeus.
// azParallax = arcsin[8.794 arcsec/(distance/au)]
// The following possibly gives more accuracy and is more general
// without the random unexplained factors of Meeus.
azParallax = arcsin[earthradius_equatorial/distance]
// TODO: Correct this for non-sphericity of Earth. For now, rho
// is taken to be 1.
altParallax = arcsin[sin[azParallax] cos[origAltitude]]
return altParallax
}
// Returns the angular separation between two bodies given their right
// ascensions and declinations.
angularSeparation[ra1, decl1, ra2, decl2] :=
{
return arccos[sin[decl1] sin[decl2] + cos[decl1] cos[decl2] cos[ra1-ra2]]
}
// Convenience function to obtain the angular separation between the sun
// and moon for a given time.
sunMoonAngularSeparation[date] :=
{
[sunRA, sunDecl] = sunApparentRADecl[date]
[moonRA, moonDecl] = moonApparentRADecl[date]
return angularSeparation[sunRA, sunDecl, moonRA, moonDecl]
}
// Guess when the moon will be near the specified azimuth.
// The azimuth should be specified in real compass coordinates,
// *not* Meeus-like coordinates. Thus, East should be 90 degrees, West 270
// degrees. This gives a very sloppy guess which should not be used
// directly, but rather as a seed to another routine, like moonrise[] or
// moonset[]
guessMoonAzimuth[date, lat, long, desiredAzimuth] :=
{
[az, alt] = airlessMoonAzimuthAltitude[date, lat, long]
realaz = (az + 180 deg)
// println["Real azimuth is " + (realaz -> degrees)]
// println["Altitude is " + (alt -> degrees)]
azcorr = (realaz - desiredAzimuth) / (circle/ day)
if (azcorr > 12 hours)
azcorr = azcorr - day
if (azcorr < -12 hours)
azcorr = azcorr + day
// println["Azimuth correction is " + (azcorr -> "hours")]
timeguess = date - azcorr
// println["Time guess is $timeguess"]
return timeguess
}
// Calculate the time when the upper limb of the moon just crosses the
// horizon.
moonrise[date, lat, long, temp = 283 K, pressure=1010 millibars ] := moonSecantAltitude[guessMoonAzimuth[date, lat, long, 90 degrees],lat, long, -16 arcmin, temp, pressure]
// Calculate the time when the upper limb of the moon just crosses the
// horizon.
moonset[date, lat, long, temp = 283 K, pressure=1010 millibars ] := moonSecantAltitude[guessMoonAzimuth[date, lat, long, 270 degrees],lat, long, -16 arcmin, temp, pressure]
// Calculate the dip angle of the horizon. This is a simplified calculation
// that just takes the height above the average sphere.
horizonDip[height] := arccos[earthradius / (earthradius+height)]
"sun.frink included Ok"
View or download sun.frink in plain text format
This is a program written in the programming language Frink.
For more information, view the Frink
Documentation or see More Sample Frink Programs.
Alan Eliasen was born 14704 days, 20 hours, 8 minutes ago.