#! /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 File::Copy;
use File::Temp;
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;

use subs qw(
    open_fields
    filter_tti_pairs
    process
    submit_to_datastore
    false_detection_rate
    false_detection_mag
    do_field_false_detections
);


my @default_de = (0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
my $default_ds_type = 'DETECTION';         # DS file type

my $help;
my $inst;
my $instance_name;
my $obshist_file;
my $create_false;
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;
my $ocnum_flag;
my $nn_flag;


GetOptions(
    'instance=s' => \$instance_name,
    help => \$help,
    'obshist_file=s' => \$obshist_file,
    create_false => \$create_false,
    '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,
    ocnum => \$ocnum_flag,
    nn => \$nn_flag,
) or pod2usage(2);                      # opts parse error
pod2usage(-verbose => 3) if $help;      # user asked for help

if (@ARGV == 0) {
	pod2usage(2);
	exit;
}

my %ocnum_map;
my %nn_map;
if ($ocnum_flag) {
    $ocnum_map{$_} = 1 foreach (@ARGV);     # create hash with OCNUMs as keys
}
elsif ($nn_flag) {
    $nn_map{$_} = 1 foreach (@ARGV);        # create hash with NNs as keys
}
else {
	pod2usage(2);
	exit;
}



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

my $gmt_offset_hours = $mops_config->{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
}



# False detection stuff.
my $_log10 = log(10);       # used to calculate log10

# Configuration.
my $field_size_deg2 = $mops_config->{site}->{field_size_deg2}
    or die("couldn't get field_size_deg2");
my $field_shape = $mops_config->{site}->{field_shape} || 'circle';

my $astrometric_error_arcseconds = $mops_config->{site}->{astrometric_error_arcseconds}
    or die("couldn't get astrometric_error_arcseconds");

# Deprecated.  Use configuration FWHM and S2N to compute false detections/deg2.
my $false_detection_s2n = $mops_config->{site}->{s2n} || 5;  # default 5-sigma

my $fwhm_arcseconds = $mops_config->{site}->{fwhm_arcseconds}
    or die("couldn't get fwhm_arcseconds");
my $false_detection_rate_per_deg2 = false_detection_rate(
    $false_detection_s2n, $fwhm_arcseconds
);
my $num_false = int($field_size_deg2 * $false_detection_rate_per_deg2);



#my $obshist_file = "/usr/local/MOPS/data/ssm/surveys/$survey/obshist.dump";
process($obshist_file, \%ocnum_map, \%nn_map);
exit;


sub open_fields
{
	my ($field_file, $oc_mapref, $nn_mapref) = @_;
    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_nn;
    my $field_ocnum;

    # Compute max MJD to consider so we can stop parsing big file.
    my $max_mjd;
    if ($nn_mapref) {
        $max_mjd = mopslib_nn2mjd(
            (sort { $b <=> $a } keys %{$nn_mapref})[0],
            $mops_config->{site}->{gmt_offset_hours}
        ) + 1;
    }
    elsif ($oc_mapref) {
        $max_mjd = mopslib_ocnum2mjd(
            (sort { $b <=> $a } keys %{$oc_mapref})[0],
        ) + 1;
    }

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

        $field_mjd = $cols{expMJD};
        $field_nn = mopslib_mjd2nn($field_mjd, $mops_config->{site}->{gmt_offset_hours});
        $field_ocnum = mopslib_mjd2ocnum($field_mjd);

        last if $field_mjd > $max_mjd;
        next if $filter_ignore_table{$cols{filter}};

        if ($cols{propProgram} =~ /AP|SS/ and (exists(${$oc_mapref}{$field_ocnum}) or exists(${$nn_mapref}{$field_nn}))) {
            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_FALSEDETECTIONS,
            };
        }   # if
    }   # while
    close $fh;
    my @sorted_fields = sort { $a->{epoch_mjd} <=> $b->{epoch_mjd} } @fields;
    printf STDERR "Found %d fields.\n", scalar @sorted_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_filename, $oc_mapref, $nn_mapref) = @_;


    # Initialize datastore.
    submit_to_datastore();


    # Convert src to a filename.
    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, $oc_mapref, $nn_mapref);

    # 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 $tmpfile = File::Temp::tempnam('/tmp', 'obsHist2DS-XXXXXXXX');
            my $p = new PS::MOPS::DET::table;
            my $epoch = ($field->{'time_start'} + $field->{'time_stop'}) / 2.0; 

            my $force = 1;      # force file overwrite
            $p->createfile($tmpfile, $force,
                'FPA_ID' => $field->{'field_id'},
                'MJD-OBS' => $epoch,
                'RA' => PS::MOPS::DET::util::DDtoHMS($field->{'ra_deg'}),
                'DEC' => PS::MOPS::DET::util::DDtoDMS($field->{'dec_deg'}),
                'EXPTIME' => $field->{'time_stop'} - $field->{'time_start'},
                'ROTANGLE' => 0.0,
                'FILTER' => $field->{'filter_id'},
                'STARPSF' => 0.0,
                'LIMITMAG' => $field->{'limiting_mag'},
                'DE1' => $field->{'de1'},
                'DE2' => $field->{'de2'},
                'DE3' => $field->{'de3'},
                'DE4' => $field->{'de4'},
                'DE5' => $field->{'de5'},
                'DE6' => $field->{'de6'},
                'DE7' => $field->{'de7'},
                'DE8' => $field->{'de8'},
                'DE9' => $field->{'de9'},
                'DE10' => $field->{'de10'},
                'OBSCODE' => $field->{'obscode'},
                'TEL_ALT' => 0.0,
                'TEL_AZ' => 0.0);

            if ($create_false) {
                do_field_false_detections(
                    $p, $field->{ra_deg}, $field->{dec_deg}, 'DUMMY_FIELD_ID', $field->{epoch_mjd}, $field->{filter_id},
                    $field->{limiting_mag}, $field->{obscode}, $num_false, $field_size_deg2
                );
            }


            $p->closefile();
            print STDERR "wrote $tmpfile ($epoch).\n";

            # Now install the file into the datastore.
            submit_to_datastore($tmpfile);
            unlink $tmpfile;
        }

        my $ocnum = mopslib_mjd2ocnum($nn);
        printf STDERR "$ocnum/$nn: inserted %d fields.\n", scalar @nightly_fields;
    }   # bin
}


