sun.frink

Download or view sun.frink in plain text format


// 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 (equatorial coordinates)
// 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 (equatorial coordinates).
//
// 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

    Also see:
    https://gea.esac.esa.int/archive/documentation/GDR2/Data_processing/chap_cu3ast/sec_cu3ast_intro/ssec_cu3ast_intro_tansforms.html#SSS1

    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 epoch.  It corrects for:
      * Proper motion (assumed to be uniform)
      * Radial velocity
      * 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"


Download or view 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 19935 days, 10 hours, 6 minutes ago.