#! /usr/bin/env perl

# $Id$

#use lib '/home/yamada/prj/panstarrs/mops/src';

use warnings;
use strict;

use Getopt::Long;
use Pod::Usage;
use FileHandle;
use Astro::SLA;

use PS::MOPS::DET::table;
use PS::MOPS::DET::util;
use PS::MOPS::Lib qw(:all);
use PS::MOPS::Constants qw(:all);
use PS::MOPS::DC::Instance;
use PS::MOPS::DC::Field;

my @default_de = (0, 0, 0, 0, 0, 0, 0, 0, 0, 0);

my $help;
my $inst;
my $instance_name;
my $survey = 'DRM1';            # default survey directory to load from
my $pairs;
my $noinsert;
my $max_exposure_time = 65;     # max exposure time allowed for TTI field
my $obscode = '566';            # Haleakala
my $limiting_mag = 22.7;        # PS-1
my $filter_replace = 'r';       # replace input filters in simulation
my $filter_ignore = 'z,y';      # ignore fields with these filters
my $mjd;
my $start_mjd;
my $end_mjd;

GetOptions(
    'instance=s' => \$instance_name,
    help => \$help,
    'survey=s' => \$survey,
    'max_exposure_time=f' => \$max_exposure_time,
    pairs => \$pairs,
    noinsert => \$noinsert,
    'obscode=s' => \$obscode,
    'limiting_mag=f' => \$limiting_mag,
    'filter_replace=s' => \$filter_replace,
    'filter_ignore=s' => \$filter_ignore,
    'mjd=f' => \$mjd,
    'start_mjd=f' => \$start_mjd,
    'end_mjd=f' => \$end_mjd,
) or pod2usage(2);                      # opts parse error
pod2usage(-verbose => 3) if $help;      # user asked for help

if ( $#ARGV < 0 ) {
	pod2usage(2);
	exit
}

$inst = PS::MOPS::DC::Instance->new(DBNAME => $instance_name);
my $cfg = $inst->getConfig;
my $gmt_offset_hours = $cfg->{site}->{gmt_offset_hours};


# Create table of filters to ingore.
my %filter_ignore_table;
if ($filter_ignore) {
    my @filts = split ',', $filter_ignore;
    $filter_ignore_table{$_} = 1 foreach @filts;    # set keys in table
}

sub open_fields
{
	my ($field_file, %ocnum_map) = @_;
    my $fh = new FileHandle $field_file or die "Cannot open $field_file";
    my @fields;
    my $line;

    # First line contains tab-delimited mysql fields.
	my $fields_line = <$fh>;
    chomp $fields_line;
    my @field_names = split /\s+/, $fields_line;
    my $go;         # whether to process this field
    my $field_mjd;
    my $field_ocnum;

    while (defined($line = <$fh>)) {
        chomp $line;
        my @vals = split /\s+/, $line;
        my %cols;
        @cols{@field_names} = @vals;

        $field_mjd = int($cols{expMJD});
        $field_ocnum = mopslib_mjd2ocnum($field_mjd);

        next if $filter_ignore_table{$cols{filter}};

        if ($cols{propProgram} =~ /AP|SS/ and $ocnum_map{$field_ocnum}) {
            push @fields, {
                field_id => $cols{obsHistID}, 
                epoch_mjd => $cols{expMJD}, 
                ra_deg => $cols{fieldRA}, 
                dec_deg => $cols{fieldDec}, 
                survey_mode => $cols{propProgram}, 
                time_start => $cols{expMJD} - $cols{expTime} / 86400 / 2,
                time_stop => $cols{expMJD} + $cols{expTime} / 86400 / 2, 
                filter_id => ($filter_replace or $cols{filter}), 
                limiting_mag => $limiting_mag,
                ra_sigma_deg => 0,
                dec_sigma_deg => 0,
                obscode => $obscode, 
                de1 => 0,
                de2 => 0,
                de3 => 0,
                de4 => 0,
                de5 => 0,
                de6 => 0,
                de7 => 0,
                de8 => 0,
                de9 => 0,
                de10 => 0,
                ocnum => mopslib_mjd2ocnum($cols{expMJD}), 
                status => $FIELD_STATUS_NEW,
            };
        }   # if
    }   # while
    close $fh;
    my @sorted_fields = sort { $a->{epoch_mjd} <=> $b->{epoch_mjd} } @fields;
    return \@sorted_fields;
}


sub _sep_rad {
    my ($f1, $f2) = @_;
    return slaDsep(
        $f1->{ra_deg} / $DEG_PER_RAD, $f1->{dec_deg} / $DEG_PER_RAD,
        $f2->{ra_deg} / $DEG_PER_RAD, $f2->{dec_deg} / $DEG_PER_RAD,
    );
}


