#!/usr/bin/perl

use strict;
use warnings;

my @vsop_planets = qw( mer ven ear mar jup sat ura nep ); # not emb

my @smart_variables = qw( psi omega phi );

my %vsop_accuracy = (
  mer => 1E-3,
  ven => 1E-3,
  ear => 1E-6,
  mar => 1E-3,
  jup => 1E-3,
  sat => 1E-2,
  ura => 1E-2,
  nep => 1E-2,
  "*" => 1E-3,
);

my $elp_accuracy_angle = 1;
my $elp_accuracy_radius = 1;

my $smart_accuracy = 1000000;

my $time_span = 100;

my %vsop_data;

sub getnum($$) {
  my ($s, $l) = @_;
  return substr($_, $s, $l);
}

for my $planet (@vsop_planets) {
  my $fn = "vsop87/VSOP87A.$planet";
  open my $f, "<", $fn or die "$fn: $!\n";
  my $nt = 0;
  my $coord;
  my $deg;
  my $errf;
  my $errt = $vsop_accuracy{$planet} || $vsop_accuracy{"*"};
  while(<$f>) {
    if($nt == 0) {
      die unless /^ VSOP87/;
      $coord = getnum(41, 1) - 1;
      $deg = getnum(59, 1);
      $nt = getnum(60, 7);
      $errf = ($time_span / 1000) ** $deg;
      if(!defined $vsop_data{$planet}[$coord]) {
	$vsop_data{$planet}[$coord] = {
	  error => 0,
	  terms => []
	}
      }
    } else {
      my $a = getnum(79, 18);
      my $b = getnum(97, 14);
      my $c = getnum(111, 20);
      my $errv = abs($a) * $errf;
      if($errv >= $errt) {
	push @{$vsop_data{$planet}[$coord]{terms}[$deg]}, [ $a, $b, $c ];
      } else {
	$vsop_data{$planet}[$coord]{error} += $errv;
      }
      $nt--;
    }
  }
}

sub PI() { 4 * atan2(1, 1); }

sub dms($$$) {
  my ($d, $m, $s) = @_;
  return ($d + $m / 60 + $s / 3600) * PI / 180;
}

# page 11 of elp82b.ps
my $elp_param_d = [
  + dms(297, 51, 0.73512),
  + dms(0, 0, 1602961601.4603),
  - dms(0, 0, 5.8681),
  + dms(0, 0, 0.006595),
  - dms(0, 0, 0.00003184),
];
my $elp_param_lp = [
  + dms(357, 31, 44.79306),
  + dms(0, 0, 129596581.0474),
  - dms(0, 0, 0.5529),
  + dms(0, 0, 0.000147),
];
my $elp_param_l = [
  + dms(134, 57, 48.28096),
  + dms(0, 0, 1717915923.4728),
  + dms(0, 0, 32.3893),
  + dms(0, 0, 0.051651),
  - dms(0, 0, 0.00024470),
];
my $elp_param_f = [
  + dms(93, 16, 19.55755),
  + dms(0, 0, 1739527263.0983),
  - dms(0, 0, 12.2505),
  - dms(0, 0, 0.001021),
  + dms(0, 0, 0.00000417),
];
my $elp_param_w1 = [
  + dms(218, 18, 59.95571),
  + dms(0, 0, 1732559443.73604),
  - dms(0, 0, 5.8883),
  + dms(0, 0, 0.006604),
  - dms(0, 0, 0.00003169),
];
my $elp_param_plt = [
  + dms(100, 27, 59.22059),
  + dms(0, 0, 129597742.2758),
  - dms(0, 0, 0.0202),
  + dms(0, 0, 0.000009),
  - dms(0, 0, 0.00000015),
];
# page 7, bottom, page 9 for constant p
my $elp_param_zeta = [ @$elp_param_w1 ];
$elp_param_zeta->[1] += dms(0, 0, 5029.0966);
# page 8
my $elp_param_plme = [ dms(252, 15, 3.25986), dms(0, 0, 538101628.68898) ];
my $elp_param_plv = [ dms(181, 58, 47.28305), dms(0, 0, 210664136.43355) ];
my $elp_param_plma = [ dms(355, 25, 59.78866), dms(0, 0, 68905077.59284) ];
my $elp_param_plj = [ dms(34, 21, 5.34212), dms(0, 0, 10925660.42861) ];
my $elp_param_pls = [ dms(50, 4, 38.89694), dms(0, 0, 4399609.65932) ];
my $elp_param_plu = [ dms(314, 3, 18.01841), dms(0, 0, 1542481.19393) ];
my $elp_param_pln = [ dms(304, 20, 55.19575), dms(0, 0, 789550.32074) ];