sub submit_to_datastore {
	my ($src_file) = @_;
	my ($destdir) = '/usr/local/MOPS_DEV/ds/dsroot/det';
    my $dbname = $inst->dbname();

	my $fh = new FileHandle;

	my ($sec,$min,$hour,$mday,$mon,$year) = gmtime(time);
	$mon++;
	$year += 1900;

    my $state = "N";

    # See if the directory exists.  If not, create it and initialize
    # the count file and index.
    my $dbdir = $destdir . "/" . $dbname;
    my $index = $destdir . "/index";
    my @dbstat = stat($dbdir);

    unless (@dbstat) {
        $mops_logger->info("Initializing datastore.");

        # Create the per-instance directory.
        mkdir($dbdir, 0777);
        chmod(0777, $dbdir);

        open($fh, ">$dbdir" . "/count");
        flock($fh, 2);
        chmod(0666, $dbdir . "/count");
        print $fh "0\n";
        close $fh;

        open($fh, ">$dbdir" . "/index");
        flock($fh, 2);
        chmod(0666, $dbdir . "/index");
        close $fh;

        # Create an entry in the root index.
        open($fh, ">>$index");
        flock($fh, 2);
        chmod(0666, $index);
        printf($fh "%s|%04d-%02d-%02dT%02d:%02d:%02d|%s\n",
            $dbname,
            $year, $mon, $mday, $hour, $min, $sec, 'INSTANCE');
        close $fh;
    }

    # If $src_file was specified, copy it to the DataStore.  This gives us a cheap
    # way of initializing the sim's datastore area by passing $src_file=undef.
    if ($src_file) {
        # Read and increment the count number.
        my $count;

        open($fh, "+<$dbdir/count");
        flock($fh, 2);
        chmod(0666, "$dbdir/count");
        $count = <$fh>;
        chomp($count);
        seek($fh, 0, 0);
        print $fh $count + 1 . "\n";
        close $fh;

        # Send the file.
        my $output_file = sprintf("%s/%.05d", $dbdir, $count);
        copy($src_file, $output_file);

        open($fh, ">>$dbdir/index");
        flock($fh, 2);
        printf($fh "%.05d|%04d-%02d-%02dT%02d:%02d:%02d|%s|%s\n",
            $count,
            $year, $mon, $mday, $hour, $min, $sec,
            $default_ds_type, $state);
        close($fh);
    }
}


