#!/usr/bin/env perl

use strict;
use warnings;

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

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

use subs qw(
    false_detection_rate
    false_detection_mag
    do_field
);

# Options.
my $inst;
my $instance_name = $ENV{'MOPS_DBINSTANCE'};
my $field_file;
my $output_prefix;
my $output_suffix = '.txt';
my $debug;
my $help;
my $no_false = 0;

# Field status management.

my $_log10 = log(10);       # used to calculate log10

# Start program here.
GetOptions(
    'instance=s' => \$instance_name,
    'fields=s' => \$field_file,
    'prefix=s' => \$output_prefix,
    'suffix=s' => \$output_suffix,
    'nofalse' => \$no_false,
    debug => \$debug,
    help => \$help,
) or pod2usage(2);                      # opts parse error
pod2usage(-verbose => 3) if $help;      # user asked for help
pod2usage(2) unless $field_file;
pod2usage(2) unless $instance_name;
pod2usage(2) unless $output_prefix;
$inst = PS::MOPS::DC::Instance->new(DBNAME => $instance_name);
my $mops_config = $inst->getConfig;

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

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->{ssm}->{false_detection_s2n} 
    || $mops_config->{ssm}->{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
);

# Start.

## Get list of fieldIds.

my $line;

my $fh = new FileHandle;

open($fh, "<$field_file") or die "Cannot open $field_file";

foreach (<$fh>) {
	$line = $_;
	chomp $line;
	my ($field_id, $epoch_mjd, $ra_deg, $dec_deg, $survey_mode, $time_start,
		$time_stop, $filter_id, $limiting_mag, $ra_sigma, $dec_sigma,
		$obscode, $de1, $de2, $de3, $de4, $de5, $de6, $de7, $de8, $de9, $de10,
		$ocnum, $status) = split(/	/, $line);

	my $num_false = int($field_size_deg2 * $false_detection_rate_per_deg2);     # total dets to create
	my $num_added;

	$num_false = -1 if $no_false;

	my $output_file_name = sprintf("%s%05d%s", $output_prefix,
		$field_id, $output_suffix);
	my $output_file = new FileHandle ">$output_file_name" or die "can't create $output_file_name";

	print $output_file $line, "\n";
	$num_added = do_field($output_file, $ra_deg, $dec_deg, $field_id,
		$epoch_mjd, $filter_id, $limiting_mag, $obscode, $num_false,
		$field_size_deg2); 
	if ($num_added > 0 and !$debug)  {
		printf "%s: Added $num_false false detections to field %d.\n",
			$output_file_name, $field_id;
	}

	close $output_file;
}

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 {
    my ($output_file, $field_ra, $field_dec, $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 ($field_dra_rad, $field_ddecl_rad);
		my $mag;
		my $ra_sigma_deg;
		my $dec_sigma_deg;
		my $mag_sigma;
		my $S2N;

		$field_dra_rad = sqrt($field_size_deg2) / $DEG_PER_RAD;    # RA size of field in RAD
		$field_ddecl_rad = sqrt($field_size_deg2) / $DEG_PER_RAD;  # DEC size of field in RAD

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

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

			($mag, $S2N) = false_detection_mag($limiting_mag);
			$mag_sigma = mopslib_photoError($mag, $S2N);

			# 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;

			my $det_id = sprintf "F%07d", $det_num++;
			my $s2n = $S2N;
			my $is_synthetic = 0;
			my $orient_deg = 0;
			my $orient_sigma_deg = 0;
			my $length_deg = 0;
			my $length_sigma_deg = 0;
			my $object_name = $MOPS_NONSYNTHETIC_OBJECT_NAME;
			my $status = 'F';
			my $ref_mag = 0;
			my $flux = 0;
			my $flux_sigma = 0;

			print $output_file
				$det_id, "\t",
				$field_id, "\t",
				$fra_deg, "\t",
				$fdecl_deg, "\t",
				$epoch, "\t",
				$mag, "\t",
				$ref_mag, "\t",
				$flux, "\t",
				$filter, "\t",
				$s2n, "\t",
				$ra_sigma_deg, "\t",
				$dec_sigma_deg, "\t",
				$mag_sigma, "\t",
				$flux_sigma, "\t",
				$orient_deg, "\t",
				$orient_sigma_deg, "\t",
				$length_deg, "\t",
				$length_sigma_deg, "\t",
				$object_name, "\t",
				$obscode, "\t",
				$is_synthetic, "\t",
				$status, "\n";

			$i++;
		}   # for
	};

    warn $@ if $@;
    return $i;      # return num dets added
}

__END__

=head1 NAME

detCreateRawFalse - Program to create a file of false detections.

=head1 SYNOPSIS

detCreateRawFalse --fields FIELDS [options] [--help]

  --fields FIELDS : file containing fields
  --instance INSTANCE : instance name
  --prefix PREFIX : write detection descriptions with prefix PREFIX
  --prefix SUFFIX : write detection descriptions with suffix SUFFIX
  --help : show man page
  --nofalse : Don't add false detections

=head1 DESCRIPTION

insertFPDetections is currently the primary method of inserting synthetic
and false detections from the MOPS SSM into the detections DC.

=head1 ENVIRONMENT VARIABLES

=over 4

=item MOPS_DBINSTANCE

Simulation instance.

=back

=cut