my %elp_params = (
  main => [ $elp_param_d, $elp_param_lp, $elp_param_l, $elp_param_f  ],
  pert => [ $elp_param_zeta, $elp_param_d, $elp_param_lp, $elp_param_l,
    $elp_param_f  ],
  pla1 => [ $elp_param_plme, $elp_param_plv, $elp_param_plt, $elp_param_plma,
    $elp_param_plj, $elp_param_pls, $elp_param_plu, $elp_param_pln,
    $elp_param_d, $elp_param_l, $elp_param_f ],
  pla2 => [ $elp_param_plme, $elp_param_plv, $elp_param_plt, $elp_param_plma,
    $elp_param_plj, $elp_param_pls, $elp_param_plu, $elp_param_d, $elp_param_lp,
    $elp_param_l, $elp_param_f ],
);

my @elp82_error;
my @elp82_data;

sub elp82_parse($$$$) {
  my ($file, $part, $coord, $deg) = @_;
  my $fn = "elp82/ELP$file";
  open my $f, "<", $fn or die "$fn: $!\n";
  my $errf = ($time_span / 100) ** ($deg || 0);
  my $errt = $coord == 2 ? $elp_accuracy_radius : $elp_accuracy_angle;
  scalar <$f>;
  my $ni;
  my $pa;
  my $la;
  my $pphi;
  my $lphi;
  if($part eq "main") {
    $ni = 4;
    $pa = 14;
    $la = 13;
  } elsif($part eq "pert") {
    $ni = 5;
    $pa = 26;
    $pphi = 16;
    $la = $lphi = 9;
  } elsif($part eq "pla1" || $part eq "pla2") {
    $ni = 11;
    $pa = 44;
    $pphi = 34;
    $la = $lphi = 9;
  } else {
    die;
  }
  $deg ||= 0;
  my @factors = @{$elp_params{$part}};
  die unless $ni == @factors;
  while(<$f>) {
    my $a = getnum $pa, $la;
    my $phi = defined $pphi ? getnum($pphi, $lphi) : $coord == 2 ? 90 : 0;
    $phi *= PI / 180;
    my $errv = abs($a) * $errf;
    if($errv >= $errt) {
      my @r;
      for my $i (0 .. $ni - 1) {
	push @r, getnum $i * 3, 3;
      }
      my @k;
      for my $i (0 .. $ni - 1) {
	my $coef = $r[$i];
	next if $coef == 0;
	my $term = $factors[$i];
	#print dms(93, 16, 19.55755), "\n";
	#print "$coef -> @$term (@$elp_param_f)\n";
	for my $j (0 .. $#$term) {
	  $k[$j] += $coef * $term->[$j];
	}
      }
      push @{$elp82_data[$coord][$deg]}, [ $a, $phi, @k ];
    } else {
      $elp82_error[$coord] += $errv;
    }
  }
}

sub elp82_parse3($$$) {
  my ($file, $part, $deg) = @_;
  for my $v (0 .. 2) {
    elp82_parse $file + $v, $part, $v, $deg;
  }
}

elp82_parse3  1, "main", undef;
elp82_parse3  4, "pert", 0;
elp82_parse3  7, "pert", 1;
elp82_parse3 10, "pla1", 0;
elp82_parse3 13, "pla1", 1;
elp82_parse3 16, "pla2", 0;
elp82_parse3 19, "pla2", 1;
elp82_parse3 22, "pert", 0;
elp82_parse3 25, "pert", 1;
elp82_parse3 28, "pert", 0;
elp82_parse3 31, "pert", 0;
elp82_parse3 34, "pert", 2;

my @smart97_data;
my @smart97_error;

