#! /bin/sh
# This is the LHEA perl script: /cvmfs/extras-fp7.egi.eu/extras/heasoft/heatools/x86_64-unknown-linux-gnu-libc2.19-0/bin/ftcoco
# The purpose of this special block is to make this script work with
# the user's local perl, regardless of where that perl is installed.
# The variable LHEAPERL is set by the initialization script to
# point to the local perl installation.
#-------------------------------------------------------------------------------
eval '
if [ "x$LHEAPERL" = x ]; then
  echo "Please run standard LHEA initialization before attempting to run /cvmfs/extras-fp7.egi.eu/extras/heasoft/heatools/x86_64-unknown-linux-gnu-libc2.19-0/bin/ftcoco."
  exit 3
elif [ "$LHEAPERL" = noperl ]; then
  echo "During LHEA initialization, no acceptable version of Perl was found."
  echo "Cannot execute script /cvmfs/extras-fp7.egi.eu/extras/heasoft/heatools/x86_64-unknown-linux-gnu-libc2.19-0/bin/ftcoco."
  exit 3
elif [ `$LHEAPERL -v < /dev/null 2> /dev/null | grep -ic "perl"` -eq 0 ]; then
  echo "LHEAPERL variable does not point to a usable perl."
  exit 3
else
  # Force Perl into 32-bit mode (to match the binaries) if necessary:
  if [ "x$HD_BUILD_ARCH_32_BIT" = xyes ]; then
    if [ `$LHEAPERL -V 2> /dev/null | grep -ic "USE_64_BIT"` -ne 0 ]; then
      VERSIONER_PERL_PREFER_32_BIT=yes
      export VERSIONER_PERL_PREFER_32_BIT
    fi
  fi
  exec $LHEAPERL -x $0 ${1+"$@"}
fi
'
if(0);
# Do not delete anything above this comment from an installed LHEA script!
#-------------------------------------------------------------------------------
#!/usr/bin/perl
#
# ftcoco
#
# Type "fhelp ftcoco" for more help.
#
# $Log: ftcoco,v $
# Revision 1.7  2011/03/23 20:22:26  craigm
# Version 2.3; request of Bill Pence and Web Hera team to have infile='' mean the same as infile='NONE'; set default chatter to chatter=2 so that console input prints out; unit tests still pass --CM
#
# Revision 1.6  2008/09/18 01:11:39  craigm
# Update to version 2.2; now output units are correctly specified (before the TUNIT# notation was triggering CFITSIO 'hierarchical' keywords; this has been fixed); unit test is unchanged --CM
#
# Revision 1.5  2006/06/01 07:15:57  craigm
# The task no longer removes the existing columns (but this also means that if RA/DEC columns already exist, they must be floating point; add TUNITs; unit tests still pass --CM
#
# Revision 1.4  2006/05/26 06:33:29  craigm
# Remove a debugging statement--CM
#
# Revision 1.3  2006/05/26 03:07:38  craigm
# Bump to version 2.1 to reflect previous change --CM
#
# Revision 1.2  2006/05/26 03:03:31  craigm
# Make sure the FK4/FK5 selector actually works;  I have no idea how this slipped through the cracks;  the unit test now works again; bump to version 2.1 --CM
#
# Revision 1.1  2006/03/26 20:35:39  craigm
# Initial commit (version 2.0), which is a translation from galeq --CM
#
#
#
#
$taskname = "ftcoco";
$taskver  = "2.3";
$id = "$taskname v$taskver";

use HEACORE::HEAINIT;
use POSIX qw(floor);
use Math::Trig qw(asin);

# ===================================
# Execute main subroutine, with error trapping
$status = 0;

eval {
  $status = headas_main(\&ftcoco);
};

# ===================================
# Check for errors and report them to the console
if ($@) {
    if ($status == 0) { $status = -1; }
    warn $@;
    exit $status;
}
exit 0;

# ===================================
# Main subroutine
sub ftcoco {
    # Makes all environment variables available
    use Env;
    use Cwd;

    # The HEAUTILS module provides access to HDpar_stamp()
    # set_toolname(), set_toolversion(), and headas_clobberfile()
    use HEACORE::HEAUTILS;
    use HEACORE::PIL;
    # include the file specification functions
    use Astro::FITS::CFITSIO qw( :shortnames :constants );

    # Use the standard HEAdas methods for registering the toolname and version number to be
    # used in error reporting and in the record of parameter values written by HDpar_stamp
    set_toolname($taskname);
    set_toolversion($version);

    eval {
	$status = &ftcoco_work();
    };

    if ($@) {
	if ($status == 0) { $status = -1; }
	warn $@;
	return $status;
    }

    return $status;
}

