vsop87.frink

View or download vsop87.frink in plain text format

//
// Routines for parsing VSOP87 coefficient files and generating Frink code
// from them (or, with small modifications, for other languages).
// The output of running this program is a partial Frink program
// that contains function definitions to find the coordinates of each planet.
//
// The full coefficients of the VSOP87 theory are available for download
// from:
//  ftp://ftp.imcce.fr/pub/ephem/planets/vsop87/
//
// From this, the series we require for use with the equations in Meeus
// are the VSOP87D.xxx files, which contains the
// heliocentric spherical variables referred to equinox and ecliptic of date.
//
// This reverses the lines in the file so that smaller coefficients are
// added first, reducing numerical error.

planets = ["Mercury", "Venus", "Earth", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune"]

for planet = planets
{
   ext = lc[left[planet, 3]]

   firstline = true
   fullvar = undef
   running = undef

   println["""
//
// $planet
//   
// This function calculates the heliocentric coordinates of $planet
// referred to the mean equinox *of the date*.  You may want to convert this
// to another coordinate system, such as FK5.
//
// arguments:
//     d: the date/time to be calculated for
//   
// returns:
//   [L, B, R]
//   
// Where 
//       L is the heliocentric longitude,
//       B is the heliocentric latitude
//       R is the distance from the sun.
${planet}HeliocentricCoordinates = {|d|

   tau = meeusT[d] / 10
   """]

   for line=lines["file:vsop87/VSOP87D.$ext"]
   {
      if [planet, varno, varnames, exponent] = line =~ %r/VSOP87\s+VERSION\s+D4\s+(\w+)\s+VARIABLE\s+(\d)\s+\((\w+)\)\s+\*T\*\*(\d)/
      {
         varno = parseInt[varno]
         varname = substrLen[varnames, varno-1, 1]
         firstline = true
         if (fullvar)
            println["   $fullvar = $running;\n"]

         fullvar = varname+exponent
         running = ""
         //      println["Planet: $planet\tvarno: $varno\tvarname: $varname\texponent: $exponent\tname: $fullvar"]
      } else
      {
         A = eval[substr[line, 79, 97]]
         B = eval[substr[line, 97, 111]]
         C = eval[substr[line, 111, 132]]

         if length[line] != 132
            println[length[line]]

         // Remove lines with coefficients of 0.0
         if (eval[A] != 0.0)
         {
            running = "$A * cos[$B + $C * tau]" + (firstline ? "" : " +\n        ") + running
            firstline = false
         }
      }
   }

   println["""
   $fullvar = $running;

   L = ((L0 + L1 tau + L2 tau^2 + L3 tau^3 + L4 tau^4 + L5 tau^5) radians) mod circle
   B = ((B0 + B1 tau + B2 tau^2 + B3 tau^3 + B4 tau^4 + B5 tau^5) radians) mod circle
   R = (R0 + R1 tau + R2 tau^2 + R3 tau^3 + R4 tau^4 + R5 tau^5) au

   return [L, B, R]
}"""]
}


View or download vsop87.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 17651 days, 5 hours, 15 minutes ago.