for my $coord (0 .. 2) {
  my $filename = "smart97/figdyn97.$smart_variables[$coord]";
  open my $f, "<", $filename or die "$filename: $!\n";
  while(<$f>) {
    my $deg = getnum 3, 2;
    my $a = getnum 85, 24;
    my $b = getnum 109, 13;
    my $c = getnum 122, 18;
    my $errv = abs($a) * ($time_span / 1000) ** $deg;
    if($errv >= $smart_accuracy) {
      push @{$smart97_data[$coord][$deg]}, [ $a, $b, $c ];
    } else {
      $smart97_error[$coord] += $errv;
    }
  }
}

print <<EOF ;
/*
This file was automatically generated from the VSOP87 and ELP82 data that
once were available on the FTP site of the Bureau des longitudes.

The format is the following:

For vsop87:

  An object, with one property per planet.

  Each property is a 3-elements array, for the 3 coordinates (x, y, z; (x,y)
  being the ecliptic plane).
  For each coordinate, an array indexed by the degree.
  For each degree, an array of terms.
  Each term is an array [ A, B, C ].
  The result for each coordinate is:
    sum t^degree A cos(B + C t)
  The unit for t is one thousand Julian years, with origin at J2000.
  The unit for the result is one astronomical unit.

For elp82

  A 3-elements array, for the 3 coordinates.
  For each coordinate, an array indexed by the degree.
  For each degree, an array of terms.
  Each term is an array [ A, phi, i[0], i[1]... ]
  The result for each coordinate is:
    sum t^degree A sin(param)
    param = sum i[j] t^j
  The unit for t is one hundred Julian years, with origin at J2000.
  The coordinates are longitude, latitude and distance.
  The unit for longitude and latitude is one arc second.
  The unit for distance is one kilimeter.
  To get absolute rectangular coordinates (the same as VSOP), the formulas are:
    x = distance cos(latitude) cos(longitude + w1)
    y = distance cos(latitude) sin(longitude + w1)
    z = distance sin(latitude)
  w1 is a polynomial in t; its coefficients are elp82_w1, in radians.

For smart97

  A 3-elements array, for the 3 angles.
  For each angle, an array indexed by the degree.
  For each degree, an array of terms.
  Each term is an array [ A, B, C ]
  The result for each angle is
    sum t^degree A cos(B + C t)
  The unit for t is one thousand Julian years, with origin at J2000.
  The unit for the result is one microarcsecond.
  The three angles are psi, omega, phi.
  The correspondance with <URL: http://en.wikipedia.org/wiki/Euler_Angles > is:
  phi -> alpha
  omega -> beta
  psi -> gamma

*/
EOF

print "var vsop87 = {\n";
for my $planet (@vsop_planets) {
  print "  $planet: [\n";
  for my $coord (0 .. 2) {
    print "    [ // error ", $vsop_data{$planet}[$coord]{error}, "\n";
    my @terms = @{$vsop_data{$planet}[$coord]{terms}};
    for my $ts (@terms) {
      print "      [\n";
      for my $t (@$ts) {
	print "        [ ", join(", ", @$t), " ],\n";
      }
      print "      ],\n";
    }
    print "    ],\n";
  }
  print "  ],\n";
}
print "};\n";

print "\n";

print "var elp82_w1 = [ ", join(", ", @$elp_param_w1), " ];\n";
print "var elp82 = [\n";
for my $coord (0 .. 2) {
  print "  [ // error ", $elp82_error[$coord], "\n";
  my $terms = $elp82_data[$coord];
  for my $deg (0 .. $#$terms) {
    print "    [\n";
    for my $term (@{$terms->[$deg]}) {
      print "      [ ", join(", ", @$term), " ],\n";
    }
    print "    ],\n";
  }
  print "  ],\n";
}
print "];\n";

print "\n";

print "var smart97 = [\n";
for my $coord (0 .. 2) {
  print "  [ // error ", $smart97_error[$coord], "\n";
  my $terms = $smart97_data[$coord];
  for my $deg (0 .. $#$terms) {
    print "    [\n";
    for my $term (@{$terms->[$deg]}) {
      print "      [ ", join(", ", @$term), " ],\n";
    }
    print "    ],\n";
  }
  print "  ],\n";
}
print "];\n";