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

# 
# Makes light curve over entire FOV for event files supplied by user. 
# If a region file is also given, the light curve will be extracted 
# over both the full FOV and 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 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 $count = 0;
my $currLine;
my @numbers;
my $num;
my @files;
my $file;
my $m1_evt; 
my $m2_evt; 
my $pn_evt; 
my $j;
my $time_bin;
my $status;
my $inst;

my $process_pn = 0;
my $process_m1 = 0;
my $process_m2 = 0;
my $lcbin = 0 ; 
my $position_file;

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_file = <> );
$event_file = $currDir."/".$event_file;

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;
}

print "Please enter the size of the time bins. ";
chomp ($lcbin = <>) ;


if (-e $event_file) {
   open(EVENTS, "<", $event_file);
   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. \n" ;
   exit;
}

#print "$pn_evt\n";
#print "$m1_evt\n";
#print "$m2_evt\n";

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

# Does the coord file exist and have proper format?

if ($coordflag eq 1) {
   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);
   }
}

#
# 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.'/lightcurve') { 
   `mkdir lightcurve`; 
} 

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



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

chdir $currDir.'/lightcurve';

if (($coordflag eq 0) or ($coordflag eq 1)) {  # make the light curve for whole FOV
   if ($process_m1 eq 1) {
      `evselect table=m1_evt_cp withrateset=yes rateset='m1_ltcrv_bin$lcbin.fits' maketimecolumn=yes timecolumn=TIME timebinsize=$lcbin makeratecolumn=yes`;
   }
   if ($process_m2 eq 1) {
      `evselect table=m2_evt_cp withrateset=yes rateset=m2_ltcrv_bin$lcbin.fits maketimecolumn=yes timecolumn=TIME timebinsize=$lcbin makeratecolumn=yes`;
   } 
   if ($process_pn eq 1) { 
      `evselect table=pn_evt_cp withrateset=yes rateset=pn_ltcrv_bin$lcbin.fits maketimecolumn=yes timecolumn=TIME timebinsize=$lcbin makeratecolumn=yes`;
   }
}


if ($coordflag eq 1) {   # make the light curves for specific sources
   if ($process_m1 eq 1) {
      for ($j=0; $j<$countlines; $j++)  {
         `evselect table=m1_evt_cp withrateset=yes rateset="m1_ltcrv_src$j\_bin$lcbin.fits" maketimecolumn=yes timecolumn=TIME timebinsize=$lcbin makeratecolumn=yes expression='(X,Y) IN circle($xval[$j],$yval[$j],$outerRad[$j])'`;
      }
   }
   if ($process_m2 eq 1) {
      for ($j=0; $j<$countlines; $j++)  {
         `evselect table=m2_evt_cp withrateset=yes rateset=m2_ltcrv_src$j\_bin$lcbin.fits maketimecolumn=yes timecolumn=TIME timebinsize=$lcbin makeratecolumn=yes expression='(X,Y) IN circle($xval[$j],$yval[$j],$outerRad[$j])'`;
      }
   }
   if ($process_pn eq 1) {
      for ($j=0; $j<$countlines; $j++)  {
         `evselect table=pn_evt_cp withrateset=yes rateset=pn_ltcrv_src$j\_bin$lcbin.fits maketimecolumn=yes timecolumn=TIME timebinsize=$lcbin makeratecolumn=yes expression='(X,Y) IN circle($xval[$j],$yval[$j],$outerRad[$j])' `;
      }
   }
}

# Let's be tidy now.

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

print "Done. \n";