#
# Subroutine ftcoco - build expressions to compute ra/dec or glon/glat
#
# This routine builds CFITSIO calculator expressions that can be used
# to convert between RA/DEC and GLON/GLAT.  The direction of the 
# conversion is determined by the $select parameter.  The input
# and output angles are defined as follows:
#                   $ai   $bi        $ao   $bo
#                  =============================
#    $select == 1: ($ra,  $dec)  -> ($glon,$glat)
#    $select == 2: ($glon,$glat) -> ($ra,  $dec)
# 
# Inputs
#    $ai = name of input longitude column (e.g. "RA_OBJ")
#    $bi = name of output latitude column (e.g. "DEC_OBJ")
#
# Returns: an array of two expressions.  $output[0] computes $ao and
#          $output[1] computes $bo.
#
sub mymin {
    my ($a,$b) = (@_);
    if ($a < $b) { return $a; }
    return $b;
}
sub mymod {
    my ($a,$b) = (@_);
    my ($c,$rem);
    $c = floor($a / $b);
    $rem = ($a - $b*$c);
    return $rem;
}
sub cfitsio_mod {
    my ($a,$b) = (@_);
    return "$a % $b";
}

sub ftcoco_calc {
    my ($ai,$bi,$select,$internal) = @_;
    
    $select = 1 unless(defined($select));
    $psi = $psi0[$select - 1];
    $phi = $phi0[$select - 1];
    $stheta = $stheta0[$select - 1];
    $ctheta = $ctheta0[$select - 1];
#    print "select=$select psi=$psi phi=$phi stheta=$stheta ctheta=$ctheta\n";

    $pi     = 3.1415926535897931;
    $twopi  = $pi * 2.0;
    $fourpi = $pi * 4.0;
    $deg_to_rad = 180.0 / $pi;

    if ($internal) {
	$min = "&mymin";
	$mod = "&mymod";
	$atan2 = "atan2";
	$asin = "asin";
    } else {
	$min = "min";
	$mod = "&cfitsio_mod";
	$atan2 = "arctan2";
	$asin = "arcsin";
    }

    $a    = "($ai/$deg_to_rad - ($phi))";
    $b    = "($bi/$deg_to_rad)";
    $sb   = "sin($b)"; $cb = "cos($b)";
    $cbsa = "($cb * sin($a))";
    $bb   = "(-($stheta) * $cbsa + ($ctheta) * $sb)";
    $bb   = "$min($bb,1.0)";
    $bo   = "($asin($bb)*$deg_to_rad)";

    $aa =  "$atan2( ($ctheta) * $cbsa + ($stheta) * $sb, $cb * cos($a) )";
    $ao = "($aa+($psi)+$fourpi)";
    if ($internal) {
	$ao = "($mod($ao,$twopi))";
    } else {
	$ao = "($ao % $twopi)";
    }
    $ao = "($ao*$deg_to_rad)";
    
    return ($ao,$bo);
}

