#!/usr/bin/perl
use strict;
use warnings;

# 
# Makes images for event files supplied by user over entire FOV. 
# If a region file is also given, an image will also be made for the region. 
#
#
# Requirements:
#    0) must be in same directory as input event files. 
#    1) text file which lists the events files to process, one per line
#    2) If the region file exists, it should have the following format: 
#
#       unit
#       x y inner_radius outer_radius 
#
#    where unit = unit the coordinates and radii are in; can be either d (decimal degrees) 
#                  or p (pixels) 
#          inner_radius = inner radius of annulus region to extract. Set this 
#                         to 0 to extract a circular region.
#          outer_radius = outer radius of annulus region to extract. This is 
#                         equal to the radius of a circular region. 

my $currDir;
chomp ($currDir = `pwd`);
my $unittype = 'N';
my $currLine;
my @numbers;
my @files;
my $file;
my $m1_evt;
my $m2_evt;
my $pn_evt;
my $position_file=$currDir."/"."filler";
my $status=0;
my $process_pn = 0;
my $process_m1 = 0;
my $process_m2 = 0;
my $j;
my $inst;
my $coordflag;

# Set the environment variables.

my $baseDir;
if ($currDir =~ /^(.*?[0-9]{10})/) {
    $baseDir = $1 ;
}

if (-d $baseDir."/ODF") { 
   $ENV{SAS_ODFPATH} = $baseDir."/ODF";
   my $sumsas = glob ($baseDir."/ODF/*SUM.SAS");
   $ENV{SAS_ODF} = $sumsas;
   $ENV{SAS_CCF} = $baseDir."/ODF/ccf.cif";
} 
if (-d $baseDir."/odf") {
   $ENV{SAS_ODFPATH} = $baseDir."/odf";
   my $sumsas = glob ($baseDir."/odf/*SUM.SAS");
   $ENV{SAS_ODF} = $sumsas;
   $ENV{SAS_CCF} = $baseDir."/odf/ccf.cif";
}

# On to the main event.


print "Enter the file listing the event files: ";
chomp (my $event_files = <> );
$event_files = $currDir."/".$event_files;

print "Should I use a region file? If so, enter the file name; otherwise, enter N. ";
chomp (my $ans = <>) ;
if (($ans eq "N") or ($ans eq "n")) {
   $coordflag = 0 ;
} else {
   $coordflag = 1;
   $position_file = $currDir."/".$ans;
}

if (-e $event_files) {
   open(EVENTS, "<", $event_files);
   while ($currLine = <EVENTS>) {
      $currLine = $currDir."/".$currLine;
      chomp ($currLine);
      @files = split (/\s/, $currLine);
      foreach $file (@files) {  # Check the headers to find out which instruments we need.
         $inst = $currDir."/"."instru.list";
         `fkeyprint $file instrume outfile=$inst`;
         open (INSTRU, "<", $inst);
         while ($currLine = <INSTRU>) {
            if ($currLine =~ m/INSTRUME= 'EPN     '/ ) {
               $pn_evt = $file;
               $process_pn = 1;
            }
            if ($currLine =~ m/INSTRUME= 'EMOS1   '/ ) {
               $m1_evt = $file;
               $process_m1 = 1;
            }
            if ($currLine =~ m/INSTRUME= 'EMOS2   '/ ) {
               $m2_evt = $file;
               $process_m2 = 1;
            }
         }
         close (INSTRU);
         `rm $inst`;
      }
   }
   close (EVENTS);
}
else {
   print "No events file found; stopping here. \n" ;
   exit;
}

my @xval;
my @yval;
my @innerRad;
my @outerRad;
my @vals;
my $countlines=0;

# Does the coord file exist and have proper format?

if ($coordflag) { 
   if (-e $position_file) {
      open(COORDS, "<", $position_file);
      while ($currLine = <COORDS>) {
         if ($currLine =~ m/D/i) {
            $unittype = 'D';              # decimal degrees
         }
         if ($currLine =~ m/P/i ) {
            $unittype = 'P';              # pixel
         }
         if (($unittype ne 'D') and ($unittype ne 'P')) {
            print "No units found in the coords file. Exiting. \n";
            exit;
         }
         if ($currLine =~ m/\d/) {
            $countlines ++;
            # count numbers to confirm the correct format assign to variables
            @numbers = split (/\s/, $currLine);
            if (@numbers != 4) {
               print "The input file does not have the correct format. \n";
               exit; 
            }
            push (@xval, $numbers[0]);
            push (@yval, $numbers[1]);
            push (@innerRad, $numbers[2]);
            push (@outerRad, $numbers[3]);
         }
      }
      close (COORDS);
   }
   else {
      print "No region file found. \n"; 
      exit; 
   } 
}

#
# If necessary, change the degree decimals to pixel values. We don't know
# which detectors the user wants, so we have to prepare for each case.
#

my $convolved = 0;

