#! /bin/sh
# This is the LHEA perl script: /cvmfs/extras-fp7.egi.eu/extras/heasoft/ftools/x86_64-unknown-linux-gnu-libc2.19-0/bin/proc2pcu
# 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/ftools/x86_64-unknown-linux-gnu-libc2.19-0/bin/proc2pcu."
  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/ftools/x86_64-unknown-linux-gnu-libc2.19-0/bin/proc2pcu."
  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
#
# File: proc2pcu
#
# Description:  Preprocess RXTE 2PCU mode
#
# Author: M. Tripicco (ADNET Systems)
# Date:   2009-08-27
#
# v1.1: added checksum/datasum update
#

use Getopt::Std;

getopts(':h');
if (defined $opt_h) {
    print <<EOHELP;

NAME
    proc2pcu - Pre-process RXTE/PCA binned-mode data files containing
               two separate channels of PCU data to produce files that
               are compatible with saextrct.

USAGE
    proc2pcu infile outfile pcugroup (clobber) (chatter) (history)

DESCRIPTION
    Several new RXTE PCA data modes were introduced in mid-2007. The new binned
    mode "B_250us_2PCU_0_249_H" divides the PCUs into two groups (1,3,5 and 2,4)
    and files containing data in this mode cannot be processed by saextrct. This
    script preprocesses these files and generates output that can be successfully
    handled by saextrct. One can output either the PCU 1,3,5 group, the PCU 2,4
    group, or the sum of both groups. The original XeCnt data column is preserved
    in the output file as OrigCntXe.

PARAMETERS

    infile [filename]
        Name of the PCA data file (must be 2PCU mode) to be processed

    outfile [filename]
        Name of the output file (suitable for processing by saextrct)

    pcugroup [string]
        Specifies which data should be written to the output file.
        Allowed values are:
           "1" or "024" or "135" (group 1; data from PCUs 0, 2, and 4)
           "2" or "13" or "24" (group 2; data from PCUs 1 and 3)
           "3" or "all" or "sum" (group 3; data from all PCUs)

    (chatter = 1) [int] range 0-5
        Verbosity level of output

    (clobber = no) [boolean]
        Overwrite output file?

    (history = yes) [boolean]
        Write standard HEADAS parameter history into output file?

EXAMPLES
    1. proc2pcu
    2. proc2pcu infile=FS37_1c330360-1c33099c.gz outfile=pcu024 pcugroup=1 clobber=yes

NOTES
    To date the relevant PCA data mode has been only used rarely, primarily for
    observations of Sco X-1 in proposal P93067.

BUGS
    Please report problems to xtehelp\@athena.gsfc.nasa.gov.

SEE ALSO
    SAEXTRCT

EOHELP
exit;
    }

use HEACORE::HEAINIT;

# ==================================================================
# Call the main task subroutine with an exception handler
$status = 0;
eval {
    $status = headas_main(\&proc2pcu);
};

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


# ===================================
# Main subroutine
sub proc2pcu {

    $taskname = "proc2pcu";
    $taskvers = "1.1";

    # 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($taskvers);

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

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


    return $status;
}

