#!/usr/bin/perl
# Use starlink perl when running at the telescope on aatlxa
#
#
# ast3_crosstalk
# ==============
# Perl script to slice and dice AST3-1 images from 2012 to look
# at cross talk
#
#
# Produce a cross-talk image, that can be multiplied by a fudge factor
# and will result in an image you can add to your data file to
# correct cross talk
#

{
   #use diagnostics;
   use PDL;
   #use PDL::Graphics::PGPLOT;
   #use PGPLOT;
   use PDL::IO::Misc;
   use PDL::IO::NDF;
   use PDL::NiceSlice;
   #use Math::Fortran qw(sign);
   use Astro::SLA;
 
   my $THIS = 'ast3_crosstalk';
   
   # Som trig stuff
   my $pi = atan2(1,1) * 4.0;
   my $piover180 = $pi/180.0;

   # Parse the provided arguments
   #
   my @command = @ARGV;

   my $VERBOSE = 0;
   #my $CHIP = '\d';
   
   my $usage="Usage: $THIS [-h] [-v] [-f fudge_factor] file [output_file] \n\n".
          "       -h display this help\n".
          "       -v be verbose and display debugging information\n".
          "       -f fudge factor that will be multipled by cross-talk image\n\n".
          "       if output file is not specified only the file 'crosstalk.fits' will be written\n".
          "       if output file is specified a corrected version of the input file will be created\n".
          "       \n";



   $fudge = 1;
   while ( $ARGV[0] =~ /^-/ ) {
      if ( $ARGV[0] =~ /^-h$/){
         shift @ARGV;
         print "$usage\n";
         exit;
         }
      elsif ( $ARGV[0] =~ /^-v$/){
         shift @ARGV;
         $VERBOSE = 1;
         }
      elsif ( $ARGV[0] =~ /^-f$/){
         shift @ARGV;
         $fudge = shift @ARGV;
         }
      else {
         print "$THIS: Unknown option $ARGV[0] ignored\n";
         shift @ARGV;
         }
      }
   
   my $in  = shift @ARGV;
   my $out = '';
   if ( @ARGV ) {
      $out = shift @ARGV;
      }

   # Make sure both direcories exist. If the output one doesn't, then create it
   if ( ! -e $in ) {
      die "$THIS: File $in cannot be found\nAborting\n\n$usage\n";
      }
   
 
   my $prefix;
   my $num;
   my $loop;
   my $chip;
   my $nloops;
   # Force output FITS files to be in REAL
   my $BITPIX = -32;
   my $notgotall = 0;
   my $try;
   my $pad;
   my @try;
   my @PDLS;
   my $l;
   my $t;
   my $p;
   my $pdl2;
   my $pdl_out;
   
   my @CHIPS = ();
   #my $tr = pdl(0.0,1.0,0.0,0.0,0.0,1.0);
   #dev('/xs');
   #pgpap(5.0,1.5);
   #env(0.0,1500.0,0.0,2650.0);

   # Mask the input image at 65535
   $pdl_in = rfits($in,{hdrcpy=>1})->float >= 65534.0;
   $pdl_in *= 65535;
   
   print "Masked $in at 65535\n" if $VERBOSE;


   # Slice the input file into 16 sub-images. Note this takes
   # advantage of perl's slicing system so that each slice is
   # simply an indexed image that points back to the original mask.

   $c = 1;
   $xsize = 1500;
   $ysize = 2650;
   $xs = 0;
   $ys = 0;
   foreach $i ( 0 .. 7 ) {
      $xe = $xs + $xsize - 1;
      $ye = $ys + $ysize - 1;
      $slice = slice $pdl_in, "$xs:$xe,$ys:$ye";
      print "Sliced $in $xs:$xe,$ys:$ye => CCD $c\n" if $VERBOSE;
      #($mean,$rms,$median,$min,$max) = $slice->stats;
      #print " - slice has\$mean,\$rms,\$median,\$min,\$max = $mean,$rms,$median,$min,$max\n";
      if ( $c == 1 ) {
         $sum = $slice->copy;
         }
      else {
         $sum += $slice;
      }
      $xs = $xs + $xsize;
      $c++;
      #$low = 0.0;
      #$high = 200000.0;
      #imag $sum, $low, $high,$tr;
      push @CHIPS, $slice;
      #pgcurs($xin, $yin, $character);
      }

   # Top row - note reverse if Y direction for the slice as it needs to be
   # flipped.
   $xs = 0;
   $ys = $ysize;
   foreach $i ( 0 .. 7 ) {
      $xe = $xs + $xsize - 1;
      $ye = $ys + $ysize - 1;
      $slice = slice $pdl_in, "$xs:$xe,$ye:$ys";
      print "Sliced $in $xs:$xe,$ys:$ye => CCD $c\n" if $VERBOSE;
      #($mean,$rms,$median,$min,$max) = $slice->stats;
      #print " - slice has\$mean,\$rms,\$median,\$min,\$max = $mean,$rms,$median,$min,$max\n";
      $sum += $slice;
      $xs = $xs + $xsize;
      $c++;
      #$low = 0.0;
      #$high = 200000.0;
      #imag $sum, $low, $high,$tr;
      push @CHIPS, $slice;
      #pgcurs($xin, $yin, $character);
      }

   wfits $sum, "SUM.fits" if $VERBOSE;

   #
   # Now we want to do a median over those 16 images to try and suppress the
   # star light and just leave the cross talk.
   # $stack = pdl(@CHIPS)->reorder(2,0,1);
   # $median = medover $stack;
   # wfits $median, "MEDIAN.fits";
   #

   # Now use the same procedure to repopulate copies of this 'summed' detector
   # back out into all the other detectors.

   foreach $c ( @CHIPS ) {
      $c .= $sum;
      }

   #
   # Apply the fudge factor if requested
   #
   if ( $fudge != 1 ) {
      $pdl_in *= $fudge;
      }

   # write the cross talk file
   wfits $pdl_in, 'crosstalk.fits';
   print "Crosstalk file with fudge factor $fudge written to crosstalk.fits\n" if $VERBOSE;

   # If an output file has been requested, apply the correction to the
   # input file and write that.
   
   if ( $out ne '' ) {
      $pdl_out = rfits($in,{hdrcpy=>1});
      $pdl_out += $pdl_in;
      wfits $pdl_out, $out;
      print "Crosstalk corrected file written to $out\n" if $VERBOSE;

      }



}
   
