#!/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

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

GetOptions(
	   'input=s'   => \$input,
	   'extname=s' => \$extname,
	   'skycell=s' => \$skycell,
	   'output=s'  => \$output,
) 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.NAME'     => { 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' },
    'MJD-OBS'      => { name => 'MJD-OBS',  type => TDOUBLE, comment => 'Time of exposure' },
	    };


# 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 PHU keywords
$outFits->create_img( 16, 0, undef, $status );
check_fitsio( $status );
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 );
}

# Write the table
$outFits->create_tbl( BINARY_TBL(), $numRows, scalar @colNames, \@colNames, \@colTypes, undef, EXTNAME,
		      $status );
check_fitsio( $status );
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));
    
    return (rad2deg($ra), rad2deg($dec));
}


__END__







        fits_create_tbl(fits->fd,
                        BINARY_TBL,
                        table->n, // number of rows in table
                        numColumns, // number of columns in table
                        (char**)columnNames->data, // names of the columns
                        (char**)columnTypes->data, // format of the columns
                        NULL, // physical unit of columns
                        NULL, // skip extension name: we set the by hand below
                        &status);


my $status;			# Status of FITSIO calls
my $fits = Astro::FITS::CFITSIO::open_file( $filename, READONLY, $status ); # FITS file handle
check_fitsio($status);
$fits->movnam_hdu(BINARY_TBL, $extname, 0, $status) and check_fitsio($status);
my $numRows;			# Number of rows in table
$fits->get_num_rows($numRows, $status) and check_fitsio($status);

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

my ($xType, $yType);		# Types for x and y
$fits->get_coltype($xCol, $xType, undef, undef, $status) and check_fitsio($status);
$fits->get_coltype($yCol, $yType, undef, undef, $status) and check_fitsio($status);

my ($x, $y);			# Coordinates read from table
$fits->read_col($xType, $xCol, 1, 1, $numRows, 0, $x, undef, $status) and check_fitsio($status);
$fits->read_col($yType, $yCol, 1, 1, $numRows, 0, $y, undef, $status) and check_fitsio($status);

my ($coordFile, $coordName) = tempfile( "/tmp/ds9_cmf_regions.XXXX", UNLINK => 1 );
for (my $i = 0; $i < $numRows; $i++) {
    print $coordFile "image; circle(" . ($$x[$i] + 1) . ',' . ($$y[$i] + 1) .
	",$radius \# color = $colour\n";
}
close $coordFile;
