# sun.frink

``` // 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 earths 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 } // Returns the longitude of the perihelion of the orbit for the eccentric // orbit calculated above.  This appears in Meeus before eq. 23.2. and is // called pi in the text. longitudeOfEarthPerihelion[date] := {    T = meeusT[date]    return (102.93735 + 1.71946 T + 0.00046 T^2) degrees } // 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 precision. // 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 changes in right ascension and declinaton due to nutation.     The calculations above calculate nutation in ecliptical latitude and     longitude, but this recasts those corrections to right ascension and     declination.  This is equation 23.1 in Meeus.     arguments:        date:  The date to be calculated     returns [deltaRA, deltaDecl] , the corrections to add to right ascension                                    and declination respectively. */ highAccuracyNutationInRADecl[date, ra, decl] := {    epsilon = meanObliquityOfEcliptic[date]    [deltapsi, deltaepsilon] = highAccuracyNutation[date]    ce = cos[epsilon]    se = sin[epsilon]    sa = sin[ra]    ca = cos[ra]    td = tan[decl]    deltaRA = (ce + se sa td) deltapsi - (ca td) deltaepsilon    deltaDecl = (se ca) deltapsi + sa deltaepsilon    return [deltaRA, deltaDecl] } /** Calculate the effect of annual aberration.  This is Meeus 23.2 and 23.3     TODO:  Possibly implement the Ron-Vondrak equations for aberration, see     Meeus p. 153     returns:        [deltaRA, deltaDecl] : the corrections to be added to right ascension                               and declination. */ annualAberration[date, ra, decl] := {    eccentricity = earthEccentricity[date]    perihelionLong = longitudeOfEarthPerihelion[date]    k = 20.49552 arcsec  // Constant of aberration    sunLong = sunTrueLongitude[date]    epsilon = meanObliquityOfEcliptic[date]    // Equations 23.3    deltaRA = -k (cos[ra] cos[sunLong] cos[epsilon] + sin[ra] sin[sunLong]) / cos[decl] + eccentricity k (cos[ra] cos[perihelionLong] cos[epsilon] + sin[ra] sin[perihelionLong]) / cos[decl]    deltaDecl = -k (cos[sunLong] cos[epsilon] (tan[epsilon] cos[decl] - sin[ra] sin[decl]) + cos[ra] sin[decl] sin[sunLong]) + eccentricity k (cos[perihelionLong] cos[epsilon] (tan[epsilon] cos[decl] - sin[ra] sin[decl]) + cos[ra] sin[decl] sin[perihelionLong])    return [deltaRA, deltaDecl] } // 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] // //   where: //     lambda:  ecliptic longitude //     beta:    ecliptic latitude //     epsilon: obliquity of the ecliptic // // 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 ecliptical coordinates of an object given right ascension // and declination. // // arguments: //  [ra, decl, epsilon] // //   where: //     ra:   right ascension (can be angle or time) //     decl: declination //     epsilon: obliquity of the ecliptic // // returns: //  [lambda, beta] // //  where: //     lambda:  ecliptic longitude //     beta:    ecliptic latitude raDeclToEcliptical[ra, decl, epsilon] := {    // We want all of this in angles... correct if necessary.    if (ra conforms hour)       ra = ra (circle/day)    lambda = arctan[sin[ra] cos[epsilon] + tan[decl] sin[epsilon], cos[ra]]    beta = arcsin[sin[decl] cos[epsilon] - cos[decl] sin[epsilon] sin[ra]]    return [lambda, beta] } /** Convert from right ascension and declination to galactic coordinates.     This implements Meeus equations 13.7 and 13.8.     Also see:     Blaauw, A., Gum, C. S., Pawsey, J. L., & Westerhout, G. (1960).     The New I.A.U. System of Galactic Coordinates (1958 Revision).     Monthly Notices of the Royal Astronomical Society, 121(2), 123–131.     doi:10.1093/mnras/121.2.123     Available at:     http://sci-hub.tw/10.1093/mnras/121.2.123     arguments:       [ra, decl]     where:        ra:   right ascension (can be angle or time)        decl: declination      both are referred to the the standard equinox of B1950.0.  If a different      mean place (e.g. 2000.0) of the star is given instead of the 1950.0 mean      place, you need to convert ra and decl to ra_1950 and decl_1950.  Use      the "apparentPosition" function below to perform this conversion.      See Meeus 21.2     returns:       [l, b]        where:        l:  galactic longitude        b:  galactic latitude      both are referred to the standard equinox of B1950.0. */ raDeclToGalactic[ra, decl] := {    // We want all of this in angles... correct if necessary.    if (ra conforms hour)       ra = ra (circle/day)    // Galactic north pole right ascension in 1950 coordinates    poleRA = 192.25 degrees    // Galactic north pole declination in 1950 coordinates    poleDecl = 27.4 degrees    x = arctan[sin[poleRA - ra],               cos[poleRA - ra] sin[poleDecl] - tan[decl] cos[poleDecl]] // 13.7    l = (303 degrees - x) mod circle    // eq. 13.8    b = arcsin[sin[decl] sin[poleDecl] + cos[decl] cos[poleDecl] cos[poleRA - ra]]    return [l, b] } /** Convert from right ascension and declination to galactic coordinates.     This implements the equations after Meeus equations 13.7 and 13.8     arguments:       [l, b]     where:        l:  galactic longitude        b:  galactic latitude      both of these should be referenced to the standard equinox of B1950.0.     returns:       [ra, decl]     where:        ra:   right ascension        decl: declination      both are referred to the standard equinox of B1950.0.    If a different      mean place (e.g. 2000.0) of the star is desired instead of the 1950.0 mean      place, you need to convert ra and decl using       the "apparentPosition" function below to perform this conversion.      See Meeus 21.2 */ galacticToRADecl[l, b] := {    // Galactic north pole right ascension in 1950 coordinates    poleRA = 192.25 degrees    // Galactic north pole declination in 1950 coordinates    poleDecl = 27.4 degrees    y = arctan[sin[l - 123 degrees],               cos[l - 123 degrees] sin[poleDecl] - tan[b] cos[poleDecl]]    alpha = (y + 12.25 degrees) mod circle    delta = arcsin[sin[b] sin[poleDecl] + cos[b] cos[poleDecl] cos[l - 123 degrees]]    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 (angle) // 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 if 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, (which is very uncharacteristic for Meeus' // extraordinarily good book,) 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) //      if (abs[correction] > 6 hours) //         println["Correction is \$correction"]              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"] //      println["Altitude 1: " + (altitude1 -> "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] = 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 chapters 25 and 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. // Returns undef if the sun does not rise or set on that date at that latitude. sunrise[date, lat, long, temp = 283 K, pressure=1010 millibars ] := {    edate = approxMidnight[date,long] + 5.5 hours //   println["edate is " + (edate->Guam)]    [ra, decl] = sunApparentRADecl[edate]    [rise, set] = approxRiseSet[edate, lat, long, ra, decl] //   println["arise is " + (rise->Guam)]    if rise == undef       return undef    else       return sunSecantAltitude[rise, lat, long, -16 arcmin, temp, pressure] } // Calculate the time when the upper limb of the sun just crosses the // horizon. // Returns undef if the sun does not rise or set on that date at that latitude. sunset[date, lat, long, temp = 283 K, pressure=1010 millibars ] := {    edate = approxMidnight[date,long] + 17.5 hours    [ra, decl] = sunApparentRADecl[edate]    [rise, set] = approxRiseSet[edate, lat, long, ra, decl]    if set == undef       return undef    else       return sunSecantAltitude[set, lat, long, -16 arcmin, temp, pressure] } // Calculate the time of the sun's transit (that is, its highest point in // the sky. sunTransit[date, lat, long, temp = 283 K, pressure=1010 millibars ] := {    approxNoon = approxMidnight[date, long] + 12 hours    [az, alt] = refractedSunAzimuthAltitude[approxNoon, lat, long]    // This is unstable when the sun is directly overhead, so we'll    // correct a bit (getting the same answer.)    if (alt > 85 degrees)    { //      println["azimuth is " + (az->"degrees")] //      println["altitude is " + (alt->"degrees")] //      println["Correcting."]       if lat > 0 degrees          lat = lat + 5 degrees       else          lat = lat - 5 degrees       [az, alt] = refractedSunAzimuthAltitude[approxNoon, lat, long]    }     //   println["azimuth is " + (az->"degrees")] //   println["altitude is " + (alt->"degrees")]    if abs[az] < 90 degrees    // North of sun, looking south.       targetAz = 0 degrees    else       targetAz = 180 degrees        return sunSecantAzimuth[approxNoon, lat, long, targetAz, 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] // Calculate the darkness of the sky at a specific time. // The (optional) symbols passed in are an array, containing values that are // returned for the sun's azimuth when it is: // //  0: astronomical twilight ( lower than -18 degrees) //  1: civil twilight        ( lower than -12 degrees) //  2: nautical twilight     ( lower than  -6 degrees) //  3: sunset                ( lower than   0 degrees) //  4: daylight              ( above 0 degrees)  // //  respectively skyDarkness[date, lat, long, syms = ["****", "***", "**", "*", ""] ] := {    [az, alt] = airlessSunAzimuthAltitude[date, lat, long]    if alt < -18 degrees       return syms@0    if alt < -12 degrees       return syms@1    if alt < -6 degrees       return syms@2    if alt < 0 degrees       return syms@3    else       return syms@4 } // Get a time object representing the approximate midnight that begins // 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 (numerically) the moon's radial velocity towards earth. moonRadialVelocity[date] := {    halfStep = .5 s    d1 = date - halfStep    d2 = date + halfStep    return (moonDistance[d1] - moonDistance[d2]) / (halfStep * 2) } // 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 } // The position angle of the moon's axis of rotation (this is *not* the // illuminated axis,) reckoned eastward from the zenith. // 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. moonAxisAngleRelativeToZenith[date, lat, long, heightAboveSeaLevel] := {    axis = moonAxisAngle[date, lat, long, heightAboveSeaLevel]    [ra, decl] = moonApparentRADecl[date]    q = parallacticAngle[date, lat, long, ra, decl]    return axis - q } // The position angle of the moon's axis of rotation (this is *not* the // illuminated axis,) 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.) // See chapter 53 in Meeus. and the moonLibrationAndAxis formula. moonAxisAngle[date, lat, long, heightAboveSeaLevel] := {    // TODO:  pass in lambda and beta for lat, long    // See Meeus chapter 53, "Topocentric Librations" section    [ra, decl] = moonApparentRADecl[date]    distance = moonDistance[date]    [ra1, decl1] = correctRADeclForParallax[date, ra, decl, lat, long, distance, heightAboveSeaLevel]    epsilon = trueObliquityOfEcliptic[date]    [lambda, beta] = raDeclToEcliptical[ra1, decl1, epsilon]    [l, b, P] = moonLibrationAndAxis[date, lambda, beta]        return P }            // Returns a GeneralPath that represents the sunlit part of a body.  This can // be either the moon or a planet. //   arguments: //     cx, cy: centerpoint of body //     radius: radius of body //     angle:  an angle, perhaps as obtained by one of the moonPositionAngle... //             functions above.  This is the angle of the center of the bright //             limb. //     illum:  the illuminated fraction of the body //     filled: A boolean value indicating if the polygon is filled or not. drawSunlitPolygon[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 } // 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.  This is the angle of the center of the bright //             limb. //     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] := {    drawSunlitPolygon[cx, cy, radius, angle, illumFraction, filled] } // Returns a GeneralPath that represents the position of the moon relative // to the zenith. // This should be added to a graphics object with the graphics.add[x] method // 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 fill polygon. drawMoonPolygonRelativeToZenith[date, lat, long, cx, cy, radius, filled] := {    return drawMoonPolygon[cx, cy, radius, moonPositionAngleRelativeToZenith[date, lat, long], moonIlluminatedFraction[date], filled] } // Returns a graphics object that represents the axis of the moon or a planet. //   arguments: //     cx, cy: centerpoint of object //     radius: radius of object (polar radius) //     angle:  an angle, as obtained by one of the AxisAngle... //             functions above.  The angle is counterclockwise from north. //     axisradius:  The length of the axis to draw as a multiple of the //                  object's radius drawRotationAxis[cx, cy, radius, angle, axisradius=1.1] := {    g = new graphics    len = radius axisradius    sa = sin[angle]    ca = cos[angle]        x1 = cx - len sa    x2 = cx - radius sa    y1 = cy - len ca    y2 = cy - radius ca    g.line[x1,y1,x2,y2]    x1 = cx + len sa    x2 = cx + radius sa    y1 = cy + len ca    y2 = cy + radius ca    g.line[x1,y1,x2,y2]    return g } // Returns a graphics object that represents the axis of the moon. //   arguments: //     cx, cy: centerpoint of moon //     radius: radius of moon //     angle:  an angle, as obtained by one of the moonAxisAngle... //             functions above. //     axisradius:  The length of the axis to draw as a multiple of the //                  moon's radius drawMoonRotationAxisRelativeToZenith[date, lat, long, heightAboveSeaLevel, cx, cy, moonradius, axisradius=1.1] := {    angle = moonAxisAngleRelativeToZenith[date, lat, long, heightAboveSeaLevel]    return drawRotationAxis[cx, cy, moonradius, angle, axisradius] } // 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 } // Parallax angle corrections (in right ascension and declination) 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. // RETURNS: //  [ra, decl]  corrected for parallax correctRADeclForParallax[date, ra, decl, lat, long, distance, heightAboveSeaLevel] := {    // This is more clear than eq. 40.1    azParallax = arcsin[earthradius_equatorial / distance]    [rsp, rcp] = rhoSinPhi[lat, heightAboveSeaLevel]        H = hourAngle[date, long, ra]    sp = sin[azParallax]    deltaRA = arctan[-rcp sp sin[H], cos[decl] - rcp sp cos[H]] // Eq 40.2    raPrime = ra + deltaRA    declPrime = arctan[(sin[decl] - rsp sp) cos[deltaRA],  // Eq 40.3                       cos[decl] - rcp sp cos[H]]    return [raPrime, declPrime] } /** Calculate the values of rho sin phi' and rho cos phi' used in parallax     calculations.  These are defined in chapter 11 of Meeus, and used in     parallax calculations (chapter 40.)     Arguments:        [lat, heightAboveSeaLevel]       where          lat: the ordinary geographical latitude (phi)     returns:      [rho sin phi', rho cos phi'] */ rhoSinPhi[lat, heightAboveSeaLevel] := {    a = earthradius_equatorial    f = earth_flattening    b = a(1-f)   // Earth polar radius    bOvera = b/a    u = arctan[bOvera tan[lat]]    rhoSinPhi = bOvera sin[u] + heightAboveSeaLevel/a sin[lat]    rhoCosPhi = cos[u] + heightAboveSeaLevel/a cos[lat]    return [rhoSinPhi, rhoCosPhi] } // Returns the angular separation between two bodies given their right // ascensions and declinations [ra1, decl1, ra2, decl2] //   Note that when using right ascension and declination, these are geocentric // coordinates and do not take parallax into account, so they are useless for // eclipses or transit predictions or close alignments with the moon and other // objects. // // It also works for altazimuth coordinates if passed: [az1, alt1, az2, alt2] // In the altazimuth case, you can correct for parallax and refraction from // a particular point on the earth, and if you did so, you may use this to // predict eclipses, transits, etc. successfully. angularSeparation[ra1, decl1, ra2, decl2] := { /*   println["\$ra1 \$decl1 \$ra2 \$decl2"]    println["sin[decl1] = " + sin[decl1]]    println["sin[decl2] = " + sin[decl2]]    println["cos[decl1] = " + cos[decl1]]    println["cos[decl2] = " + cos[decl2]]    println["cos[ra1-ra2] = " + cos[ra1-ra2]]    println["sin[decl1] sin[decl2] + cos[decl1] cos[decl2] cos[ra1-ra2] = " + sin[decl1] sin[decl2] + cos[decl1] cos[decl2] cos[ra1-ra2]]*/        arg = sin[decl1] sin[decl2] + cos[decl1] cos[decl2] cos[ra1-ra2]    // We want to condition the argument to allow intervals, as typical    // intervals wildly overestimate to the point that the values aren't within    // the domain of arccos.  We'll constrain arguments to the domain of arccos,    // [-1,1]    constrainedArg = intersection[arg, new interval[-1,1]]    return arccos[constrainedArg] } // Convenience function to obtain the angular separation between the sun // and moon for a given time. // WARNING:  This is the GEOCENTRIC angular separation, that is, as seen from // the center of the earth and is not corrected for refraction or parallax, // so it is useless for, say, predicting eclipses or transits! // In other words, this is probably not the function you're looking for. // You most likely want to use the function below that requires latitude and // longitude. geocentricSunMoonAngularSeparation[date] := {    [sunRA, sunDecl] = sunApparentRADecl[date]    [moonRA, moonDecl] = moonApparentRADecl[date]    return angularSeparation[sunRA, sunDecl, moonRA, moonDecl] } // This returns the angular separation of the sun and moon as seen from // a particular point on the earth.  It is corrected for parallax. sunMoonAngularSeparation[date, lat, long] := {    [sunaz, sunalt]   =  refractedSunAzimuthAltitude[date, lat, long]    [moonaz, moonalt] = refractedMoonAzimuthAltitude[date, lat, long]    return angularSeparation[sunaz, sunalt, moonaz, moonalt] } // Calculates the approximate transit time of a body which has the specified // right ascension and declination.  The resultant value will be close // to the specified date.  See Section 15.2 of Meeus. // This value should not be used directly, but should be used as a seed in // something like moonTransit[] approxTransit[date, long, ra, decl] := {    m0 = (ra - apparentLocalSiderealAngle[date, long]) / (360 degrees / day)    m0 = m0 mod (24 hours) //   println["m0 is \$m0"]    if m0 > 12 hours       m0 = m0 - 24 hours    if m0 < -12 hours       m0 = m0 + 24 hours    return date + m0 } // Calculate the approximate offset of rise/set for a body with the specified // right ascension and declination as seen from the specified latitude. // h0 is the "standard altitude" with refraction.  The default value is // accurate for point objects, or -0.8333 degrees should be used for the sun. // //  Returns hour angle or undef if the object is circumpolar. //  THINK ABOUT:  Throw an exception? calcHourAngle[lat, decl, h0 = -0.5667 degrees] := {    cosH0 = (sin[h0] - sin[lat] sin[decl]) / (cos[lat] cos[decl])    if abs[cosH0] > 1       return undef    else       return arccos[cosH0] mod (180 degrees) } // Calculates the approximate rise and set time of a body which has the // specified right ascension and declination.  The resultant value will be // close to the specified date.  See Equations 15.1 - 15.2 of Meeus. // This value should not be used directly, but should be used as a seed in // something like moonTransit[] // //  returns: //    [rise, set] //  Return values are undef if the object is circumpolar for that date. approxRiseSet[date, lat, long, ra, decl, h0 = -0.5667 degrees] := {    transit = approxTransit[date, long, ra, decl]    H0 = calcHourAngle[lat, decl, h0]        if H0 == undef       return [undef, undef]        Htime = H0 / (360 degrees/day) //   println["transit is \$transit"] //   println["H0 is \$H0"] //   println["Htime is " + (Htime -> "hours")]    rise = transit - Htime    set = transit + Htime    return [rise,set] } // Calculate the time when the upper limb of the moon just crosses the // horizon. // Returns undef if the sun does not rise or set on that date at that latitude. moonrise[date, lat, long, temp = 283 K, pressure=1010 millibars ] := {    [ra, decl] = moonApparentRADecl[date]    [rise, set] = approxRiseSet[date, lat, long, ra, decl] //   println["Rise guess is \$rise"]    if rise == undef       return undef    else       return moonSecantAltitude[rise, lat, long, -16 arcmin, temp, pressure] } // Calculate the time when the upper limb of the moon just crosses the // horizon. // Returns undef if the sun does not rise or set on that date at that latitude. moonset[date, lat, long, temp = 283 K, pressure=1010 millibars ] := {    [ra, decl] = moonApparentRADecl[date]    [rise, set] = approxRiseSet[date, lat, long, ra, decl] //   println["Set guess is \$set"]    if set == undef       return undef    else       return moonSecantAltitude[set, lat, long, -16 arcmin, temp, pressure] } // Calculate the time of the moon's transit (that is, its highest point in // the sky. moonTransit[date, lat, long, temp = 283 K, pressure=1010 millibars ] := {    [moonRA, moonDecl] = moonApparentRADecl[date]    approxTransit = approxTransit[date, long, moonRA, moonDecl]    //   println["Approximate transit is \$approxTransit"]    [az, alt] = refractedMoonAzimuthAltitude[approxTransit, lat, long]    // This is unstable when the moon is directly overhead, so we'll    // correct a bit (getting the same answer.)    if (alt > 85 degrees)    { //      println["azimuth is " + (az->"degrees")] //      println["altitude is " + (alt->"degrees")] //      println["Correcting."]       if lat > 0 degrees          lat = lat + 5 degrees       else          lat = lat - 5 degrees       [az, alt] = refractedMoonAzimuthAltitude[approxTransit, lat, long]    }     //   println["azimuth is " + (az->"degrees")] //   println["altitude is " + (alt->"degrees")]    if abs[az] < 90 degrees    // North of moon, looking south.       targetAz = 0 degrees    else       targetAz = 180 degrees        return moonSecantAzimuth[approxTransit, lat, long, targetAz, 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)] // Calculate the time of perigee of the moon closest to the specified date. //  Returns:  The date of perigee. moonPerigee[date] := {    return moonMeanApogeeOrPerigeeInternal[date, 0] } // Calculate the time of apogee of the moon closest to the specified date. //  Returns:  The date of apogee. moonApogee[date] := {    return moonMeanApogeeOrPerigeeInternal[date, 0.5] } // Estimate the apogee or perigee of the moon closest to a particular // date.  See Meeus chapter 50.  This should not be called directly by the // user.  Its arguments should be preconditioned by calling // moonApogee[date] or moonPerigee[date] moonMeanApogeeOrPerigeeInternal[date, correctionFactor] := {    epoch = #1999-12-22 UTC#   // Definition of k=0 in Meeus    // Round to closest integer or half-integer.    k = round[(((date - epoch) / year) * 13.2555) + correctionFactor, 1] - correctionFactor    //   println["k is \$k"]        T = k/1325.55    // Meeus 50.3    // println["T is \$T"]    jde = JDE[2451534.6698 + 27.55454989 k - 0.0006691 T^2 - 0.000001098 T^3 + 0.0000000052 T^4]    // println["jde is \$jde " + (jde->JDE)]    // We now use the secant method to find when the moon's radial velocity    // relative to earth becomes zero around this point.  This doesn't use the    // rest of the procedure found in Meeus ch. 50, but rather allows us to    // later plug in a more accurate moon position (such as the *FULL* ELP-2000    // model) and this will automatically get more accurate.    d1 = jde    d2 = jde + .1 s    v1 = moonRadialVelocity[d1]    v2 = moonRadialVelocity[d2]    do    {       correction = (v1 * (d1 - d2)) / (v1 - v2)       dnew = d1 - correction       //   println["Correction is \$correction"]       if abs[correction] < 1 second          return d1       d2 = d1       d1 = dnew       v2 = v1       v1 = moonRadialVelocity[d1]    } while true } /** Calculate the longitude of the mean ascending node (Omega).     This is equation 47.7 in Meeus. */ moonMeanLongitudeOfAscendingNode[date] := {    T = meeusT[date]    return (125.0445479 - 1934.1362891 T + 0.0020754 T^2 +            T^3 / 467441 - T^4 / 60616000) degrees } /** Calculate the optical libration of the moon.   This is chapter 53 in     Meeus.  It also calculates the position angle of the moon's axis because     it requires all of these calculations     The parameters lambda and beta should be left as undef if you want the     geocentric libration value (That is, as seen from the center of the     earth.)     They are allowed to be passed in if you want the libration for     an exact position on the earth.  See "Topocentric Librations" in chapter     53 to see how to do this.  In this case, lambda should be the apparent     geocentric longitude and beta the apparent geocentric latitude of the moon.     returns: [l, b, P]       l:  physical libration in longitude       b:  physical libration in latitude       P:  the position angle of the Moon's axis of rotation  (See the diagram           in Chapter 42 of Meeus) */ moonLibrationAndAxis[date, lambda=undef, beta=undef] := {    // Inclination of the mean lunar equator to the ecliptic    I = 1.54242 degrees    var lambdaMinusNotation    [deltaPsi, deltaEpsilon] = highAccuracyNutation[date]        if (lambda==undef)    {       // This is the moon longitude *not corrected for nutation* because       // in this case we're going to just subtract out the nutation.       lambdaMinusNutation = moonLongitude[date]    } else    {       // We expect the apparent longitude to have been passed in (that is,       // longitude corrected for nutation) so we subtract out the nutation.       lambdaMinusNutation = lambda - deltaPsi    }    if beta == undef       beta = moonLatitude[date]    // First, calculate the optical librations.  These are the larger source    // of libration for the moon, ranging up to several degrees.    F = moonArgumentOfLatitude[date]    Omega = moonMeanLongitudeOfAscendingNode[date]    // The following are equations 53.1 in Meeus    W = lambdaMinusNutation - Omega    cb = cos[beta]    sb = sin[beta]    sw = sin[W]    ci = cos[I]    si = sin[I]    A = arctan[sw cb cos[I] - sb si, cos[W] cb]    lprime = A - F    // Normalize to small angles    if lprime > 180 degrees       lprime = lprime - 360 degrees    else       if lprime < -180 degrees          lprime = lprime + 360 degrees           bprime = arcsin[-sw cb si - sb ci]    T = meeusT[date]    K1 = (119.75 + 131.849 T) degrees    K2 =  (72.56 +  20.186 T) degrees    D = moonMeanElongation[date]    M = moonCalcSunMeanAnomaly[date]    Mprime = moonMeanAnomaly[date]    E = moonCalcEarthEccentricity[date]    rho = (-0.02752 cos[Mprime] - 0.02245 sin[F] + 0.00684 cos[Mprime - 2 F] +       -0.00293 cos[2 F] - 0.00085 cos[2 F - 2 D] - 0.00054 cos[Mprime - 2 D] +       -0.00020 sin[Mprime + F] - 0.00020 cos[Mprime + 2 F] +       -0.00020 cos[Mprime - F] + 0.00014 cos[Mprime + 2 F - 2 D]) degree    sigma = (-0.02816 sin[Mprime] + 0.02244 cos[F] - 0.00682 sin[Mprime - 2 F] +       -0.00279 sin[2 F] - 0.00083 sin[2 F - 2 D] + 0.00069 sin[Mprime - 2 D] +        0.00040 cos[Mprime + F] - 0.00025 sin[2 Mprime] +       -0.00023 sin[Mprime + 2 F] + 0.00020 cos[Mprime - F] +        0.00019 sin[Mprime - F] + 0.00013 sin[Mprime + 2 F - 2 D] +       -0.00010 cos[Mprime - 3 F]) degree    tau = (0.02520 E sin[M] + 0.00473 sin[2 Mprime - 2 F] - 0.00467 sin[Mprime]+       0.00396 sin[K1] + 0.00276 sin[2 Mprime - 2 D] + 0.00196 sin[Omega] +      -0.00183 cos[Mprime - F] + 0.00115 sin[Mprime - 2 D] +      -0.00096 sin[Mprime - D] + 0.00046 sin[2 F - 2 D] +      -0.00039 sin[Mprime - F] - 0.00032 sin[Mprime - M - D] +       0.00027 sin[2 Mprime - M - 2 D] + 0.00023 sin[K2] +      -0.00014 sin[2 D] + 0.00014 cos[2 Mprime - 2 F] +      -0.00012 sin[Mprime -  2 F] - 0.00012 sin[2 Mprime] +       0.00011 sin[2 Mprime - 2 M - 2 D]) degree    l2prime = -tau + (rho cos[A] + sigma sin[A]) tan[bprime]    b2prime = sigma cos[A] - rho sin[A]    // println["lprime=" + (lprime->"deg") +  ", l2prime=" + (l2prime->"deg") + ", bprime=" + (bprime->"deg") + ", b2prime=" (b2prime->deg)]    l = lprime + l2prime    b = bprime + b2prime           // Now calculate the position angle P of the moon's axis (chapter 53)    V = Omega + deltaPsi + sigma/si    // println["deltaPsi=" + (deltaPsi->"deg")]    // println["V=" + (V->"deg")]    // println["rho=" + (rho->"deg")]    sir = sin[I + rho]    X = sir sin[V]           epsilon = trueObliquityOfEcliptic[date]    // println["epsilon is " + (epsilon->"deg")]    Y = sir cos[V] cos[epsilon] - cos[I + rho] sin[epsilon]    omega = arctan[X, Y]           [ra, decl] = moonApparentRADecl[date]    P = arcsin[sqrt[X^2 + Y^2] cos[ra - omega] / cos[b]] mod circle    return [l, b, P] } /** Calculate precession corrections for the specified time.     These are the "rigorous method" calculations for precession from Meeus,     chapter 21.  The following is equation 21.2.     arguments:     targetEpoch:  The target time to calculate the precession.     initialEpoch: The initial epoch in which coordinates are specified.     Since star positions are usually referenced to the epoch J2000.0, this     is given as the default value for initialEpoch.     returns [zeta, z, theta] */ precessionCorrections[targetEpoch, initialEpoch = # 2000-01-01 12:00 TD #] := {    JD0 = JD[initialEpoch] / days    JD  = JD[targetEpoch]  / days    T = (JD0 - 2451545.0) / 36525    t = (JD - JD0) / 36525        zeta = ((2306.2181 + 1.39656 T - 0.000139 T^2) t +            (0.30188 - 0.000344 T) t^2 + 0.017998 t^3) arcsec    z =    ((2306.2181 + 1.39656 T - 0.000139 T^2) t +            (1.09468 + 0.000066 T) t^2 + 0.018203 t^3) arcsec    theta = ((2004.3109 - 0.85330 T - 0.000217 T^2) t -              (0.42665 + 0.000217 T) t^2 - 0.041833 t^3) arcsec    return [zeta, z, theta] } /** Precesses the coordinates to the target epoch and to the equinox     of the target date.  These implement equations 21.4 in Meeus.     arguments:       date:  The target date       ra0:   The right ascension referred to the given refEpoch, as an angle              or a time.       decl0: The declination referred to the given refEpoch       refEpoch:  The epoch of the original coordinates.          Since star positions are usually referenced to the epoch J2000.0, this          is given as the default value for refEpoch.     returns:     [ra, decl] : The precessed right ascension and declination of the                  coordinates related to the equinox of date. */ precess[date, ra0, decl0, refEpoch = # 2000-01-01 12:00 TD #] := {    // Correct ra to angle if necessary    if ra0 conforms time       ra0 = ra0 * (circle/day)        [zeta, z, theta] = precessionCorrections[date, refEpoch] //   println["\nPrecession corrections:"] //   println["zeta:  " + (zeta -> "degrees")] //   println["z: " + (z -> "degrees")] //   println["theta: " + (theta -> "degrees")]    // These are equations 21.4    cd0 = cos[decl0]    sd0 = sin[decl0]    aoc = ra0 + zeta    saoc = sin[aoc]    caoc = cos[aoc]    ctheta = cos[theta]    stheta = sin[theta]    A = cd0 saoc    B = ctheta cd0 caoc - stheta sd0    C = stheta cd0 caoc + ctheta sd0    raMinusZ = arctan[A, B]    ra = raMinusZ + z    decl = arcsin[C]    return [ra, decl] } /** Calculates the apparent position of a star to high accuracy using     coordinates from a specified equinox.  It corrects for:       * Proper motion (assumed to be uniform)       * Precession       * Nutation       * Annual aberration     It does *NOT* correct for refraction, and these calculations are     geocentric, that is, referred to the center of the earth.     You will want to use something like refractedAzimuthAltitude to find its     refracted azimuth and altitude for a specific spot on the earth.     See chapter 23 in Meeus for more explanation.    arguments:     date:  The target date for calculation     ra0:   The right ascension in some reference epoch (as angle or time)     decl0: The declination in some reference epoch (as an angle)     raV:   The proper motion in right ascension (as angle or time) per time     declV: The proper motion in declination (as angle or time) per time     distance:  The distance to the star in some reference epoch (or undef if                unknown)     radialV:  The radial velocity of the star relative to the solar system                (positive means distance increasing) or undef if unknown     refEpoch: The reference epoch in which the coordinates are specified.       Since star positions are often referenced to the epoch J2000.0, this       is given as the default value for initialEpoch.     distance and radialV are rarely important.  Even a close star like Sirius     with high radial velocity gives an error of less than 0.04 arcsecond     between 1900 and 2100, but may be necessary when extrapolating thousands     of years into the past or future.     returns:     [ra, decl] in the epoch and mean equinox of the date */ apparentPosition[date, ra0, decl0, raV, declV, distance = undef, radialV = undef, refEpoch = # 2000-01-01 12:00 TD #] := {    // Correct ra to angle if necessary    if ra0 conforms time       ra0 = ra0 * (circle/day)    // Correct velocity in right ascension to angle/time    if raV conforms (time/time)       raV = raV * (circle/day)    if (distance == undef OR radialV == undef)    {       // Correct for proper motion without radial velocity       timeOffset = date - refEpoch       dra = raV * timeOffset       ddecl = declV * timeOffset       /*       println["\nProper motion correction:"]       println["ra  : " + (dra/(circle/day)-> "s")]       println["decl: " + (ddecl-> "arcsec")]       */       ra0 = ra0 + dra       decl0 = decl0 + ddecl    } else    {       // Correct for proper motion with known radial velocity and distance       // Meeus p. 141       // Since Frink tracks the units of measure, you don't have to go through       // the convolutions of turning things into parsecs/year, radians/year,       // etc.        x = distance cos[decl0] cos[ra0]       y = distance cos[decl0] sin[ra0]       z = distance sin[decl0]       dx = (x/distance) radialV - z declV cos[ra0] - y raV       dy = (y/distance) radialV - z declV sin[ra0] + x raV       dz = (z/distance) radialV + distance declV cos[decl0]       timeOffset = date - refEpoch       xp = x + timeOffset dx       yp = y + timeOffset dy       zp = z + timeOffset dz       ra0 = arctan[yp, xp]       decl0 = arctan[zp, sqrt[xp^2 + yp^2]]    }    /*    println["\nAfter correcting for proper motion:"]    println["ra  : " + (ra0-> "degrees")]    println["decl: " + (decl0-> "degrees")]    */    // Now correct for precession.  [ra0, decl0] are now the coordinates for the    // for the mean equinox of the reference time, but for the target epoch.    [ra, decl] = precess[date, ra0, decl0, refEpoch]    /*    println["\nAfter correcting for precession:"]    println["ra  : " + (ra-> "degrees") + "\t" + (ra / (circle/day) -> HMS)]    println["decl: " + (decl-> "degrees") + "\t" + (decl -> DMS)]    */    // Now correct for nutation.  (See Meeus p. 150).    [deltaRA1, deltaDecl1] = highAccuracyNutationInRADecl[date, ra, decl]    /*    println["\nCorrections for nutation:"]    println["ra  : " + (deltaRA1-> "arcsec")]    println["decl: " + (deltaDecl1 -> "arcsec")]    */    // Now correct for annual aberration.  See Meeus eq. 23.2    [deltaRA2, deltaDecl2] = annualAberration[date, ra, decl]    /*    println["\nCorrections for aberration:"]    println["ra  : " + (deltaRA2-> "arcsec")]    println["decl: " + (deltaDecl2 -> "arcsec")]    */        ra = ra + deltaRA1 + deltaRA2    decl = decl + deltaDecl1 + deltaDecl2    return[ra, decl] } "sun.frink included Ok" ```

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 18618 days, 15 hours, 41 minutes ago.