sub filter_tti_pairs {
    my @inp_fields = @_;
    my @paired_fields;
    my ($f, $f1, $match, $idx);     # field objects
    my ($exp_mjd, $exp1_mjd);       # exposure times

    my $MAX_EXPOSURE_TIME_MJD = $max_exposure_time / 86400;         # max exposure time in MJD
    my $ONE_HOUR_MJD = 3600 / 86400;                # one hour in MJD
    my $MAX_SEP_RAD = 0.1 / $DEG_PER_RAD;           # .1 degree max spatial separation of TTI pair

    while (@inp_fields) {
        $f = shift @inp_fields;     # get first in list

        $exp_mjd = $f->{time_stop} - $f->{time_start};
        next if $exp_mjd == 0 or $exp_mjd > $MAX_EXPOSURE_TIME_MJD;    # toss bogus exp time or MD fields

        $idx = 0;
        $match = undef;
        foreach $f1 (@inp_fields) {
            $exp1_mjd = $f1->{time_stop} - $f1->{time_start};
            if (
                $f1->{filter_id} eq $f->{filter_id}     # filters must match
                and $exp1_mjd != 0 
                and abs($f->{epoch_mjd} - $f1->{epoch_mjd}) > $ONE_HOUR_MJD
                and (_sep_rad($f, $f1) < $MAX_SEP_RAD)
            ) {
                # If we've made it this far, the field is acceptable.  Pair $f and $f1 together, and remove
                # $f1 from the current list.
                $match = $f1;
                last;   # outta here
            };
            $idx++;
        }
        if ($match) {
            splice @inp_fields, $idx, 1;        # extract item, compressing list
            push @paired_fields, $f, $match;
        }
    }

    return @paired_fields;
}


sub process
{
	my ($src, %ocnum_map) = @_;

    # Convert src to a filename.
    my $src_filename = "$ENV{MOPS_HOME}/data/ssm/surveys/$src/obshist.dump";
    die "can't find file $src_filename" unless -f $src_filename;
	print STDERR "Reading all fields from $src_filename.\n";

	my $fields_list = open_fields($src_filename, %ocnum_map);
    printf STDERR "Found %d fields for OC list %s.\n", scalar @{$fields_list}, join(' ', sort keys %ocnum_map);

    # Bin our data into nights.  Then for each night, locate TTI pairs.
    my %seen_fields;            # record of what epoch/RA/DECs we've inserted

    while (@{$fields_list}) {
        # Now filter so that we have only TTI pairs.
        my @nightly_fields;
        my $f;
        my $nn;
        $f = shift @{$fields_list};         # first item
        $nn = mopslib_mjd2nn($f->{epoch_mjd}, $gmt_offset_hours);
        push @nightly_fields, $f;
        while (@{$fields_list} and mopslib_mjd2nn($fields_list->[0]->{epoch_mjd}, $gmt_offset_hours) == $nn) {
            push @nightly_fields, shift @{$fields_list};        # transfer $fields_list->[0] => $nightly_fields
        }

        # If pairs is specified, filter out only apparent TTI pairs.
        if ($pairs) {
            @nightly_fields = filter_tti_pairs(@nightly_fields);
            printf STDERR "Found %d fields in TTI pairs for night $nn.\n", scalar @nightly_fields;
        }

        next if $noinsert;

        # Insert into DB.
        foreach my $field (@nightly_fields) {
            my $id = $field->{field_id};
#            my $dest = (sprintf "%012d", $field->{field_id}) . ".fits";
#            my $p = new PS::MOPS::DET::table;
            my $epoch = ($field->{'time_start'} + $field->{'time_stop'}) / 2.0; 

            # See if this field has been added already.  The scheduler's obsHistory
            # reports a field multiple times if it satisfies multiple proposals.
            my $key = join(':', $epoch, $field->{ra_deg}, $field->{dec_deg});
            if ($seen_fields{$key}) {
                printf STDERR "Found duplicate field: $key.\n";
                next;
            }
            $seen_fields{$key} = 1;         # record that we've seen this

            my $mops_field = PS::MOPS::DC::Field->new($inst, 
                epoch => $epoch,
                ra => $field->{ra_deg},
                dec => $field->{dec_deg},
                timeStart => $field->{time_start},
                timeStop => $field->{time_stop},
                filter => $field->{filter_id},
                surveyMode => $field->{survey_mode},
                limitingMag => $field->{limiting_mag},
                raSigma => 0,
                decSigma => 0,
                obscode => $field->{obscode},
                de => \@default_de
            );
            $id = $mops_field->insert;
#            printf STDERR "created field id $id (%.6f).\n", $field->epoch;
        }
        my $ocnum = mopslib_mjd2ocnum($nn);
        printf STDERR "$ocnum/$nn: inserted %d fields.\n", scalar @nightly_fields;
    }   # bin
}

my %ocnum_map;
$ocnum_map{$_} = 1 foreach (@ARGV);     # create hash with OCNUMs as keys
process($survey, %ocnum_map);
exit;

__END__

=head1 NAME

insertObsHist - Insert fields from obsHist dump into current MOPS simulation.

=head1 SYNOPSIS

insertObsHist [options] OCNUM

  --survey NAME : get survey data from specified filename, default $MOPS_HOME/data/ssm/surveys/drm1-default/obshist.dump
  --max_exposure_time T : don't include fields with exposure time > T sec
  --pairs : only import fields that constitute TTI pairs
  --obscode OBSCODE : set field observatory codes to OBSCODE, default 566 (PS1)
  --limiting_mag MAG : set limiting magnitude, default 22.7 (PS1)
  --filter_replace FILT : replace all filter values (grizyw) with FILT , default "r"
  --filter_ignore FILT : ignore fields with filters in FILT (default "z,y")
  --noinsert : don't insert anything into the database, just scan input file
  --help : show help page

insertObsHist 83        # insert all fields from OCNUM 83

insertObsHist --survey=spacewatch --max_exposure_time=120 --pairs --obscode=691

=head1 DESCRIPTION

Insert ObsHistory database dump from scheduler into MOPS instance.  If --survey=NAME is specified,
it is used to select a the file $MOPS_HOME/data/ssm/surveys/NAME/obshist.dump.

=head1 ENVIRONMENT VARIABLES

=head1 TO DO

=head1 BUGS

=cut
