#!/usr/bin/env perl

# given a table file that contains rows from a P2 detections table (either from the
# results of a PSPS query # or ipp2psps P2 batch file)
# build a set of delta columns which give the difference between the p2 columns and
# the expected value.
# The list of columns and the formulas are given in this file by way of a perl here document.
# optionally check that the differences are small

use strict;
use warnings;

use Getopt::Long qw( GetOptions :config auto_help auto_version gnu_getopt );
use Pod::Usage qw( pod2usage );
use File::Temp qw( tempdir  tempfile );


my $verbose;
my $checkdeltas;
my $save_temps;

GetOptions(
    "check"             =>  \$checkdeltas,
    "verbose|v"         =>  \$verbose,
    "save-temps"        =>  \$save_temps,
) or pod2usage( 2 );


my $in = shift;
my $out = shift;

die "usage $0 <input matched detections table> <output> [--check]\n" unless $in and $out;

my $jardir = $ENV{JARDIR};
$jardir =  '/home/panstarrs/bills/jars' unless $jardir;

my $stiltsJar =  "$jardir/stilts.jar";
my $stilts = "java -jar $stiltsJar";

my $p2columns = getp2columns();

my @lines = split "\n", $p2columns;

my $cmd="'addcol -after zp zpFactor 3630.78*pow(10,-0.4*zp);";

# foreach my $line (<C>) {
my @deltas;
my @columns;
foreach my $line (@lines) {
    next if $line =~ '^#';
    next if !$line;
    chomp $line;
    my ($data, $comment) = split '#', $line;
    next if !$data;
    my ($p2, $ipp, $conversion) = split ' ', $data;

    next if !$ipp;
    $conversion = "" unless $conversion;
    my $equation = "$ipp$conversion";
    $equation = "$p2-$equation";
    my $dcol = "d_$p2";
    $cmd .=  "addcol -after $p2 $dcol $equation;";
    push @deltas, $dcol;
    push @columns, $p2;
}
$cmd .= "'";

my $command = "$stilts tpipe in=$in out=$out cmd=$cmd";

print "$command\n" if $verbose;

my $rc = system $command;
if ($rc) {
    exit $rc >> 8;
}

if ($checkdeltas) {
    my ($manafile, $mananame) = tempfile( "/var/tmp/manacmds.XXXX", UNLINK => !$save_temps);
    print "$mananame\n" if $save_temps;
    print $manafile (<<HERE);
\$limit = 1e-9
output /dev/null
macro checkdelta
    \$dcol = \$1
    \$col = \$2
    \$thislimit = \$limit
    set goodmask = \$col != -999
    subset good = \$dcol where goodmask == 1
    \$length = good[]
    for i 0 10
        subset d = good where abs(good) > \$thislimit
        \$noutliers = d[]
        if (\$noutliers > 0) 
            echo \$dcol has \$noutliers out of \$length with difference > \$thislimit
            \$thislimit = \$thislimit * 10
        else
            return
        end
    end
end
HERE

    print $manafile "\nmacro dotests\n";
    my $n = scalar @deltas;
    for (my $i=0; $i < $n; $i++) {
        my $col = $columns[$i];
        my $dcol = $deltas[$i];
        print $manafile "    checkdelta $dcol $col\n";
    }
    print $manafile "end\n";

    print $manafile "datafile $out\n";
    print $manafile "read -fits Joined";
    foreach my $col (@deltas) {
        print $manafile " $col";
    }
    print $manafile "\n";
    print $manafile "read -fits Joined";
    foreach my $col (@columns) {
        print $manafile " $col";
    }
    print $manafile "\n";


    print $manafile "output stdout\n";
    print $manafile "dotests\n";
    close $manafile;

    system "echo input $mananame | mana";
}
exit 0;


sub getp2columns {

my $list = <<END_OF_LIST;
# Note zpFactor = 3630.78*pow(10, -0.4*zp)

objID
uniquePspsP2id
detectID
ippObjID
ippDetectID     IPP_IDET
filterID
surveyID
imageID
randomDetID
dvoRegionID
obsTime
xPos            X_PSF
yPos            Y_PSF
xPosErr         X_PSF_SIG
yPosErr         Y_PSF_SIG
pltScale        PLTSCALE
posAngle        POSANGLE
ra
dec
raErr           X_PSF_SIG               *abs(PLTSCALE)
decErr          Y_PSF_SIG               *abs(PLTSCALE)
extNSigma       EXT_NSIGMA
zp
telluricExt
expTime
airMass
psfFlux         PSF_INST_FLUX 		*zpFactor/expTime/PSF_NDOF
psfFluxErr      PSF_INST_FLUX_SIG	*zpFactor/expTime
psfMajorFWHM    PSF_FWHM_MAJ            *abs(PLTSCALE)
psfMinorFWHM    PSF_FWHM_MIN            *abs(PLTSCALE)
psfTheta        PSF_THETA
psfCore         PSF_CORE
psfQf           PSF_QF
psfQfPerfect    PSF_QF_PERFECT
psfChiSq        PSF_CHISQ
psfLikelihood                 # actually psfLikelihood(EXT_NSIGMA). We assume the formula is ok.
momentXX        MOMENTS_XX		*PLTSCALE*PLTSCALE
momentXY        MOMENTS_XY		*PLTSCALE*PLTSCALE
momentYY        MOMENTS_YY		*PLTSCALE*PLTSCALE
momentR1        MOMENTS_R1              *abs(PLTSCALE)
momentRH        MOMENTS_RH              *sqrt(abs(PLTSCALE))
momentM3C       MOMENTS_M3C		*PLTSCALE*PLTSCALE
momentM3S       MOMENTS_M3S		*PLTSCALE*PLTSCALE
momentM4C       MOMENTS_M4C		*PLTSCALE*PLTSCALE
momentM4S       MOMENTS_M4S		*PLTSCALE*PLTSCALE
apFlux          AP_FLUX		        *zpFactor/expTime
apFluxErr       AP_FLUX_SIG		*zpFactor/expTime
apFillF         AP_NPIX                 /(3.14159265359*AP_MAG_RADIUS*AP_MAG_RADIUS)
apRadius        AP_MAG_RADIUS           *abs(PLTSCALE)
kronFlux        KRON_FLUX		*zpFactor/expTime
kronFluxErr     KRON_FLUX_ERR		*zpFactor/expTime
kronRad         MOMENTS_R1              *2.5*abs(PLTSCALE)
sky             SKY                     /PLTSCALE/PLTSCALE*zpFactor/expTime
skyErr          SKY_SIGMA               /PLTSCALE/PLTSCALE*zpFactor/expTime
infoFlag        FLAGS
infoFlag2       FLAGS2
infoFlag3
processingVersion
END_OF_LIST
    return $list;
}
