#!/usr/bin/env perl

use warnings;
use strict;

use Astro::FITS::CFITSIO qw( :constants );
Astro::FITS::CFITSIO::PerlyUnpacking(1);
use Getopt::Long qw( GetOptions :config auto_help auto_version gnu_getopt );
use Pod::Usage qw( pod2usage );
use Math::Trig;
use Data::Dumper;

use constant EXTNAME => 'MOPS_TRANSIENT_DETECTIONS'; # Extension name for output table
use constant POSITION_ERROR => 0.2 / 3600; # Error in positions, degrees
use constant ZERO_POINT => 25;  # Magnitude zero point
use constant MAG_ERROR => 0.1;  # Error in magnitudes
use constant OBSERVATORY_CODE => 566; # IAU Observatory Code
use constant FAKE_LIMITING_MAG => 23.0; # Fake limiting magnitude to report
use constant FAKE_DETECTION_EFFICIENCY => 0.0; # Fake detection efficiency to report

my ( $input,                    # Name of input file with IPP photometry
     $extname,                  # Name of extension containing photometry
     $skycell,                  # Skycell file with WCS
     $output,                   # Name of output file
     $save_temps,               # Save temporary files?
     );

GetOptions(
           'input=s'    => \$input,
           'extname=s'  => \$extname,
           'skycell=s'  => \$skycell,
           'output=s'   => \$output,
           'save-temps' => \$save_temps,
) or pod2usage( 2 );

pod2usage( -msg => "Unknown option: @ARGV", -exitval => 2 ) if @ARGV;
pod2usage( -msg => "Required options: --input --extname --skycell --output",
           -exitval => 3)
    unless defined $input
    and defined $extname
    and defined $skycell
    and defined $output;

# Specification of columns to write
my $columns = [ { name => 'RA_DEG',   type => 'D' }, # Right ascension
                { name => 'DEC_DEG',  type => 'D' }, # Declination
                { name => 'RA_SIG',   type => 'D' }, # Error in right ascension
                { name => 'DEC_SIG',  type => 'D' }, # Error in declination
                { name => 'FLUX',     type => 'D' }, # Flux
                { name => 'FLUX_SIG', type => 'D' }, # Error in flux
                { name => 'ANG',      type => 'D' }, # Angle
                { name => 'ANG_SIG',  type => 'D' }, # Error in angle
                { name => 'LEN',      type => 'D' }, # Length
                { name => 'LEN_SIG',  type => 'D' }, # Error in length
                ];

# Header translation table
my $headers = {
    'FPA.RA'       => { name => 'RA',       type => TSTRING, comment => 'Right ascension of boresight'},
    'FPA.DEC'      => { name => 'DEC',      type => TSTRING, comment => 'Declination of boresight' },
    'FPA.OBS'      => { name => 'FPA_ID',   type => TSTRING, comment => 'Exposure identifier' },
    'FPA.FILTER'   => { name => 'FILTER',   type => TSTRING, comment => 'Filter name' },
    'EXPTIME'      => { name => 'EXPTIME',  type => TDOUBLE, comment => 'Exposure time' },
    'FPA.POSANGLE' => { name => 'ROTANGLE', type => TDOUBLE, comment => 'Position angle' },
    'FPA.ALT'      => { name => 'TEL_ALT',  type => TDOUBLE, comment => 'Telescope altitude' },
    'FPA.AZ'       => { name => 'TEL_AZ',   type => TDOUBLE, comment => 'Telescope azimuth' },
            };


# Read the skycell header for the WCS
my ($wcsheader, $status) = Astro::FITS::CFITSIO::fits_read_header( $skycell );

# Get the useful header keywords
my $naxis1 = $$wcsheader{'NAXIS1'} or die("Can't find NAXIS1");
my $naxis2 = $$wcsheader{'NAXIS2'} or die("Can't find NAXIS2");
my $ctype1 = $$wcsheader{'CTYPE1'} or die("Can't find CTYPE1");
my $ctype2 = $$wcsheader{'CTYPE2'} or die("Can't find CTYPE2");
my $cdelt1 = $$wcsheader{'CDELT1'} or die("Can't find CDELT1");
my $cdelt2 = $$wcsheader{'CDELT2'} or die("Can't find CDELT2");
my $crval1 = deg2rad($$wcsheader{'CRVAL1'}) or die("Can't find CRVAL1");
my $crval2 = deg2rad($$wcsheader{'CRVAL2'}) or die("Can't find CRVAL2");
my $crpix1 = $$wcsheader{'CRPIX1'} or die("Can't find CRPIX1");
my $crpix2 = $$wcsheader{'CRPIX2'} or die("Can't find CRPIX2");
my $pc11 = $$wcsheader{'PC001001'} or die("Can't find PC001001");
my $pc12 = $$wcsheader{'PC001002'} or die("Can't find PC001002");
my $pc21 = $$wcsheader{'PC002001'} or die("Can't find PC002001");
my $pc22 = $$wcsheader{'PC002002'} or die("Can't find PC002002");
my $crota1 = $$wcsheader{'CROTA1'};
my $crota2 = $$wcsheader{'CROTA2'};