if ($unittype eq 'D' ) {
   if (($process_m1 == 1) and ($convolved == 0)) {
      `evselect table=$m1_evt withimageset=yes imageset=m1_img_hold.fits xcolumn=X ycolumn=Y imagebinning=binSize ximagebinsize=32 yimagebinsize=32`;
      for ($j=0; $j<$countlines; $j++) {
         `ecoordconv imageset=m1_img_hold.fits x=$xval[$j] y=$yval[$j] coordtype=eqpos > hold.txt`;
         open (HOLD, "<", "hold.txt");
         while ($currLine = <HOLD>)  {
            if ($currLine =~ m/ X: Y: / ) {
               @vals = split (/\s/, $currLine);
               $xval[$j] = $vals[3];
               $yval[$j] = $vals[4];
            }
         }
         close (HOLD);
         $innerRad[$j] = $innerRad[$j]/1.38889e-5;
         $outerRad[$j] = $outerRad[$j]/1.38889e-5;
      }
      `rm hold.txt`;
      `rm m1_img_hold.fits`;
      $convolved = 1;
   }
   if (($process_m2 == 1) and ($convolved == 0)) {
      `evselect table=$m2_evt withimageset=yes imageset=m2_img_hold.fits xcolumn=X ycolumn=Y imagebinning=binSize ximagebinsize=32 yimagebinsize=32`;
      for ($j=0; $j<$countlines; $j++) {
         `ecoordconv imageset=m2_img_hold.fits x=$xval[$j] y=$yval[$j] coordtype=eqpos > hold.txt`;
         open (HOLD, "<", "hold.txt");
         while ($currLine = <HOLD>)  {
            if ($currLine =~ m/ X: Y: / ) {
               @vals = split (/\s/, $currLine);
               $xval[$j] = $vals[3];
               $yval[$j] = $vals[4];
            }
         }
         close (HOLD);
         $innerRad[$j] = $innerRad[$j]/1.38889e-5;
         $outerRad[$j] = $outerRad[$j]/1.38889e-5;
      }
      `rm hold.txt`;
      `rm m2_img_hold.fits`;
      $convolved = 1;
   }
   if (($process_pn == 1) and ($convolved == 0)) {
      `evselect table=$pn_evt withimageset=yes imageset=pn_img_hold.fits xcolumn=X ycolumn=Y imagebinning=binSize ximagebinsize=32 yimagebinsize=32`;
      for ($j=0; $j<$countlines; $j++) {
         `ecoordconv imageset=pn_img_hold.fits x=$xval[$j] y=$yval[$j] coordtype=eqpos > hold.txt`;
         open (HOLD, "<", "hold.txt");
         while ($currLine = <HOLD>)  {
            if ($currLine =~ m/ X: Y: / ) {
               @vals = split (/\s/, $currLine);
               $xval[$j] = $vals[3];
               $yval[$j] = $vals[4];
            }
         }
         close (HOLD);
         $innerRad[$j] = $innerRad[$j]/1.38889e-5;
         $outerRad[$j] = $outerRad[$j]/1.38889e-5;
      }
      `rm hold.txt`;
      `rm pn_img_hold.fits`;
   }
}

if (! -d $currDir.'/image') { 
   `mkdir image`; 
} 

my $m1_evt_cp;
my $m2_evt_cp;
my $pn_evt_cp;


if ($process_m1 eq 1) { 
   $m1_evt_cp = "image/m1_evt_cp";
   `cp $m1_evt $m1_evt_cp`;
}
if ($process_m2 eq 1) {
   $m2_evt_cp = "image/m2_evt_cp";
   `cp $m2_evt $m2_evt_cp`;
}
if ($process_pn eq 1) {
   $pn_evt_cp = "image/pn_evt_cp";
   `cp $pn_evt $pn_evt_cp`;
}

chdir $currDir.'/image';

if ((! $coordflag) or ($coordflag)) {
   if ($process_m1) {
      print "Making the M1 FOV image... \n";
      `evselect table=m1_evt_cp withimageset=yes imageset=m1_img.fits xcolumn=X ycolumn=Y imagebinning=binSize ximagebinsize=32 yimagebinsize=32`;
   }
   if ($process_m2) {
      print "Making the M2 FOV image... \n";
      `evselect table=m2_evt_cp withimageset=yes imageset=m2_img.fits xcolumn=X ycolumn=Y imagebinning=binSize ximagebinsize=32 yimagebinsize=32`;
   } 
   if ($process_pn) { 
      print "Making the PN FOV image... \n";
      `evselect table=pn_evt_cp withimageset=yes imageset=pn_img.fits xcolumn=X ycolumn=Y imagebinning=binSize ximagebinsize=32 yimagebinsize=32`;
   }
}

if ($coordflag) { 
   if ($process_m1) {
      for ($j=0; $j<$countlines; $j++) { 
         print "Making the M1 image $j... \n";
         `evselect table=m1_evt_cp withimageset=yes imageset=m1_img_src$j.fits xcolumn=X ycolumn=Y imagebinning=binSize ximagebinsize=32 yimagebinsize=32 expression='(X,Y) IN circle($xval[$j],$yval[$j],$outerRad[$j])'`;
      }
   }
   if ($process_m2) {
      for ($j=0; $j<$countlines; $j++) { 
         print "Making the M2 image $j... \n";
         `evselect table=m2_evt_cp withimageset=yes imageset=m2_img_src$j.fits xcolumn=X ycolumn=Y imagebinning=binSize ximagebinsize=32 yimagebinsize=32 expression='(X,Y) IN circle($xval[$j],$yval[$j],$outerRad[$j])'`;
      }
   }
   if ($process_pn) {
      for ($j=0; $j<$countlines; $j++) {
         print "Making the PN image $j... \n";
         `evselect table=pn_evt_cp withimageset=yes imageset=pn_img_src$j.fits xcolumn=X ycolumn=Y imagebinning=binSize ximagebinsize=32 yimagebinsize=32 expression='(X,Y) IN circle($xval[$j],$yval[$j],$outerRad[$j])' `;
      }
   }
}

# Next to godliness.

if ($process_m1) {
   `rm m1_evt_cp`;
}
if ($process_m2) {
   `rm m2_evt_cp`;
}
if ($process_pn) {
   `rm pn_evt_cp`;
}

print "Done. \n";