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 17592 days, 15 hours, 38 minutes ago.