die("Unexpected projection: $ctype1 and $ctype2.") unless $ctype1 =~ /^\'RA---TAN\s*\'$/ and
    $ctype2 =~ /^\'DEC--TAN\s*\'$/;
die("Can't determine size of skycell ($naxis1,$naxis2)") unless $naxis1 > 0 and $naxis2 > 0;
die("Can't determine scale of skycell ($cdelt1,$cdelt2)") if not defined $cdelt1 or $cdelt1 == 0 or
    not defined $cdelt2 or $cdelt2 == 0;
die("We don't know how to handle rotations ($crota1,$crota2)") if defined $crota1 or defined $crota2;


# Read the input table
my $inFits = Astro::FITS::CFITSIO::open_file( $input, READONLY, $status ); # FITS file handle
my $inHeader = $inFits->read_header(); # Header for input
check_fitsio($status);
$inFits->movnam_hdu(BINARY_TBL, $extname, 0, $status) and check_fitsio($status);
my $numRows;                    # Number of rows in table
$inFits->get_num_rows($numRows, $status) and check_fitsio($status);

my ($xCol, $yCol, $magCol);             # Column numbers for x and y
$inFits->get_colnum(0, 'X_PSF', $xCol, $status) and check_fitsio($status);
$inFits->get_colnum(0, 'Y_PSF', $yCol, $status) and check_fitsio($status);
$inFits->get_colnum(0, 'PSF_INST_MAG', $magCol, $status) and check_fitsio($status);

my ($xType, $yType, $magType);  # Column types
$inFits->get_coltype($xCol, $xType, undef, undef, $status) and check_fitsio($status);
$inFits->get_coltype($yCol, $yType, undef, undef, $status) and check_fitsio($status);
$inFits->get_coltype($magCol, $magType, undef, undef, $status) and check_fitsio($status);

my ($x, $y, $mag);              # Sources arrays
$inFits->read_col($xType, $xCol, 1, 1, $numRows, 0, $x, undef, $status) and check_fitsio($status);
$inFits->read_col($yType, $yCol, 1, 1, $numRows, 0, $y, undef, $status) and check_fitsio($status);
$inFits->read_col($magType, $magCol, 1, 1, $numRows, 0, $mag, undef, $status) and check_fitsio($status);

$inFits->close_file( $status );

# Parse the list of columns
my @colNames;                   # Names of columns
my @colTypes;                   # Types of columns
my %colData;                    # Data for each column
foreach my $colSpec ( @$columns) {
    push @colNames, $colSpec->{name};
    push @colTypes, $colSpec->{type};
    $colData{$colSpec->{name}} = [];
}


# Convert the input data into the output formats
for (my $i = 0; $i < $numRows; $i++) {
    my ($ra, $dec) = pixels_to_sky($$x[$i], $$y[$i]); # Sky coordinates

    push @{$colData{'RA_DEG'}}, $ra;
    push @{$colData{'DEC_DEG'}}, $dec;
    push @{$colData{'RA_SIG'}}, POSITION_ERROR;
    push @{$colData{'DEC_SIG'}}, POSITION_ERROR;
    push @{$colData{'FLUX'}}, $$mag[$i] + ZERO_POINT;
    push @{$colData{'FLUX_SIG'}}, MAG_ERROR;
    push @{$colData{'ANG'}}, 0.0;
    push @{$colData{'ANG_SIG'}}, 0.0;
    push @{$colData{'LEN'}}, 0.0;
    push @{$colData{'LEN_SIG'}}, 0.0;
}


# Write the output
unlink "$output" if -e "$output";
my $outFits = Astro::FITS::CFITSIO::create_file( $output, $status ); # Output file handle
check_fitsio( $status );

# Write the table
$outFits->create_tbl( BINARY_TBL(), $numRows, scalar @colNames, \@colNames, \@colTypes, undef, EXTNAME,
                      $status );
check_fitsio( $status );

# Write the header keywords
foreach my $keyword ( keys %$headers ) {
    my $value = $inHeader->{$keyword}; # Header keyword value
    unless (defined $value) {
        print "Can't find header keyword $keyword\n";
        next;
    }
    $value =~ s/\'//g;
    my $name = $headers->{$keyword}->{name}; # New name
    my $type = $headers->{$keyword}->{type}; # Type
    my $comment = $headers->{$keyword}->{comment}; # Comment
    $outFits->write_key( $type, $name, $value, $comment, $status );
    check_fitsio( $status );
}