# ==================================================================
# Working subroutine
sub proc2pcu_work {

    # 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 :longnames);
    # User defined module which contains the Perl-CFITSIO interface functions 
    use SimpleFITS;

    my ($infile, $outfile, $chatter);
    my ($fits1, $fits2, $handle1, $handle2);
    my ($mode, $tddes5);

    ($status = PILGetFname('infile', $infile)) == 0 || die "error getting infile parameter";
    ($status = PILGetFname('outfile', $outfile)) == 0 || die "error getting outfile parameter";
    ($status = PILGetString('pcugroup', $pcugroup)) == 0 || die "error getting pcugroup parameer";
    ($status = PILGetString('clobber', $clobber)) == 0 || die "error getting clobber parameter";
    ($status = PILGetInt('chatter', $chatter)) == 0 || die "error getting chatter parameter";

    print "Running $taskname v$taskvers\n" if ($chatter >= 1);
    print "----------------------------------------\n" if ($chatter >= 2);

    if ($pcugroup =~ /(024|135)/ || $pcugroup eq "1"){
      $pcugroup = "XeCnt = OrigCntXe[1];";
      $newmode = "#DATAMODE = 'B_250us_1M_0_249_PCU135';";
      $descriptor = '#TDDES#  = "D[0^2^4] & E[X1L^X1R^X2L^X2R^X3L^X3R] & C[0~249]";';
      $desckey = 'D\[0\^2\^4\]';
      $longdescr = "D[0^2^4] & E[X1L^X1R^X2L^X2R^X3L^X3R] & C[0~249] & T[0.0;0.000244140625;8192]";
    }elsif ($pcugroup =~ /(24|13)/ || $pcugroup eq "2"){
      $pcugroup = "XeCnt = OrigCntXe[2];";
      $newmode = "#DATAMODE = 'B_250us_1M_0_249_PCU24';";
      $descriptor = '#TDDES#  = "D[1^3] & E[X1L^X1R^X2L^X2R^X3L^X3R] & C[0~249]";';
      $desckey = 'D\[1\^3\]';
      $longdescr = "D[1^3] & E[X1L^X1R^X2L^X2R^X3L^X3R] & C[0~249] & T[0.0;0.000244140625;8192]";
    }elsif ($pcugroup =~ /(all|sum)/i || $pcugroup eq "3"){
      $pcugroup = "XeCnt = OrigCntXe[1] + OrigCntXe[2];";
      $newmode = "#DATAMODE = 'B_250us_1M_0_249_H';";
      $descriptor = '#TDDES#  = "D[0~4] & E[X1L^X1R^X2L^X2R^X3L^X3R] & C[0~249]";';
      $desckey = 'D\[0\~4\]';
      $longdescr = "D[0~4] & E[X1L^X1R^X2L^X2R^X3L^X3R] & C[0~249] & T[0.0;0.000244140625;8192]";
    }else {die "ERROR: cannot parse pcugroup parameter"}

    @colfilt1 = qw(
		   *;
		   OrigCntXe == XeCnt;
		   -TDDES#;
		  );
    @colfilt2 = qw(
		   #TUNIT#  = 'count   ';
		   #TDISP#  = 'I5      ';
		  );
    @colfilt3 = qw(
		   #TDIM#   = '(8192,1)';
		   #1CTYP#  = 'TIME    ';
		   #1CUNI#  = 's       ';
		   #1CRPX#  = 1;
		   #1CRVL#  = 0.0;
		   #1CRVR#  = 'Time    ';
		   #1CDLT#  = 0.000244140625;
		   #2CTYP#  = 'CHANNEL ';
		   #2CPIX#  = '0~249';
		  );
    @colfilt=(@colfilt1,$pcugroup,$newmode,@colfilt2,$descriptor,@colfilt3);

    $fits1 = SimpleFITS->open("<$infile",ext=>2);
    die "ERROR: could not open $infile" if (! $fits1);
    $mode = $fits1 -> readkey("DATAMODE");
    unless ($mode){
      print "ERROR: could not read DATAMODE\n";
      $status = -1;
      goto CLEANUP;
    }
    if ($mode !~ /_2PCU_/){
      print "ERROR: Not a 2PCU DATAMODE ($mode)\n";
      $status = -1;
      goto CLEANUP;
    }
    print "Input file DATAMODE is $mode\n" if ($chatter > 4);
    $fits1->setstatus(0)->close();

    $infiltered = "$infile\[col ".join(" ",@colfilt)."\]";

    $fits1 = SimpleFITS->open("<$infiltered");
    die "ERROR: could not open $infile" if (! $fits1);
    $handle1 = $fits1->handle();
    print "Opened $infiltered\n" if ($chatter >= 3);

    &headas_clobberfile($outfile);
    if (-e $outfile){
      print "ERROR: output file $outfile already exists\n";
      $status = -1;
      goto CLEANUP;
    }
    $fits2 = SimpleFITS->open("$outfile",access=>"create");
    if (! $fits2){
      print "ERROR: could not open $outfile";
      goto CLEANUP;
    }
    $handle2 = $fits2->handle();
    print "Created $outfile\n" if ($chatter >= 3);

    &fits_copy_file($handle1,$handle2,1,1,1,$status);
    if ($status){
      print "ERROR: problem copying to $outfile\n";
      goto CLEANUP;
    }

    # Update new TDDES keyword (only way to deal with long string)
    $fits2 -> move(2);
    # Should always be column5 but check the descriptor first (else leave it alone)
    $tddes5 = $fits2 -> readkey("TDDES5");
    if ($tddes5 =~ /$desckey/){&ffmkls($handle2,'TDDES5',$longdescr,'&',$status)}

    # Update checksum/datasum for CHDU
    &fits_write_chksum($handle2,$status);

    $fits2 -> move(1);
    if (! $fits2){
      print "ERROR: could not move to extension 1\n";
      goto CLEANUP;
    }

    &HDpar_stamp($handle2,0,$status);

    # Update checksum/datasum for primary
    &fits_write_chksum($handle2,$status);

  CLEANUP:
    if ($fits1){$fits1->setstatus(0)->close()}
    if ($fits2){$fits2->setstatus(0)->close()}

    return $status;
}