#!/usr/bin/env perl

use strict;
use warnings;

=head1 NAME

findast - Find detectinos for numbered asteroid in MOPS DB

=head1 SYNOPSIS

findast NUMBER

  --instance=INSTANE_NAME : specify simulation to use
  --radius_arcsec=DIST_ARCSEC : search for detections with DIST arcseconds of prediction
  NUMBER : asteroid number

=head1 DESCRIPTION

Queries AstDys for ephemerides of specified asteroid at all field
times in the specified simulation, then brute-force searches for 
detections near the predictions within a specified radius.

=cut

use warnings;
use strict;

use Getopt::Long;
use Pod::Usage;
use FileHandle;
use Data::Dumper;
use LWP::UserAgent;
use Astro::Time;

use PS::MOPS::DC::Instance;
use PS::MOPS::Lib qw(:all);
use PS::MOPS::Constants qw(:all);
use PS::MOPS::DC::Field;
use PS::MOPS::DC::Detection;

use subs qw(
    get_prediction
    find_detections
    fmtf
);


my $instance_name;
my $inst;
my $radius_arcsec = 2.0;        # 2-arcsec radius default
my $show_unfound;
my $one;
my $help;

GetOptions(
    'instance=s' => \$instance_name,
    'radius_arcsec=f' => \$radius_arcsec,
    show_unfound => \$show_unfound,
    1 => \$one,
    help => \$help,
) or pod2usage(2);                      # opts parse error
pod2usage(-verbose => 3) if $help;      # user asked for help

$inst = PS::MOPS::DC::Instance->new(DBNAME => $instance_name);
my $mops_config = $inst->getConfig();
my $mops_logger = $inst->getLogger();

my $v2filt = $mops_config->{site}->{v2filt};

# Now check some args.
pod2usage(-msg => 'Asteroid number not specified.') unless @ARGV;
my $anum = shift;   # asteroid number


# First thing is to get a list of fields we need ephemerides for.
my $field_i = modcf_retrieve($inst, all => 1);
my $field;
my @fields;
while ($field = $field_i->next) {
    push @fields, $field;
    last if $one;       # hack for just one field for now
}
printf STDERR "Found %d fields.\n", scalar @fields;


# Now get the prediction for each of the fields.
print STDERR "Fetching ephemerides for $anum from AstDys...";
my %predictions;
foreach $field (@fields) {
    $predictions{$field->fieldId} = get_prediction($anum, $field);
    print STDERR '.';
}
print STDERR "done.\n";

# Dump the predictions.
if (0) {
    foreach my $key (sort {$a <=> $b} keys %predictions) { 
        my $p = $predictions{$key};
        print STDERR join(' ', 
            sprintf("%.6f", $p->{FIELD}->epoch),
            sprintf("%.6f", $p->{RA_DEG}),
            sprintf("%.6f", $p->{DEC_DEG}),
            sprintf("%.3f", $p->{MAG}),
        ), "\n";
    }
}


# Now search each field for detections matching the prediction.
print <<"HEADER";
ASTEROID FIELD_ID EPOCH_MJD DET_ID RA_DEG DEC_DEG MAG S2N PRED_RA PRED_DEC PRED_MAG DELTA_RA_ARCSEC DELTA_DEC_ARCSEC DIST_ARCSEC
HEADER
foreach my $key (sort {$a <=> $b} keys %predictions) { 
    my $p = $predictions{$key};
    my $result = find_detection($p, $radius_arcsec);
    if ($result) {
        print join(' ', 
            $anum,
            $p->{FIELD}->fieldId,
            sprintf("%.6f", $p->{FIELD}->epoch),
            @{$result}{qw(
                DET_ID
                RA_DEG
                DEC_DEG
                MAG
                S2N
                PRED_RA
                PRED_DEC
                PRED_MAG
                DELTA_RA_ARCSEC
                DELTA_DEC_ARCSEC
            )},
            $result->{DIST_ARCSEC},
        ), "\n";
    }
}

exit;