# Adjust the time from start to mid-point
{
    my $mjd = $inHeader->{'MJD-OBS'} or die "Can't find MJD.\n"; # Modified Julian Date
    my $exptime = $inHeader->{'EXPTIME'} or die "Can't find EXPTIME.\n"; # Exposure time, seconds
    $mjd += $exptime / 2.0 / 3600 / 24;     # start --> mid-point
    $outFits->write_key( TDOUBLE, 'MJD-OBS', $mjd, 'Time of exposure mid-point', $status);
    check_fitsio( $status );
}

# Stellar PSF
my $fwhm = 0.5 * 3600 * ($inHeader->{'FWHM_X'} * $cdelt1 + $inHeader->{'FWHM_Y'} * $cdelt2); # FWHM for star
$outFits->write_key( TSTRING, 'STARPSF', $fwhm, 'Stellar PSF (arcsec)', $status );
check_fitsio( $status );

# Observatory code
$outFits->write_key( TINT, 'OBSCODE', OBSERVATORY_CODE, 'IAU observatory code', $status );

# Limiting magnitude
$outFits->write_key( TINT, 'LIMITMAG', FAKE_LIMITING_MAG, 'Limiting magnitude (FAKE)', $status );

# Detection efficiency
$outFits->write_key( TINT, 'DE1', FAKE_DETECTION_EFFICIENCY, 'Detection efficiency (FAKE)', $status );
$outFits->write_key( TINT, 'DE2', FAKE_DETECTION_EFFICIENCY, 'Detection efficiency (FAKE)', $status );
$outFits->write_key( TINT, 'DE3', FAKE_DETECTION_EFFICIENCY, 'Detection efficiency (FAKE)', $status );
$outFits->write_key( TINT, 'DE4', FAKE_DETECTION_EFFICIENCY, 'Detection efficiency (FAKE)', $status );
$outFits->write_key( TINT, 'DE5', FAKE_DETECTION_EFFICIENCY, 'Detection efficiency (FAKE)', $status );
$outFits->write_key( TINT, 'DE6', FAKE_DETECTION_EFFICIENCY, 'Detection efficiency (FAKE)', $status );
$outFits->write_key( TINT, 'DE7', FAKE_DETECTION_EFFICIENCY, 'Detection efficiency (FAKE)', $status );
$outFits->write_key( TINT, 'DE8', FAKE_DETECTION_EFFICIENCY, 'Detection efficiency (FAKE)', $status );
$outFits->write_key( TINT, 'DE9', FAKE_DETECTION_EFFICIENCY, 'Detection efficiency (FAKE)', $status );
$outFits->write_key( TINT, 'DE10', FAKE_DETECTION_EFFICIENCY, 'Detection efficiency (FAKE)', $status );

# Write the data
for (my $i = 0; $i < scalar @colNames; $i++) {
    my $colName = $colNames[$i];# Column name
    $outFits->write_col( TDOUBLE, $i + 1, 1, 1, $numRows, $colData{$colName}, $status );
    check_fitsio( $status );
}

$outFits->close_file( $status );

### Pau.



# From Astro::FITS::CFITSIO demo
sub check_fitsio
{
    my $status = shift;         # Status of FITSIO calls

    if ($status != 0) {
        my $msg;                # Message to output
        Astro::FITS::CFITSIO::fits_get_errstatus( $status , $msg );
        die "CFITSIO error: $msg\n";
    }
}

# Convert pixel coordinates to sky coordinates
sub pixels_to_sky
{
    my ($x, $y) = @_;           # Coordinates

    # Pixel coordinate relative to reference
    $x -= $crpix1;
    $y -= $crpix2;

    # Coordinates on tangent plane
    my $xi = $pc11 * ($x) + $pc12 * ($y);
    my $eta = $pc21 * ($x) + $pc22 * ($y);
    $xi *= $cdelt1;
    $eta *= $cdelt2;

    # Coordinates on rotated celestial sphere
    my $phi = atan2($eta,$xi) + pi/2;
    my $theta = atan(180 / pi / sqrt($xi**2 + $eta**2));

    # Coordinates on celestial sphere
    my $ra = $crval1 + atan2(cos($theta) * sin($phi),
                             sin($theta) * cos($crval2) + cos($theta) * sin($crval2) * cos($phi));
    my $dec = asin(sin($theta) * sin($crval2) - cos($theta) * cos($crval2) * cos($phi));

    $ra = rad2deg($ra);
    $dec = rad2deg($dec);

    $ra += 360 if $ra < 0;
    $dec += 360 if $dec < 0;
    $ra -= 360 if $ra >= 360;
    $dec -= 360 if $dec >= 360;

    return ($ra, $dec);
}


__END__

