#!/usr/bin/env perl

use Carp;
use warnings;
use strict;

## report the program and machine
use Sys::Hostname;
my $host = hostname();
print "\n\n";
print "Starting script $0 on $host\n\n";

use vars qw( $VERSION );
$VERSION = '0.01';

use IPC::Cmd 0.36 qw( can_run run );
use PS::IPP::Metadata::Config;
use PS::IPP::Metadata::List qw( parse_md_list );

use Astro::FITS::CFITSIO qw( :constants );
Astro::FITS::CFITSIO::PerlyUnpacking(1);

use Math::Trig;
use File::Temp qw( tempfile );

use PS::IPP::Config qw($PS_EXIT_SUCCESS
		       $PS_EXIT_UNKNOWN_ERROR
		       $PS_EXIT_SYS_ERROR
		       $PS_EXIT_CONFIG_ERROR
		       $PS_EXIT_PROG_ERROR
		       $PS_EXIT_DATA_ERROR
		       $PS_EXIT_TIMEOUT_ERROR
		       caturi
		       );
my $ipprc = PS::IPP::Config->new(); # IPP configuration

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

use constant MAX_FIELDS => 4;	# Maximum number of fields to be in a node

# Parse the command-line arguments
my ($magic_id, $tess_id, $camera, $ra0, $dec0, $dbname, $outroot, $save_temps, $verbose, $no_update, $no_op);
GetOptions(
	   'magic_id=s'    => \$magic_id,   # Magic identifier
	   'tess_id=s'     => \$tess_id,    # Tessellation identifier
	   'camera=s'      => \$camera,	    # Camera name
	   'ra=f'          => \$ra0,        # Boresight right ascension, radians
	   'dec=f'         => \$dec0,       # Boresight declination, radians
	   'dbname=s'      => \$dbname,     # Database name
	   'outroot=s'     => \$outroot,    # Output root name
	   'save-temps'    => \$save_temps, # Save temporary files?
	   'verbose'       => \$verbose,    # Print stuff?
	   'no-update'     => \$no_update,  # Don't update the database?
	   'no-op'         => \$no_op,	    # Don't do any operations?
	   ) or pod2usage( 2 );

pod2usage( -msg => "Unknown option: @ARGV", -exitval => 2 ) if @ARGV;
pod2usage( -msg => "Required options: --magic_id --camera --outroot",
	   -exitval => 3) unless
    defined $magic_id and
    defined $tess_id and
    defined $ra0 and
    defined $dec0 and
    defined $camera and
    defined $outroot;

$ipprc->define_camera($camera);


# Look for programs we need
my $missing_tools;
my $magictool = can_run('magictool') or (warn "Can't find magictool" and $missing_tools = 1);
my $ppImage = can_run('ppImage') or (warn "Can't find ppImage" and $missing_tools = 1);
if ($missing_tools) { 
    warn("Can't find required tools.");
    exit($PS_EXIT_CONFIG_ERROR); 
}

my $mdcParser = PS::IPP::Metadata::Config->new;	# Parser for metadata config files

### Get a list of skycells
my @skycells;			# List of skycells
{
    my $command = "$magictool -inputskyfile -magic_id $magic_id"; # Command to run
    $command .= " -dbname $dbname" if defined $dbname;
    my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
	run(command => $command, verbose => $verbose);
    unless ($success) {
	$error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR);
	&my_die("Unable to perform magictool -inputfile: $error_code", $magic_id, $error_code);
    }

    my $metadata = $mdcParser->parse(join "", @$stdout_buf) or
	&my_die("Unable to parse metadata config doc", $magic_id, $PS_EXIT_PROG_ERROR);

    my $inputs = parse_md_list($metadata) or
	&my_die("Unable to parse metadata list", $magic_id, $PS_EXIT_PROG_ERROR);

    foreach my $input ( @$inputs ) {
	push @skycells, $input->{skycell_id};
    }
}