=head1 NAME

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

=head1 SYNOPSIS

insertObsHist [options] OCNUM

  --obshist_file FILE : get survey data from specified filename, default $MOPS_VAR/run/obshist
  --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







sub false_detection_mag {
=pod

As per our lunch discussion yesterday, I've worked out how to generate fake 
fluxes for the false detections...

According to eq. 153 from Kaiser's PSDC-200-010-00 the number distn of false 
detections goes like

f(x) ~ x * exp( -x*x/2 )

where x=S/N.

Let F(y) be the normalized fraction of objects in the range [x_min,y] or:

       \int_{x_min}^{y}    x * exp( -x*x/2 ) dx
F(y) = -------------------------------------
       \int_{x_min}^{\inf} x * exp( -x*x/2 ) dx

     = 1 - exp[ ( - ( y^2 - x_min^2 ) / 2 ] 

so, if we generate a random number p in the range [0,1] then the corresponding 
x=S/N is given by

x = S/N = sqrt( x_min^2 - 2 ln( 1 - p ) )

We want x_min=5.0 at V=24.0 so to convert from x to V use:

V = 24.0 - 2.5 log( x/5 )


SUMMARY:

For each randomly generated false detection on the image:

1) generate p randomly in range [0,1)
2) get corresponding value of S/N (= sqrt( 25 - 2 ln( 1 - p ) ))
  (let F(y) = p; solve for y)
3) convert to V = 24.0 - 2.5 log( x/5 )


I'm sure one of you will let me know if there is a problem...

=cut

	my ($limiting_mag) = @_;
    my $p;
    my $SN;
    my $V;

    $p = rand(1);       # uniform random number in [0, 1)
#    $SN = sqrt($x_min * $x_min - 2 * log(1 - $p));   # 25 = 5 sigma ^2 (?)
    $SN = sqrt(25 - 2 * log(1 - $p));   # 25 = 5 sigma ^2 (?)
    $V = $limiting_mag - 2.5 * log($SN / 5) / $_log10;
    return ($V, $SN);
}


sub false_detection_rate {
=pod

    PANSTARRS#S2N-FALSE-DETECTIONS  FWHM=0.7

    * Returns number of false detections per deg^2 of S/N=x in PS images
    * according to eq. 153 from Kaiser's PSDC-200-010-00.

    * FWHM is the PSF guassian FWHM in arcseconds.

    * To use FWHM instead of SIGMA use SIG2FWHM=5.2255
    SIGMAG=[FWHM]/2.35

    * AS2PERDEG2=1.296e7 = ARC-SECOND^2 PER DEGREE^2
    AS2PERDEG2=1.296e7

    * Hokey-Correction Factor (HCF) to get 250 detections/deg^2 @ S/N=5
    * NOTE: is HCF due to using sigma instead of real FWHM as input parameter?
    HCF=2.9

    * 2^(5/2)*pi^(3/2)=31.49921989
    C=[HCF]*[AS2PERDEG2]/31.499/[SIGMAG]/[SIGMAG]

    NPERDEG2 = [C]*S2N*exp(-S2N*S2N/2.0)

    RETURN

=cut

    my ($s2n, $fwhm) = @_;
    my $SIGMAG = $fwhm / 2.35;
#    $HCF = 2.9;                     # hokey correction factor
#    $ARCSEC2_PER_DEG2 = 1.296e7;
#    $SOMETHING = 31.49921989;
#    $c = $HCF * $ARCSEC2_PER_DEG2 / $SOMETHING / ($SIGMAG * $SIGMAG);
    my $c = 2.9 * 1.296e7 / 31.49921989 / ($SIGMAG * $SIGMAG);
    return int($c * $s2n * exp(-$s2n * $s2n / 2));
}