sub ftcoco_work {
    use Astro::FITS::CFITSIO qw(:longnames :constants);
    use HEACORE::HEAUTILS;
    use HEACORE::PIL;
    $status = 0;

    if (PILGetString('infile', $infile) ||
	PILGetString('outfile', $outfile) ||
	PILGetString('incoord', $incoord) ||
	PILGetString('outcoord', $outcoord) ||
	PILGetString('lon1', $lon1) ||
	PILGetString('lat1', $lat1) ||
	PILGetString('radecsys', $radecsys) ||
	PILGetBool('clobber',$clobber) || 
	PILGetInt('chatter',$chatter) ||
	PILGetBool('history',$history)) {
	die "ERROR: could not read task parameters";
    }

    if ($infile =~ m/^none$/i || $infile eq "") { 
	$internal = 1; 
	print "Internal calculations\n" if ($chatter >= 5);
	    
    } else {
	$internal = 0;
	print "CFITSIO calculations\n" if ($chatter >= 5);
	if (PILGetString('lon2', $lon2) ||
	    PILGetString('lat2', $lat2)) {
	    die "ERROR: could not read lon2/lat2 parameters";
	}
    }	

    if ($clobber) { $clobber = "YES"; } else { $clobber = "NO"; }
    if ($history) { $history = "YES"; } else { $history = "NO"; }

    if ($outfile =~ m/^infile/i) {
	$outfile = "$infile";
    }

    $inequ = 0;  $ingal = 0;  $inecl = 0;
    $outequ = 0; $outgal = 0; $outecl = 0;
    $inequ = 1 if ($incoord =~ m/^r/i);
    $ingal = 1 if ($incoord =~ m/^g/i);
    $inecl = 1 if ($incoord =~ m/^e/i);

    $outequ = 1 if ($outcoord =~ m/^r/i);
    $outgal = 1 if ($outcoord =~ m/^g/i);
    $outecl = 1 if ($outcoord =~ m/^e/i);

    if    ($inequ && $outequ) { $select = 0; }  # Identity
    elsif ($ingal && $outgal) { $select = 0; }  # Identity
    elsif ($inecl && $outecl) { $select = 0; }  # Identity 
    elsif ($inequ && $outgal) { $select = 1; }  # EQU -> GAL
    elsif ($ingal && $outequ) { $select = 2; }  # GAL -> EQU
    elsif ($inequ && $outecl) { $select = 3; }  # EQU -> ECL
    elsif ($inecl && $outequ) { $select = 4; }  # ECL -> EQU
    elsif ($inecl && $outgal) { $select = 5; }  # ECL -> GAL
    elsif ($ingal && $outecl) { $select = 6; }  # GAL -> ECL
    else {
	die "ERROR: incoord and outcoord must be either 'R' or 'G'";
    }
    print "Output selector=$select radecsys=/$radecsys/ \n" if ($chatter >= 5);

    # FK4 SYSTEM
    @psi0_fk4    = ( 0.57595865315, 4.9261918136, 0.00000000000, 0.0000000000,
		     0.11129056012, 4.7005372834 );
    @stheta0_fk4 = ( 0.88781538514,-0.88781538514, 0.39788119938,-0.39788119938,
		     0.86766174755,-0.86766174755);
    @ctheta0_fk4 = ( 0.46019978478, 0.46019978478, 0.91743694670, 0.91743694670,
		     0.49715499774, 0.49715499774);
    @phi0_fk4    = ( 4.9261918136,  0.57595865315, 0.0000000000, 0.00000000000,
		     4.7005372834, 0.11129056012);

    # FK5 SYSTEM
    @psi0_fk5    = ( 0.57477043300, 4.9368292465, 
		     0.00000000000, 0.0000000000, 
		     0.11142137093, 4.71279419371);
    @stheta0_fk5 = ( 0.88998808748,-0.88998808748,
		     0.39777715593,-0.39777715593,
		     0.86766622025,-0.86766622025);
    @ctheta0_fk5 = ( 0.45598377618, 0.45598377618,
		     0.91748206207, 0.91748206207,
		     0.49714719172, 0.49714719172);
    @phi0_fk5    = ( 4.9368292465,  0.57477043300,
		     0.0000000000, 0.00000000000,
		     4.71279419371, 0.11142137093);

    if ($radecsys =~ m/^fk4$/i) {
	@psi0 = @psi0_fk4;
	@phi0 = @phi0_fk4;
	@stheta0 = @stheta0_fk4;
	@ctheta0 = @ctheta0_fk4;
    } elsif ($radecsys =~ m/^fk5$/i) {
	@psi0 = @psi0_fk5;
	@phi0 = @phi0_fk5;
	@stheta0 = @stheta0_fk5;
	@ctheta0 = @ctheta0_fk5;
    } else {
	die "ERROR: radecsys must be either FK4 or FK5";
    }

    if ($select > 0) {
	@exprs = &ftcoco_calc("$lon1", "$lat1", $select, $internal);
    } else {
	@exprs = ("$lon1", "$lat1");
    }

    if ($chatter >= 5) {
	print "$out1 = $exprs[0]\n";
	print "$out2 = $exprs[1]\n";
    }

    if ($internal) {
	# Compute things internally within perl
	$cmd1 = "\$lon2 = $exprs[0];"; print "$cmd1\n" if ($chatter >= 5);
	$cmd2 = "\$lat2 = $exprs[1];"; print "$cmd2\n" if ($chatter >= 5);
	eval "$cmd1";
	if (! $@) { eval "$cmd2"; }
	if (! $@) {
	    $status = PILPutString('lon2',"$lon2");
	    # && warn "ERROR: could not save lon2";
	    $status = PILPutString('lat2',"$lat2");
	    # && warn "ERROR: could not save lat2";
	    print "$lon2 $lat2\n" if ($chatter >= 2);
	    return 0;
	} 
	warn $@ if $@;
	return -1;

    } else {
	$outlon2 = "$lon2"."(D)=$exprs[0]; #TUNIT#( physical unit of field) = \"deg\";";
	$outlat2 = "$lat2"."(D)=$exprs[1]; #TUNIT#( physical unit of field) = \"deg\";";

	# NOTE: the following expression is commented out.
	# Following expressions make sure that the output columns are
	# fully erased before trying to create new ones.  It avoids the
	# situation where we reuse a column called GLON_OBJ which is an
	# integer.  It still doesn't seem to work though.
	# NOTE: this is commented out.
# 	$secexpr = "$lon2=0;$lat2=0;-$lon2;-$lat2;";

	$secexpr = "";

	$cmd = "ftcopy '$infile"."[col *;$secexpr $outlon2 $outlat2]' '$outfile' clobber=$clobber history=$history chatter=$chatter";
	print "$cmd\n" if ($chatter >= 5);

	return system($cmd) / 256;
	
    }

}

# ====================================================================
# ====================================================================
# ====================================================================
# ================ End of MAIN =======================================