### For each skycell, project centre of skycell onto tangent plane of boresight
my @fields;
foreach my $skycell_id ( @skycells ) {
    my $skyfile = $ipprc->filename("SKYCELL.TEMPLATE", $outroot, $skycell_id );
    $ipprc->skycell_file($tess_id, $skycell_id, $skyFile, $verbose) or &my_die("Unable to generate skycells $skycell_id", $magic_id, $PS_EXIT_PROG_ERROR);
    my $skyfileResolved = $ipprc->file_resolve( $skyfile );
    my ($header, $status) = Astro::FITS::CFITSIO::fits_read_header( $skyfileResolved );
    &my_die("Unable to read skycell header: $status", $magic_id, $PS_EXIT_SYS_ERROR) if $status;
    
    # Get the useful header keywords
    my $naxis1 = $$header{'NAXIS1'} or &my_die("Can't find NAXIS1", $magic_id, $PS_EXIT_SYS_ERROR);
    my $naxis2 = $$header{'NAXIS2'} or &my_die("Can't find NAXIS2", $magic_id, $PS_EXIT_SYS_ERROR);
    my $ctype1 = $$header{'CTYPE1'} or &my_die("Can't find CTYPE1", $magic_id, $PS_EXIT_SYS_ERROR);
    my $ctype2 = $$header{'CTYPE2'} or &my_die("Can't find CTYPE2", $magic_id, $PS_EXIT_SYS_ERROR);
    my $cdelt1 = $$header{'CDELT1'} or &my_die("Can't find CDELT1", $magic_id, $PS_EXIT_SYS_ERROR);
    my $cdelt2 = $$header{'CDELT2'} or &my_die("Can't find CDELT2", $magic_id, $PS_EXIT_SYS_ERROR);
    my $crval1 = $$header{'CRVAL1'} or &my_die("Can't find CRVAL1", $magic_id, $PS_EXIT_SYS_ERROR);
    my $crval2 = $$header{'CRVAL2'} or &my_die("Can't find CRVAL2", $magic_id, $PS_EXIT_SYS_ERROR);
    my $crpix1 = $$header{'CRPIX1'} or &my_die("Can't find CRPIX1", $magic_id, $PS_EXIT_SYS_ERROR);
    my $crpix2 = $$header{'CRPIX2'} or &my_die("Can't find CRPIX2", $magic_id, $PS_EXIT_SYS_ERROR);
    my $pc11 = $$header{'PC001001'} or &my_die("Can't find PC001001", $magic_id, $PS_EXIT_SYS_ERROR);
    my $pc12 = $$header{'PC001002'} or &my_die("Can't find PC001002", $magic_id, $PS_EXIT_SYS_ERROR);
    my $pc21 = $$header{'PC002001'} or &my_die("Can't find PC002001", $magic_id, $PS_EXIT_SYS_ERROR);
    my $pc22 = $$header{'PC002002'} or &my_die("Can't find PC002002", $magic_id, $PS_EXIT_SYS_ERROR);
    my $crota1 = $$header{'CROTA1'};
    my $crota2 = $$header{'CROTA2'};

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

    # Relative coordinates of centre of the field
    my $x = $naxis1 - $crpix1;
    my $y = $naxis2 - $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
    $crval1 = deg2rad($crval1);
    $crval2 = deg2rad($crval2);
    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));

    # Rotate to boresight
    my $phi_new = atan2(cos($dec) * sin($ra - $ra0), 
			sin($dec) * cos($dec0) + cos($dec) * sin($dec0) * cos($ra - $ra0));
    my $theta_new = asin(sin($dec) * sin($dec0) - cos($dec) * cos($dec0) * cos($ra - $ra0));

    # Project
    my $rad = 180 / pi * cot($theta_new);
    my $xi_new = $rad * sin($phi);
    my $eta_new = - $rad * cos($phi);

    my $field = { id => $skycell_id,
		  xi => $xi_new,
		  eta => $eta_new,
	      };

    push @fields, $field;
}

### Subdivide list of positions into kd-tree
my $root = {			# Root node of tree
    contents => \@fields,	# Contents of node
    position => 'root',		# Position in tree
    children => {},		# Children of node
};
my @tasks = ( $root );
while (scalar @tasks > 0) {
    my $node = shift @tasks;
    divide_node($node, \@tasks);
}