sub do_field_false_detections {
    my ($fits_obj, $field_ra_deg, $field_dec_deg, $field_id, $epoch, $filter,
		$limiting_mag, $obscode, $num, $field_size_deg2) = @_;
    my $i = 0;          # count inserted fields
    my $det_num = 0;    # night-unique false det ID

	# For the fieldId provided, generate $num uniformly-distributed
	# false detections in the field.
	eval {
		my ($tpra_rad, $tpdecl_rad);  # tangent plane RA/DEC in RAD
		my ($fra_rad, $fdecl_rad);  # false RA, DEC in RAD
		my ($fra_deg, $fdecl_deg);  # false RA, DEC in DEG
		my $tp_size_rad;            # size of tangent plane area in rad
		my $mag;
		my $ra_sigma_deg;
		my $dec_sigma_deg;
		my $mag_sigma;
		my $S2N;
        my $flux;
        my $flux_sigma;
        my $orient_deg = 0;
        my $orient_sigma_deg = 0;
        my $length_deg = 0;
        my $length_sigma_deg = 0;
        my $in_field;


        my $fov_size_deg;
        if ($field_shape eq 'circle') {
            $fov_size_deg = 2 * sqrt($field_size_deg2 / $PI);  # width of field
        }
        else {
            $fov_size_deg = sqrt($field_size_deg2);            # width of field
        }
        $tp_size_rad = $fov_size_deg / $DEG_PER_RAD;

		for my $n (0..$num) {
			$tpra_rad = (rand() - 0.5) * $tp_size_rad;         # tangent plane coords in rad
			$tpdecl_rad = (rand() - 0.5) * $tp_size_rad;     # tangent plane coords

			slaDtp2s(
				$tpra_rad, $tpdecl_rad, 
				$field_ra_deg / $DEG_PER_RAD, $field_dec_deg / $DEG_PER_RAD,
				$fra_rad, $fdecl_rad
			);
			$fra_deg = $fra_rad * $DEG_PER_RAD;
			$fdecl_deg = $fdecl_rad * $DEG_PER_RAD;

            # Now make sure the detection is actually in the field. If so, insert.
            $in_field = 0;
            if ($field_shape eq 'circle') {
                if (mopslib_inField($field_ra_deg, $field_dec_deg, $fov_size_deg, $fra_deg, $fdecl_deg)) {
                    $in_field = 1;
                }
            }
            else {
                if (mopslib_inSquareField($field_ra_deg, $field_dec_deg, $fov_size_deg, $fra_deg, $fdecl_deg)) {
                    $in_field = 1;
                }
            }


            if ($in_field) {
                ($mag, $S2N) = false_detection_mag($limiting_mag);
                ($flux, $flux_sigma) = mopslib_hsn2flux($mag, $S2N, $filter);
                next if $flux == 0 or $flux_sigma == 0;

                # Compute uncertainty in position.  This is the baseline PS uncertainty +
                # uncertainty due to photometry.
                $ra_sigma_deg = ($astrometric_error_arcseconds + mopslib_astroError($mag, $S2N)) / 3600;
                $dec_sigma_deg = $ra_sigma_deg;

                $fits_obj->writerecord(
                    ($fra_deg, $fdecl_deg, 
                    $ra_sigma_deg, $dec_sigma_deg, 
                    $flux, $flux_sigma, 
                    $orient_deg, $orient_sigma_deg,
                    $length_deg, $length_sigma_deg,
                    )
                );
                $i++;
            }

		}   # for
	};

    $mops_logger->info("Inserted $i false detections.");
    warn $@ if $@;
    return $i;      # return num dets added
}

__END__
