#!/usr/bin/perl


# This script reads in a FITS file with extensions and computes an offset
# list (suitable for use by IRAF's imcombine) from the WCS parameters.
# It will also accept a magnification factor on the command line if you
# want to magnify the offsets (useful if you've magnified the images).
#
# C. Engelbracht 2002-11-15


# set parameters

$pi = 4 * atan2 (1, 1);


# do some Astro::FITS::CFITSIO setup

use Astro::FITS::CFITSIO;
$readonly = Astro::FITS::CFITSIO::READONLY();


# get input filename

if ($ARGV[0]) {$file = $ARGV[0];}
else {
   print "Please enter name of FITS file: ";
   chop ($file = <STDIN>);
}


# get magnification factor

if (defined ($ARGV[1])) {$mag = $ARGV[1];}
else {$mag = 1;}


# open input FITS file

Astro::FITS::CFITSIO::fits_open_file ($fptr, $file, $readonly, $status);
if ($status) {die "Error opening $file:  status = $status";}


# get R.A. and DEC for image center in first extension

Astro::FITS::CFITSIO::fits_movabs_hdu ($fptr, 2, $hdutype, $status);
if ($status) {die "Error moving to first extension:  status = $status";}
Astro::FITS::CFITSIO::fits_read_keyword ($fptr, "CRVAL1", $crval1, $comment,
   $status);
if ($status) {die "Error reading CRVAL1 in extension 1:  status = $status";}
Astro::FITS::CFITSIO::fits_read_keyword ($fptr, "CRVAL2", $crval2, $comment,
   $status);
if ($status) {die "Error reading CRVAL2 in extension 1:  status = $status";}
Astro::FITS::CFITSIO::fits_read_keyword ($fptr, "CRPIX1", $crpix1, $comment,
   $status);
if ($status) {die "Error reading CRPIX1 in extension 1:  status = $status";}
Astro::FITS::CFITSIO::fits_read_keyword ($fptr, "CRPIX2", $crpix2, $comment,
   $status);
if ($status) {die "Error reading CRPIX2 in extension 1:  status = $status";}
Astro::FITS::CFITSIO::fits_get_img_parm ($fptr, undef, $naxis, $naxes,
   $status);
if ($status) {die "Error reading size of extension 1:  status = $status";}
Astro::FITS::CFITSIO::fits_read_keyword ($fptr, "CD1_1", $cd1_1, $comment,
   $status);
$status = 0;
Astro::FITS::CFITSIO::fits_read_keyword ($fptr, "CDELT1", $cdelt1, $comment,
   $status);
$status = 0;
Astro::FITS::CFITSIO::fits_read_keyword ($fptr, "CD2_1", $cd2_1, $comment,
   $status);
$status = 0;

$pscale = sqrt ($cd1_1 ** 2 + $cd2_1 ** 2) * 3600;
if ($cdelt1) {$pscale = abs ($cdelt1 * 3600);}
if (!$pscale) {die "No CDELTn or CD matrix found";}

($ra0, $dec0) = &coordoffset ($crval1, $crval2,
   -(${$naxes}[0] / 2 - $crpix1) * $pscale / 3600,
   (${$naxes}[1] / 2 - $crpix2) * $pscale / 3600);


# determine number of extensions

Astro::FITS::CFITSIO::fits_get_num_hdus ($fptr, $next, $status);
if ($status) {die "Error getting number of extensions:  status = $status";}


# loop through extensions and determine offsets

for ($ext = 1; $ext < $next; $ext++) {

   Astro::FITS::CFITSIO::fits_movabs_hdu ($fptr, $ext+1, $hdutype, $status);
   if ($status) {die "Error moving to ext. $ext:  status = $status";}
   Astro::FITS::CFITSIO::fits_read_keyword ($fptr, "CRVAL1", $crval1, $comment,
      $status);
   if ($status) {die "Error reading CRVAL1 in ext. $ext:  status = $status";}
   Astro::FITS::CFITSIO::fits_read_keyword ($fptr, "CRVAL2", $crval2, $comment,
      $status);
   if ($status) {die "Error reading CRVAL2 in ext. $ext:  status = $status";}
   Astro::FITS::CFITSIO::fits_read_keyword ($fptr, "CRPIX1", $crpix1, $comment,
      $status);
   if ($status) {die "Error reading CRPIX1 in ext. $ext:  status = $status";}
   Astro::FITS::CFITSIO::fits_read_keyword ($fptr, "CRPIX2", $crpix2, $comment,
      $status);
   if ($status) {die "Error reading CRPIX2 in ext. $ext:  status = $status";}
   Astro::FITS::CFITSIO::fits_get_img_parm ($fptr, undef, $naxis, $naxes,
      $status);
   if ($status) {die "Error reading size of ext. $ext:  status = $status";}

   ($ra, $dec) = &coordoffset ($crval1, $crval2,
      -(${$naxes}[0] / 2 - $crpix1) * $pscale / 3600,
      (${$naxes}[1] / 2 - $crpix2) * $pscale / 3600);

   printf "%0.0f %0.0f\n",
      ($ra0 - $ra) * 3600 / $pscale * cos($dec0*$pi/180) * $mag,
      ($dec - $dec0) * 3600 / $pscale * $mag;

}




# SUBROUTINES

# add offset to R.A. and DEC. (in decimal degrees)

sub coordoffset {

   my ($old_ra, $old_dec, $ra_offset, $dec_offset) = @_;
   my ($new_ra, $new_dec);
   my $pi = 4 * atan2 (1, 1);

   $new_dec = $old_dec + $dec_offset;

   if ($old_dec == 90) {
      print "DEC = 90, so I won't update R.A.\n";
      $new_ra = $old_ra;
   }
   else {$new_ra = $old_ra + $ra_offset / cos (abs ($old_dec) * $pi / 180);}

   return ($new_ra, $new_dec);
}
