Index: trunk/ippScripts/scripts/ipp_mops_translate.pl
===================================================================
--- trunk/ippScripts/scripts/ipp_mops_translate.pl	(revision 17323)
+++ trunk/ippScripts/scripts/ipp_mops_translate.pl	(revision 17323)
@@ -0,0 +1,269 @@
+#!/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;