Index: trunk/ippScripts/scripts/magic_tree.pl
===================================================================
--- trunk/ippScripts/scripts/magic_tree.pl	(revision 15681)
+++ trunk/ippScripts/scripts/magic_tree.pl	(revision 15681)
@@ -0,0 +1,311 @@
+#!/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 PS::IPP::Operations qw( skycell_file
+			    generate_skycells
+			    );
+
+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, $workdir, $save_temps, $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
+	   'workdir=s'     => \$workdir,    # Working directory, for output files
+	   'save-temps'    => \$save_temps, # Save temporary files?
+	   '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 --workdir",
+	   -exitval => 3) unless
+    defined $magic_id and
+    defined $tess_id and
+    defined $ra0 and
+    defined $dec0 and
+    defined $camera and
+    defined $workdir;
+
+$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 => 1);
+    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};
+    }
+}
+
+generate_skycells($ipprc, $tess_id, \@skycells, $workdir) or
+    &my_die("Unable to generate skycells", $magic_id, $PS_EXIT_PROG_ERROR);
+
+### For each skycell, project centre of skycell onto tangent plane of boresight
+my @fields;
+foreach my $skycell_id ( @skycells ) {
+    my $skyfile = skycell_file($tess_id, $skycell_id, $workdir) or 
+	&my_die("Unable to get name of skycell file", $magic_id, $PS_EXIT_PROG_ERROR);
+    &my_die("Unable to find skycell file", $magic_id, $PS_EXIT_SYS_ERROR) unless
+	$skyfile->exists( $skyfile );
+
+    my ($header, $status) = Astro::FITS::CFITSIO::fits_read_header( $skyfile );
+    &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
+};
+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 => 1);
+	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 ($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 node into upper and lower nodes
+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 $index = (not defined $node->{index} or $node->{index} eq 'xi') ? 'eta' : 'xi'; # Property to sort
+    my @sorted = sort { $$a{$index} <=> $$b{$index} } @$contents; # Sorted list
+    my $median = int(scalar @sorted / 2); # Median point of list
+    my @top = splice(@sorted, $median);	# Top of the sorted list
+
+    my $upper = {		# Upper node
+	contents => \@top,
+	index => $index,
+	position => $node->{position} . 'U',
+    };
+    my $lower = { 		# Lower node
+	contents => \@sorted,
+	index => $index,
+	position => $node->{position} . 'L',
+    };
+
+    $node->{upper} = $upper;
+    $node->{lower} = $lower;
+    $node->{contents} = undef;
+
+    push @$tasks, $upper if scalar @top > MAX_FIELDS;
+    push @$tasks, $lower if scalar @sorted > MAX_FIELDS;
+
+    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 {
+	$output .= "$position\t\tSTR\t${position}L\n";
+	$output .= "$position\t\tSTR\t${position}U\n";
+	$output .= print_node($node->{lower});
+	$output .= print_node($node->{upper});
+    }
+
+    return $output;
+}
+
+END {
+    my $status = $?;
+    system("sync") == 0
+        or die "failed to execute sync: $!" ;
+    $? = $status;
+}
+
+__END__