### Format tree for magictool
my $mdcTree = print_node($root); # The tree in MDC format
my ($treeFile, $treeName) = tempfile( "magictree.${magic_id}.XXXX", UNLINK => !$save_temps );
print $treeFile, $mdcTree;
close $treeFile;

### Input tree into database
{
    my $command = "$magictool -inputtree";
    $command   .= " -magic_id $magic_id";
    $command   .= " -dep_file $treeName";
    $command   .= " -dbname $dbname" if defined $dbname;

    # Add the processed file to the database
    unless ($no_update) {
	my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
	    run(command => $command, verbose => $verbose);
	unless ($success) {
	    $error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR);
	    warn("Unable to perform magictool -inputtree: $error_code");
	    exit($error_code);
	}
    } else {
	print "Skipping command: $command\n";
    }
}

### Pau.

sub my_die
{
    my $msg = shift;		# Warning message on die
    my $magic_id = shift;	# Magic identifier
    my $exit_code = shift;	# Exit code to add

    carp($msg);
    if (defined $magic_id and not $no_update) {
	my $command = "$magictool -inputtree";
	$command .= " -magic_id $magic_id";
	$command .= " -code $exit_code";
	$command .= " -dbname $dbname" if defined $dbname;
	system($command);
    }
    exit $exit_code;
}

# Divide a list into two, returning the lower and upper parts
sub divide_list
{
    my $list = shift;		# List to divide
    my $index = shift;		# Name of index for sorting

    my @sorted = sort { $$a{$index} <=> $$b{$index} } @$list; # Sorted list
    my $median = int(scalar @sorted / 2); # Median point of list
    my @upper = splice(@sorted, $median); # Upper part of the sorted list

    return (\@sorted, \@upper);
}

# Create a new node, add it to the parent, and add it to the task list if required
sub new_node
{
    my $parent = shift;		# The parent node
    my $contents = shift;	# Contents of the new node
    my $position = shift;	# Position description
    my $tasks = shift;		# Tasks to do

    my $node = {
	contents => $contents,
	position => $parent->{position} . '_' . $position,
	children => {},
    };

    $parent->{children}->{$position} = $node;

    push @$tasks, $node if scalar @$contents > 4;

    return $node;
}

# Divide a node
sub divide_node
{
    my $node = shift;		# Node to divide
    my $tasks = shift;		# Tasks to do

    my $position = $node->{position};

    my $contents = $node->{contents} or die "Can't find contents of node."; # Contents of node

    my ($lower, $upper) = divide_list($contents, 'xi');

    if (scalar @$lower > 4) {
	my ($ll, $lr) = divide_list($lower, 'eta');
	new_node($node, $ll, 'll', $tasks);
	new_node($node, $lr, 'lr', $tasks);
    } else {
	new_node($node, $lower, 'L', $tasks);
    }

    if (scalar @$upper > 4) {
	my ($ul, $ur) = divide_list($upper, 'eta');
	new_node($node, $ul, 'ul', $tasks);
	new_node($node, $ur, 'ur', $tasks);
    } else {
	new_node($node, $upper, 'U', $tasks);
    }

    $node->{contents} = undef;

    return $node;
}

# Print the contents of a node
sub print_node
{
    my $node = shift;		# Node to print

    my $position = $node->{position}; # Position of node

    my $output = "$position\t\tMULTI\n"; # Output text

    if (defined $node->{contents}) {
	foreach my $field ( @{$node->{contents}} ) {
	    my $skycell_id = $field->{id};	# Skycell name
	    $output .= "$position\t\tSTR\t$skycell_id\n";
	    $output .= "$skycell_id\t\tSTR\tNULL\t\# $field->{xi},$field->{eta}\n";
	}
    } else {
	foreach my $div ( keys %{$node->{children}} ) {
	    $output .= "$position\t\tSTR\t${position}_$div\n";
	    $output .= print_node($node->{children}->{$div});
	}
    }

    return $output;
}

END {
    my $status = $?;
    system("sync") == 0
        or die "failed to execute sync: $!" ;
    $? = $status;
}

__END__