sub get_prediction {
    my ($anum, $field) = @_;
    #my $URI_TEMPLATE = 'http://hamilton.dm.unipi.it/astdys/index.php?pc=1.1.3.1&n=113763&oc=566&y0=2009&m0=6&d0=2&h0=18&mi0=59&y1=2009&m1=6&d1=2&h1=19&mi1=59&ti=1.0&tiu=days'; 
    my $uri = 'http://hamilton.dm.unipi.it/astdys/index.php?pc=1.1.3.1&n=ANUM&oc=OBSCODE&y0=YEAR1&m0=MONTH1&d0=DAY1&h0=HOUR1&mi0=MINUTE1&y1=YEAR2&m1=MONTH2&d1=DAY2&h1=HOUR2&mi1=MINUTE2&ti=1.0&tiu=days'; 

    my $ua = LWP::UserAgent->new;
    my $obscode = $field->obscode;

    # Set up time variables.
    my ($day1, $month1, $year1, $ut1) = mjd2cal($field->epoch);
    my ($sign1, $hour1, $minute1, $sec1) = time2hms($ut1, 'H', 2);

    my ($day2, $month2, $year2, $ut2) = mjd2cal($field->epoch + 0.5);
    my ($sign2, $hour2, $minute2, $sec2) = time2hms($ut2, 'H', 2);

    # Substitute our stuff.
    $uri =~ s/ANUM/$anum/g;
    $uri =~ s/OBSCODE/$obscode/g;

    $uri =~ s/YEAR1/$year1/g;
    $uri =~ s/MONTH1/$month1/g;
    $uri =~ s/DAY1/$day1/g;
    $uri =~ s/HOUR1/$hour1/g;
    $uri =~ s/MINUTE1/$minute1/g;

    $uri =~ s/YEAR2/$year2/g;
    $uri =~ s/MONTH2/$month2/g;
    $uri =~ s/DAY2/$day2/g;
    $uri =~ s/HOUR2/$hour2/g;
    $uri =~ s/MINUTE2/$minute2/g;

    my $req = HTTP::Request->new(GET => $uri);
    my $res = $ua->request($req);
    die "HTTP request failed" unless $res->is_success;

    # Now scrub for ephemeris.
    my $content = $res->content;
    if ($content =~ m|<pre.*?>(.*)</pre>|s) {
        # This is totally a hack, but if the year is present on a line, treat that as a result line.
        my @lines = grep { /\b$year1\b/ } split /\n/, $1;
        die "can't find result in HTTP response" unless @lines;
        my @crud = split /\s+/, $lines[0];
        die "empty result string" unless @crud;

        my ($dd, $dm, $ds) = @crud[8, 9, 10];

        my ($ra_deg, $dec_deg, $mag);
        $ra_deg = Astro::Time::turn2deg(Astro::Time::str2turn(join(' ', @crud[5, 6, 7]), 'H'));
        $dec_deg = Astro::Time::turn2deg(Astro::Time::str2turn(join(' ', @crud[8, 9, 10]), 'D'));
        $mag = $crud[11];

        return { FIELD => $field, RA_DEG => $ra_deg, DEC_DEG => $dec_deg, MAG => $mag };
    }
    else {
        die sprintf "did not get ephemeris for epoch %f", $field->epoch;
    }
}


sub find_detection {
    my ($p, $radius_arcsec) = @_;

    # Given the field and prediction in $p, query the database for nearby detections within
    # the specified search radius.
    my $radius_rad = $radius_arcsec / 3600 / $DEG_PER_RAD;
    my $dbh = $inst->dbh;

    my $dist_str = sprintf <<"DIST", $p->{DEC_DEG}, $p->{DEC_DEG}, $p->{RA_DEG};
/* spherical distance between fields */
acos(least(1.0, /* handle slight roundoff error causing arg > 1.0 */
    sin(radians(d.dec_deg)) * sin(radians(%f))
    + cos(radians(d.dec_deg)) * cos(radians(%f)) * cos(radians(d.ra_deg - %f))
)) 
DIST


    my $sql = sprintf <<"SQL";
select 
    d.det_id, 
    d.ra_deg, 
    d.dec_deg, 
    d.ra_sigma_deg, 
    d.dec_sigma_deg, 
    d.mag, 
    d.mag_sigma, 
    d.s2n,
    d.filter_id,
    $dist_str ddist
from detections d join fields f using (field_id) 
where $dist_str < $radius_rad and f.field_id=?
order by ddist 
limit 1
SQL

    my $aref = $dbh->selectall_arrayref($sql, undef, $p->{FIELD}->fieldId) or die $dbh->errstr;
    if (@$aref) {
        my $rref = $aref->[0];
        return {
            DET_ID => $rref->[0],
            RA_DEG => fmtf(6, $rref->[1]),
            DEC_DEG => fmtf(6, $rref->[2]),
            RA_SIGMA_DEG => fmtf(6, $rref->[3]),
            DEC_SIGMA_DEG => fmtf(6, $rref->[4]),
            MAG => fmtf(3, mopslib_filt2V($rref->[5], $rref->[8], $v2filt)),
            MAG_SIGMA => fmtf(3, $rref->[6]),
            S2N => fmtf(3, $rref->[7]),
            PRED_RA => fmtf(6, $p->{RA_DEG}),
            PRED_DEC => fmtf(6, $p->{DEC_DEG}),
            PRED_MAG => fmtf(3, $p->{MAG}),
            DELTA_RA_ARCSEC => fmtf(3, ($rref->[1] - $p->{RA_DEG}) * 3600 * cos($p->{DEC_DEG} / $DEG_PER_RAD)),
            DELTA_DEC_ARCSEC => fmtf(3, ($rref->[2] - $p->{DEC_DEG}) * 3600),
            DIST_ARCSEC => fmtf(3, $rref->[9] * $DEG_PER_RAD * 3600),
        };
    }
    else {
        return $show_unfound ? {
            DET_ID => '---',
            RA_DEG => '---',
            DEC_DEG => '---',
            RA_SIGMA_DEG => '---',
            DEC_SIGMA_DEG => '---',
            MAG => '---',
            MAG_SIGMA => '---',
            S2N => '---',
            PRED_RA => fmtf(6, $p->{RA_DEG}),
            PRED_DEC => fmtf(6, $p->{DEC_DEG}),
            PRED_MAG => fmtf(3, $p->{MAG}),
            DELTA_RA_ARCSEC => '---',
            DELTA_DEC_ARCSEC => '---',
            DIST_ARCSEC => '---',
        } : undef;
    }
}


sub fmtf {
    my ($decimals, $f) = @_;
    return sprintf("%.${decimals}f", $f);
}
