#!perl
use 5.010;
use strict;
use warnings;
use Getopt::Long;

# Compiles the eight VSOP87D files (Bretagnon and Francou 1988, CDS VI/81)
# into include/tks/tks_vsop87d.h: every term, for every planet, coordinate
# (L, B, R) and power of time, as the A, B and C of A cos(B + C t) with the
# file's own decimal digits, in the file's order (largest first). Shipped so a
# reader can regenerate the table; never installed.
#
#     perl tools/vsopcompile --dir ../plan_time_keeper_sky/sources/vsop87 --out include/tks/tks_vsop87d.h

my ($dir, $out) = ('', 'include/tks/tks_vsop87d.h');
GetOptions('dir=s' => \$dir, 'out=s' => \$out) or die "usage: vsopcompile --dir <directory of VSOP87D.*> [--out <header>]\n";
die "vsopcompile: --dir <directory> is needed\n" unless $dir;

my @planets = (
    [mer => 'MERCURY'], [ven => 'VENUS'], [ear => 'EARTH'], [mar => 'MARS'],
    [jup => 'JUPITER'], [sat => 'SATURN'], [ura => 'URANUS'], [nep => 'NEPTUNE'],
);
my @coord = qw(l b r);
my $num = qr/-?(?:\d+\.\d*|\.\d+)/;

my (@arrays, @index, $total);
for my $pi (0 .. $#planets) {
    my ($ext, $name) = @{ $planets[$pi] };
    my $file = "$dir/VSOP87D.$ext";
    open my $fh, '<', $file or die "vsopcompile: $file: $!\n";
    my ($cur, $want, @terms, %series);
    my $close = sub {
        return unless $cur;
        die "vsopcompile: $file $cur: $want terms announced, " . scalar(@terms) . " read\n" unless @terms == $want;
        $series{$cur} = [@terms];
        @terms = ();
    };
    while (my $line = <$fh>) {
        if ($line =~ /^\s*VSOP87 VERSION D4\s+(\S+)\s+VARIABLE (\d) \(LBR\)\s+\*T\*\*(\d)\s+(\d+) TERMS?\b/) {
            my ($body, $var, $pow, $n) = ($1, $2, $3, $4);
            die "vsopcompile: $file names $body, not $name\n" unless $body eq $name;
            $close->();
            $cur = "$var.$pow";
            $want = $n;
            next;
        }
        next unless $line =~ /\S/;
        my @f = split ' ', $line;
        die "vsopcompile: $file: a term before any header\n" unless $cur;
        my ($a, $b, $c) = @f[-3 .. -1];
        die "vsopcompile: $file: '$line' does not end in A B C\n" unless "$a $b $c" =~ /^$num $num $num$/;
        push @terms, [$a, $b, $c];
    }
    $close->();
    close $fh;
    for my $var (1 .. 3) {
        for my $pow (0 .. 5) {
            my $t = $series{"$var.$pow"};
            my $id = "tks_vsop_${ext}_$coord[$var - 1]$pow";
            if ($t) {
                push @arrays, "static const tks_vsop_term ${id}[" . scalar(@$t) . "] = {\n"
                    . join('', map { "    { $_->[0], $_->[1], $_->[2] },\n" } @$t) . "};\n";
                push @index, "{ $id, " . scalar(@$t) . " }";
                $total += @$t;
            }
            else {
                push @index, '{ NULL, 0 }';
            }
        }
    }
}

open my $oh, '>', $out or die "vsopcompile: $out: $!\n";
print $oh <<"EOT";
#ifndef TKS_VSOP87D_H
#define TKS_VSOP87D_H

/* Generated by tools/vsopcompile from the eight VSOP87D files of CDS
 * catalogue VI/81 (Bretagnon and Francou 1988): every term of the
 * heliocentric ecliptic longitude, latitude and radius of Mercury to Neptune,
 * referred to the ecliptic and equinox of date, $total terms in all. Each is
 * A cos(B + C t), A in radians or AU, B in radians, C in radians a Julian
 * millennium, t in Julian millennia of TDB from J2000.0; the digits are the
 * files' own. Do not edit: regenerate. */

#define TKS_VSOP_TERMS $total

EOT
print $oh @arrays;
print $oh "\n/* [planet][coordinate L, B, R][power of t, 0 to 5] */\n";
print $oh "static const tks_vsop_series tks_vsop87d[8][3][6] = {\n";
for my $pi (0 .. 7) {
    print $oh "    {\n";
    for my $var (0 .. 2) {
        print $oh "        { ", join(', ', @index[$pi * 18 + $var * 6 .. $pi * 18 + $var * 6 + 5]), " },\n";
    }
    print $oh "    },\n";
}
print $oh "};\n\n#endif\n";
close $oh;
printf "%s: %d terms\n", $out, $total;
