#!/usr1/local/bin/perl5 -w

# Attempt to do binned version

# Binary phase resolved spectral fits to Fermi LAT data

###################################
# Robin Corbet (UMBC), 2014-08-29 #
###################################

use strict 'subs';
use strict 'refs';

use POSIX ();

$VERSION = "1.02";

$MODEL_MAKER = "make4FGLxml";

$MODEL_MAKER_FLAG = "";

# counter if saving files for each bin
$BIN_COUNTER = 0;

#4FGL DR3
$DEFAULT_CATALOG = "gll_psc_v28.fit";


$help = "
like_bphase.pl

Modified version of like_lc.pl to do binary phase selected fits.

*** This script is _not_ suitable for *pulse* phase resolved fits.

Typing 
like_bphase.pl -help 
will display this help information.


like_bphase.pl -batch
runs in batch mode. (i.e. no prompting for parameters)

like_bphase.pl -file like_bphase.par

Reads parameters from the file specified (\"like_bphase.par\" in the example above).

like_bphase.pl -makepar like_bphase.par

This will create a file of the specified name (\"like_bphase.par\" in the example above) 
and write the current default parameters to it.

Perl script that generates phase-resolved Fermi LAT light curves using likelihood analysis. 
This is not an official part of the Fermi analysis software.

** Note, from version 1.01 and later *binned* likelihood is used by default, while
previously it was unbinned. It is possible to specify that unbinned analysis is used instead.

Caution. The models and other parameters used here are for demonstration purposes.
You should carefully evaluate the choice of parameters for your own particular needs!

To use this script you must be set up to 
run both the Fermi science tools and HEASOFT.
This script checks the environment variables LHEASOFT,
PFILES, and FERMI_DIR as a check on this.

In addition, if you are using the auto-generation option to automatically
make a model file from the FGL catalog, the environment variable MY_FERMI_DIR
needs to be set. This is where the FGL catalog and the
$MODEL_MAKER.py script must be located if used.

The model file generated by $MODEL_MAKER.py uses the parameters
in the FGL catalog. In this example only the parameters of the
source of interest are left free, others are held fixed at their
catalog values. Spectra are power laws. This may very well
not be optimum!

Among the inputs prompted for are two files: 
slist.dat
	This should contain a list of source names and their coordinates
	in the form:
	4FGLJ0240.5+6113 40.1317, 61.2281
	i.e. the 4FGL catalog name must be given because if you will use
	$MODEL_MAKER.py to generate the model file.
	Note that there is _no_ space in the source name.
	i.e. \"4FGLJ0240.5+6113\" rather than \"4FGL J0240.5+6113\"

Source names should not be duplicated because the output files are
constructed from the source names.

Optionally, an entry in the slist.dat file my have the form:
A_Source 40.1317, 61.2281 super_model.xml
where the 4th entry is the name of the model file to use for that source.
The source name (e.g. \"A_Source\") must be contained in the model file specified.
	

plist.dat
	This contains a list of photon files, one file per line.
It is also necessary to provide a single spacecraft (FT2) file.

Several temporary files are created while running this script.
At least most are deleted at the end.

The output file for each source is an ASCII file called
XXX_like_lc.dat, where \"XXX\" is the source name as given in slist.dat
This file contains:
time (bin center, phase),  flux (cts/cm^2/s), flux error, time bin half width (phase), TS value, [spectral parameters/errors]

The spectral parameters depend on the model used and are as taken directly from the gtlike output file.
There will thus generally be a repetition of the TS value and the flux at the end of the line.
Other values such as the \"ROI distance\" and \"Npred\" will also be included.
*** Output format changed in Version 0.21 and later! ***
The first of the spectral parameters is the normalization (or similar parameter). This was not
output in earlier versions. Other parameters such as Npred were also not previously output.
These changes were to allow the parameters of any model to be output. The change also
resulted in the repeat of TS and flux values.


Other command line options:


-start_phase [phase1]
	Specifies starting phase
-stop_phase [phase2]
	Specifies end phase
-start [mjdstart]
	The light curve will start at the specified time in units of MJD.
-stop  [mjdstop]
	The light curve will end at the specified time in units of MJD.
-debug
	Increase chatter level to \"2\" from \"0\"
-version
	Print version number of script then exit.
-batch
	Runs with minimal output and no prompting for inputs.
-makepar [filename]
	Generates a default parameter file with the name [filename]. If no name
	is given the file generated is like_lc.par
-file [filename]
	Reads parameter values from a file with the name [filename].
-prefix [prefix]
	Defines a new default prefix for output file names.
-unbinned
	Use unbinned likelihood analysis by default

CHANGES:
1.02 Use 4FGL DR3 catalog as default.
1.01 Use binned likelihood!
1.00 Select a phase range which enables pseudo-multiprocessing, and bumped up version number!
0.27 Use 4FGL DR2 catalog as default, updated IRFs to P8R3_SOURCE_V3_v1
0.26 Use 4FGL catalog as default, new relevant background files, P8R3 data
0.25 Added model file to parameter file (used if not created from e.g. 3FGL)
0.24 Option to save products (exposure cube etc.) to allow refitting
0.23 Use negative bin size to specify number of bins
0.22 Fixes for batch/file modes. Use temporary pfiles directory (avoid collisions)
0.21 Can output all spectral parameters for any spectral model
0.20 Support for Pass 8 data (default), and batch mode/parameter file usage
0.14 Read/write parameter files
0.13 Kludge to allow use of 3FGL catalog
0.12 Change to phase selection algorithm
0.1  First version of like_phase.pl - derived from like_lc.pl version 1.45

";



##########################################
# maximum and minimum functions
# taken from somewhere on web

sub max ($$) { $_[$_[0] < $_[1]] }
sub min ($$) { $_[$_[0] > $_[1]] }

$TRUE = 1;
$FALSE = 0;

$BATCH = $FALSE;

$DEBUG = $FALSE;


# default phase range
$start_phase = 0.0;
$stop_phase = 1.0;

# default minimum start time
$start_time = 0;
# default maximum stop time
$stop_time = 9e99;

# minimal output from tools
#$chatter = 0;
# standard output from tools
#$chatter = 2;
$DEFAULT_CHATTER = 0;


# Conversion from MET (s) to MJD (days, Terrestrial Time)
# $MJDREF = 51910.0 + 7.428703703703703E-4;

$DEFAULT_BIN_SIZE = 0.1;
$DEFAULT_PERIOD = -1;
$DEFAULT_T0 = -1;
$DEFAULT_ZENITH_LIMIT = 105;
$DEFAULT_ROCK = 90;
$DEFAULT_BORE = 180;
$DEFAULT_EVENT_CLASS_MIN = 3;
$DEFAULT_MODEL_FILE = "outmodel.xml";
$DEFAULT_SOURCE_LIST = "slist.dat";
$DEFAULT_FT1_LIST = "plist.dat";
$DEFAULT_IRF_CODE = 9; # 4FGL + new Galactic diffuse


$DEFAULT_SUN_MINIMUM = 5.0;

$DEFAULT_SAVE_PRODUCTS = $FALSE;

$DEFAULT_AUTO_MODEL = $TRUE;

$DEFAULT_USE_BINNED = $TRUE;


$DEFAULT_TS_THRESHOLD = 2;

# Region size
$DEFAULT_ROI = 10;

# Default, "hidden" parameters
$DEFAULT_EMIN = 100;
$DEFAULT_EMAX = 500000;

# prefix for output files
$DEFAULT_PREFIX = "";

$DEFAULT_PARAMETER_FILE = "like_bphase.par";

$catalog = $DEFAULT_CATALOG;

$DEFAULT_FT2 = "lat_spacecraft_merged.fits";


# Number of time slots to select with gtmktime
# for each call.
# (This is to reduce the number of calls to gtmktime.)
$nslots = 15;

print "like_bphase.pl: V$VERSION (R.H.D. Corbet)\n";

# Check command line parameters

$numArgs = $#ARGV + 1;
#if($numArgs > 0){print "$numArgs command-line arguments:\n"};

foreach $argnum (0 .. $#ARGV) {

   print "$ARGV[$argnum]\n";
	if($ARGV[$argnum] =~ /-help|-h|--help/) {
		print $help;
		exit;
	}
	elsif($ARGV[$argnum] eq "-version") {
		print "like_bphase.pl: Version $VERSION\n";
		exit;
	}


       elsif($ARGV[$argnum] =~ /-unbinned/) {
                print "Default set to unbinned analysis\n";
		$DEFAULT_USE_BINNED = $FALSE;
        }


        elsif($ARGV[$argnum] eq "-batch") {
                print "Batch mode - no parameter prompting, limted screen output\n";
                $BATCH = $TRUE;
		$chatter = 0;
                }


# needs to be last argument specified if the default is to be used
	elsif($ARGV[$argnum] eq "-file") {
		if($argnum+1 < $numArgs){
			$parameter_filename = $ARGV[$argnum+1];
			print "parameter_filename = $parameter_filename\n";
			set_defaults();
			print "default period = $DEFAULT_PERIOD\n";
		}
		else{
			print "No input parameter file name given, using default of $DEFAULT_PARAMETER_FILE\n";
			$parameter_filename = $DEFAULT_PARAMETER_FILE;
			set_defaults();
		}

	}

	elsif($ARGV[$argnum] eq "-makepar") {
		if($argnum+1 < $numArgs){
			$output_parameter_filename = $ARGV[$argnum+1];
			print "output parameter_filename = $output_parameter_filename\n";
			write_defaults();
		}
		else{
			print "No output parameter file name given, using default of $DEFAULT_PARAMETER_FILE\n";
			$output_parameter_filename = $DEFAULT_PARAMETER_FILE;
			write_defaults();
		}

	}
	elsif($ARGV[$argnum] eq "-start") {
		if($argnum+1 < $numArgs){
			$start_time = mjd2met($ARGV[$argnum+1]);
			print "start time = MET $start_time\n";
		}
		else{
			print "No input start time value given\n";
		}
	}
	elsif($ARGV[$argnum] eq "-stop") {
		if($argnum+1 < $numArgs){
			$stop_time = mjd2met($ARGV[$argnum+1]);
			print "stop time = MET $stop_time\n";
		}
		else{
			print "No input stop time value given\n";
		}

	}

	elsif($ARGV[$argnum] eq "-start_phase") {
		if($argnum+1 < $numArgs){
			$start_phase = $ARGV[$argnum+1];
			print "start phase = $start_phase\n";
		}
		else{
			print "No input start phase value given\n";
		}
	}
	elsif($ARGV[$argnum] eq "-stop_phase") {
		if($argnum+1 < $numArgs){
			$stop_phase = $ARGV[$argnum+1];
			print "stop phase = $stop_phase\n";
		}
		else{
			print "No input stop phase value given\n";
		}

	}






	elsif($ARGV[$argnum] =~ /-prefix/) {
		if($argnum+1 < $numArgs){
			$DEFAULT_PREFIX = $ARGV[$argnum+1];
			print "Default prefix = $DEFAULT_PREFIX\n";
		}
		else{
			print "No input start time value given\n";
		}

	}
	elsif($ARGV[$argnum] =~ /-debug/) {
		$chatter = 2;
		$DEBUG = $TRUE;

	}
}


##########################################################
# Check if FTOOLS are available
#################################

#if (not exists($ENV{'LHEASOFT'}) ||  not exists($ENV{'PFILES'})  ) 
#if($ENV{'LHEASOFT'} !~/\S/ || $ENV{'PFILES'} !~/\S/  ) 
#if (!defined($status = $ENV{'LHEASOFT'})  ||  !defined($status = $ENV{'PFILES'})  )
if (!defined($status = $ENV{'HEADAS'})  ||  !defined($status = $ENV{'PFILES'})  )
{
    print "\n Error: You need to set up HEASOFT to use this script.\n";
    print " e.g. typically you need to type \"ftools\".\n";
    exit(0);
}


# Check if environment variable FERMI_DIR set
#if($ENV{'FERMI_DIR'} !~/\S/ ) 
if (!defined($status = $ENV{'FERMI_DIR'} ))
{
    print "\n Error: You need to set up the Fermi Science Tools to use this script.\n";
    print " e.g. typically you need to type \"fermi\".\n";
    exit(0);
}



$chatter = get_value("Give chatter level (0 = none, 2=nominal, 4 = maximum)", $DEFAULT_CHATTER, 0, 4);


$irf_code = get_value("Which IRFs? (1 = P6_V3, 2 = P6_V11, 3 = P7_V6, 4 = P7REP, 5 = P7REP/newbackground, 6 = P8, 7 = P8R3 new iso, 8 = 4FGL, 9 = v2.0)", $DEFAULT_IRF_CODE, 1, 9);

if($irf_code == 1){
	$irf_prefix = "P6_V3_";
	$gll_file = "galdiffuse/gll_iem_v02.fit";
	$isotropic_file = "galdiffuse/isotropic_iem_v02.txt";
	$galactic_name = "gal_v02";
	$extra_galactic_name = "eg_v02";
	}

elsif($irf_code == 2){
	$irf_prefix = "P6_V11_";
	$gll_file = "galdiffuse/gll_iem_v02_P6_V11_DIFFUSE.fit";
	$isotropic_file = "galdiffuse/isotropic_iem_v02_P6_V11_DIFFUSE.txt";
	$galactic_name = "gal_v02";
	$extra_galactic_name = "eg_v02";
	}

elsif($irf_code == 3){
	$irf_prefix = "P7_V6_";
        $gll_file = "diffuseModels/gal_2yearp7v6_v0.fits";
	$gll_file = "galdiffuse/gal_2yearp7v6_v0.fits";
# This hardwires in "source" for now. Other alternative would be "clean".
	$isotropic_file = "galdiffuse/iso_p7v6source.txt";
	$galactic_name = "gal_2yearp7v6_v0";
	$extra_galactic_name = "iso_p7v6source";
	}

elsif($irf_code == 4){
	$irf_prefix = "P7REP_";
	$gll_file = "galdiffuse/gll_iem_v05.fits";

# This hardwires in "source" for now. Other alternative would be "clean".
	$isotropic_file = "galdiffuse/iso_source_v05.txt";
	$galactic_name = "gll_iem_v05";
	$extra_galactic_name = "iso_source_v05";
	}

# IRF code is "overloaded" to also specify the fixed versions
# of the background files
elsif($irf_code == 5){
	$irf_prefix = "P7REP_";
	$gll_file = "galdiffuse/gll_iem_v05_rev1.fit";

# This hardwires in "source" for now. Other alternative would be "clean".
	$isotropic_file = "galdiffuse/iso_source_v05_rev1.txt";
	$galactic_name = "gll_iem_v05_rev1";
	$extra_galactic_name = "iso_source_v05";
	}

elsif($irf_code == 6){
	$irf_prefix = "P8R2_";

	$gll_file = "galdiffuse/gll_iem_v06.fits";

	$isotropic_file = "galdiffuse/iso_P8R2_SOURCE_V6_v06.txt";

	$galactic_name = "gll_iem_v06";
	$extra_galactic_name = "iso_source_v06";

	$MODEL_MAKER_FLAG = "oldNames=True";
}

elsif($irf_code == 7){
	$irf_prefix = "P8R3_";

	$MODEL_MAKER = "makeFL8Yxml";

# Galactic interstellar
	$gll_file = "galdiffuse/gll_iem_v06.fits";


# isotropic
	$isotropic_file = "galdiffuse/iso_P8R3_SOURCE_V2.txt";

	$galactic_name = "gll_iem_v06";
	$extra_galactic_name = "iso_P8R3_SOURCE_V2";

#	$MODEL_MAKER_FLAG = "oldNames=True";
# don't rename extended sources but use 3FGL names
# requires hacked version of make3fglxml
#	$MODEL_MAKER_FLAG = "E3FGL=True, oldNames=True";
	$MODEL_MAKER_FLAG = "E2CAT=True, oldNames=True";

}

elsif($irf_code == 8){
	$irf_prefix = "P8R3_";

# Galactic interstellar
	$gll_file = "galdiffuse/gll_iem_v07.fits";

# isotropic
	$isotropic_file = "galdiffuse/iso_P8R3_SOURCE_V2_v1.txt";

	$galactic_name = "gll_iem_v07";
	$extra_galactic_name = "iso_P8R3_SOURCE_V2_v1";

	$MODEL_MAKER_FLAG = "E2CAT=True, oldNames=True";
}

elsif($irf_code == 9){
        $irf_prefix = "P8R3_";

# Galactic interstellar
         $gll_file = "galdiffuse/gll_iem_v07.fits";

# isotropic
          $isotropic_file = "galdiffuse/iso_P8R3_SOURCE_V3_v1.txt";

          $galactic_name = "gll_iem_v07";
          $extra_galactic_name = "iso_P8R3_SOURCE_V3_v1";

                                         #       $MODEL_MAKER_FLAG = "oldNames=True";
                                         # don't rename extended sources but use 3FGL names
                                         #       $MODEL_MAKER_FLAG = "E3FGL=True, oldNames=True";
          $MODEL_MAKER_FLAG = "E2CAT=True, oldNames=True";

}


else{
	print "Got an invalid IRF code ($irf_code) somehow!\n";
} 


$prefix = get_string("Give any prefix for output file name(s) [optional]", $DEFAULT_PREFIX);

# if phase range is selected add this to the prefix (it's a feature, not a bug!)
if($start_phase > 0.0 || $stop_phase < 1.0){
	$prefix = $prefix . $start_phase . "-" . $stop_phase . "_";

	print "Prefix set to $prefix\n";
}


$save_products = get_yn("Save products (e.g. exposure cubes) for each step?", $DEFAULT_SAVE_PRODUCTS);

$auto_model = get_yn("Auto-generate model files from FGL catalog?", $DEFAULT_AUTO_MODEL);


##############################
# Check if MY_FERMI_DIR has been set if doing auto-model option
##############################

if($auto_model){

if (!defined($status = $ENV{'MY_FERMI_DIR'} ))
{
    	print "\n Error: You need to set the environment variable MY_FERMI_DIR.\n";
        print " This is where the FGL catalog FITS file and the python script \"$MODEL_MAKER.py\" are located.\n\n";
    	exit(0);
	}
}



# Official Fermi directory
$fermi_dir = $ENV{'FERMI_DIR'};

# My own Fermi directory. FGL catalog and makexFGLxml.py must be here.
$my_fermi_dir = $ENV{'MY_FERMI_DIR'};

########################################################
# Prompt for inputs
#######################


if($auto_model){
	qprint ("auto model\n");
	$catalog = get_string("Give source catalog file name", $DEFAULT_CATALOG);
}
else{
	qprint ("NOT auto model\n");
	$default_model_file = get_file("Give XML model file", $DEFAULT_MODEL_FILE);
}


$use_binned = get_yn("Use binned likelihood?", $DEFAULT_USE_BINNED);

$ts_threshold = get_value("Give minimum TS value for including fluxes in output", $DEFAULT_TS_THRESHOLD, "INDEF", "INDEF");


# File containing source names and parameters and open the file
$source_list = get_file_open("Give source parameter file", *SOURCE_LIST, $DEFAULT_SOURCE_LIST);

# Get names of files which contain lists of the
# photon files and a single spacecraft file
# that must cover the time range of all the photon files

$ft1_file_list = get_file("Give photon file name list", $DEFAULT_FT1_LIST);


# Look for spacecraft files in the current directory.
opendir(DIR, ".");
@files = grep(/_SC00\.fits$/,readdir(DIR));
closedir(DIR);
$files = @files;
#print "no. of SC files found = $files\n";

# If there's exactly one spacecraft file in the current directory suggest that as
# default

if($files == 1){
	$ft2_file = get_file("Give spacecraft (FT2) file name", $files[0]);
}
else{
        $ft2_file = get_file("Give spacecraft (FT2) file name", $DEFAULT_FT2);

}

# spacecraft files should be sorted but they weren't sometimes
# uncomment the lines below if the problem recurs
#$command = "fsort  infile=\"$ft2_file\[1\]\" columns=START method=\"heap\"";
#	doit($command);



$roi = get_value("Give aperture radius (degrees)", $DEFAULT_ROI, 0, 180);

# source radius is 10 degrees bigger than ROI
$srcrad = $roi + 10;

$emin_to_use = get_value("Give Emin (MeV)", $DEFAULT_EMIN, 0, 1E10);


$emax_to_use = get_value("Give Emax (MeV)", $DEFAULT_EMAX, $emin_to_use, 1E10);


$zenith_limit = get_value("Give Zenith limit (degrees)", $DEFAULT_ZENITH_LIMIT, 0, "INDEF");



# Pass 6 and Pass 7 cause gtselect to require different key words.

if($irf_code == 1 || $irf_code == 2){

	$event_class_max = 10; #big value deliberate

	$event_class_min = get_value("Give EVENT_CLASS minimum (use 0 for simulated data)", 
				$DEFAULT_EVENT_CLASS_MIN, 0, 3);

	if($event_class_min == 3 || $event_class_min == 0){
		$irfs = $irf_prefix."DIFFUSE";
		}
	elsif($event_class_min == 2){
		$irfs = $irf_prefix."SOURCE";
		}
	elsif($event_class_min == 1){
		$irfs = $irf_prefix."TRANSIENT";
		}
	else{
		print "event class minimum not allowed\n";
	}

}
elsif($irf_code == 3){
	qprint ("For P7 only SOURCE event class is currently provided\n");
	$event_class_min = 0; 
	$event_class_max = 0;
	$event_class = 2;
	$irfs = "P7SOURCE_V6";
}
elsif($irf_code == 4){
	qprint ("For P7REP only SOURCE event class is currently provided\n");
	$event_class_min = 0; 
	$event_class_max = 0;
	$event_class = 2;
	$irfs = "P7REP_SOURCE_V15";
}

elsif($irf_code == 5){
	qprint ("For P7REP only SOURCE event class is currently provided\n");
	$event_class_min = 0; 
	$event_class_max = 0;
	$event_class = 2;
	$irfs = "CALDB";
}

elsif($irf_code == 6|| $irf_code == 7 || $irf_code == 8 || $irf_code == 9){
	qprint ("For P8 only SOURCE event class is currently provided\n");
	$event_class = 128; 
# clean would be 256

# evtype: 3 = front and back, 1 = front only, 2 = back only.
	$event_type = 3;

	$irfs = "CALDB";
	$irfs_specific = "P8R3_SOURCE_V3";
}

else{
	print "Shouldn't be here - got bad value of irf_code ($irf_code)\n";
}


$rock_limit = get_value("Give rock angle limit (degrees)", $DEFAULT_ROCK, -90, 90);

$bore_limit = get_value("Give Bore limit (degrees)", $DEFAULT_BORE, 0, 360);

$bin_size = get_value("Give bin size (phase; -ve = number of bins)", $DEFAULT_BIN_SIZE, -1000, 1);
if($bin_size < 0.0){
	$bin_size = 1.0/abs($bin_size);
	print "Converted to bin size of $bin_size\n";
}

# Convert to seconds for use in script
#$bin_size = $bin_size * 86400;

$period = get_value("Give Period (days)", $DEFAULT_PERIOD, 0, 1E10);
# but gtpphase wants seconds
$period = $period * 86400;

$t0 = get_value("Give time that defines phase = 0 (MJD)", $DEFAULT_T0, 0, 1E10);

# Convert time zero to MET

$t0 = mjd2met($t0);

# Length of time between phase bins
#$dead_time_length = (1 - $bin_size) * $period;


$sun_minimum = get_value("Give minimum solar distance (degrees)", $DEFAULT_SUN_MINIMUM, 0, 180);

# PFILES directory

# https://heasarc.gsfc.nasa.gov/lheasoft/scripting.html
$pfiles_tmp = "/tmp/pfiles_like_bphase_$$.tmp";
$command = "mkdir $pfiles_tmp";
	doit($command);
	$ENV{'PFILES'} = $pfiles_tmp.":".$ENV{'PFILES'};
#	print "pfiles: $ENV{'PFILES'}\n";
# I don't think this should be necessary. But if don't do
# reverts to using default pfiles directory	
$command = "cp $ENV{'FERMI_DIR'}/syspfiles/*.par $pfiles_tmp/.";
	doit($command);
$command = "cp $ENV{'HEADAS'}/syspfiles/f*.par $pfiles_tmp/.";
	doit($command);



############################################
# Start of source loop
###############

qprint ("Reading source list:\n");


#while(chomp($line = <SOURCE_LIST>)){
while($line = <SOURCE_LIST>){
    chomp($line);

    print "$line\n";

# remove any commas
$line =~ s/\,/ /g;

# Need to have at least source name, RA and dec
# in input file

#@fields = split(/\s+/, $line);
@fields = split(' ', $line);
        $length = @fields;

$source = $fields[0];
$ra = $fields[1];
$dec = $fields[2];


# Check we have at least 3 columns in the source list file
if($length < 3){
	print "*** Error in input file!\n";
	print "The source list file ($source_list) contains a line with less than 3 entries:\n";
	print "$line\n";
	print "Each line in this file must contain (at least) source name, RA, dec.\n";
	exit(0);
}


# We have $default_model_file, $model_file, $model_file_in, $model_file_out
# which is probably too many and we should be able to simplify!
# Check if a model file was specified in the entry
$model_given = $FALSE;
if($length > 3) {
	$model_file = $fields[3];
	open(TEMP, "$model_file") || 
               die "Can't open the specified model file, $model_file\n";
	close(TEMP);
	$model_given = $TRUE;	
}
else {
	$model_file = $default_model_file;
}


# change any plus sign in source name to "p"
# to avoid FTOOLS problems (may be better way to work
# around this!
# Use source name as is so can be used with xml file directly
$src = $source;
$src =~ s/\+/p/g;
# need to escape plus signs to do matching
$source_find = $source;
$source_find =~ s/\+/\\\+/g;


# Create output file to contain light curve
$outfile = $prefix . $src."_like_phase.dat";

open(OUTFILE, ">$outfile") || 
               die "Can't open output results file: $outfile\n";


# Make names of temporary files
$eventfile0 = make_temp("eventfile0.fits");
$eventfile1 = make_temp("eventfile1.fits");
$eventfile2 = make_temp("eventfile2.fits");
$eventfile3 = make_temp("eventfile3.fits");
$expcube = make_temp("expcube.fits");
$ltcube = make_temp("ltcube.fits");

$expmap = make_temp("expmap.fits");
$counts_map = make_temp("counts_map.fits");
$counts_cube = make_temp("counts_cube.fits");

$source_map = make_temp("source_map.fits");

$results_file = make_temp("results.dat");
$specfile = make_temp("counts_spectra.fits");

#@tidy_file_list = ($eventfile1, $eventfile2, $expcube, $expmap, $results_file, $specfile);

@tidy_file_list = ($eventfile1, $eventfile2, $expcube, $ltcube, $expmap, $counts_map, $counts_cube, $source_map,  $results_file, $specfile);

# print "tidy_list = \n @tidy_file_list\n";

# Combine all input photon files together

if($irf_code == 6 || $irf_code == 7 || $irf_code == 8  || $irf_code == 9){
$command = "gtselect chatter=$chatter zmax=$zenith_limit emin=$emin_to_use emax=$emax_to_use infile=\"\@$ft1_file_list\" outfile=\'$eventfile0\' ra=$ra dec=$dec rad=$roi evclass=$event_class evtype=$event_type tmin=0 tmax=0";
    doit($command);
}
elsif($irf_code >= 3 && $irf_code <= 5){
$command = "gtselect chatter=$chatter zmax=$zenith_limit emin=$emin_to_use emax=$emax_to_use infile=\"\@$ft1_file_list\" outfile=\'$eventfile0\' ra=$ra dec=$dec rad=$roi evclass=$event_class evclsmin=$event_class_min evclsmax=$event_class_max tmin=0 tmax=0";
    doit($command);
}
else{
$command = "gtselect chatter=$chatter zmax=$zenith_limit emin=$emin_to_use emax=$emax_to_use infile=\"\@$ft1_file_list\" outfile=\'$eventfile0\' ra=$ra dec=$dec rad=$roi evclsmin=$event_class_min evclsmax=$event_class_max tmin=0 tmax=0";
    doit($command);
}




if($auto_model && !$model_given){
$py_file = make_temp("pyfile.py");

# make model file (using Tyrel Johnson's script)

# Make temporary python script to run makexFGLxml

	$model_file_in = make_temp("like_bphase_xmlmodel.xml");
	$model_file_out = "edit_like_bphase_xmlmodel.xml";


	open(PY_FILE, ">$py_file");
	print PY_FILE "import sys\n";
	print PY_FILE "sys.path.append(\"$my_fermi_dir\")\n";
	print PY_FILE "from $MODEL_MAKER import *\n";
	print PY_FILE "s1=srcList(\'$my_fermi_dir/$catalog\',\'$eventfile0\',\'$model_file_in\')\n";
# Force extended sources to be treated as point-like.
	print PY_FILE "s1.makeModel(\'$fermi_dir/refdata/fermi/$gll_file\',\'$galactic_name\',\'$fermi_dir/refdata/fermi/$isotropic_file\',\'$extra_galactic_name\',psForce=True, $MODEL_MAKER_FLAG)\n";

	close(PY_FILE);

	$command = "python $py_file";
		doit($command);

# Edit output of makexFGLxml to make everything constant except source of interest

	open(MODEL_FILE_IN, $model_file_in) || 
               die "Can't open input model file: $model_file_in\n";

	open(MODEL_FILE_OUT, ">$model_file_out") || 
               die "Can't open output model file: $model_file_out\n";

#$tempfile = "tempfile.dat";
#debug
#	open(TEMPFILE, ">$tempfile") || 
#               die "Can't open temp file $tempfile\n";


	$MY_SOURCE = $FALSE;

	qprint ("looking for source: $source ");

	while(chomp($line = <MODEL_FILE_IN>)){
#	print "$line\n";
# Check if it's a source name
#		if($line =~ /source name/){
		if($line =~ /source/){
#			print "$line\n";
# Check to see if line matches source name
			if($line =~ /$source_find/){
#				print "and found the source.\n";
				$MY_SOURCE = $TRUE;
			}
			else {
				$MY_SOURCE = $FALSE;
#				print "but didn't find the source.\n";
			}
		}

# Freeze parameters for all sources at cataloged values except
# for 
# This is for demonstration primarily, it is not necessarily recommended.
# You might to have either more or less parameters frozen depending on
# the number of photons in each time bin.
# Time taken to perform likelihood analysis depends on how
# many parameters are free.

#	if($MY_SOURCE){print TEMPFILE "TRUE\n";}
#	else{print TEMPFILE "FALSE\n";}

#	print TEMPFILE "before test\n $line\n";


		if(!$MY_SOURCE){
#debug
#			print "$line\n";
			$line =~ s/parameter free=\"1\"/parameter free=\"0\"/g;
#			print "$line\n";
			}



#	print TEMPFILE "$line\n";


		print MODEL_FILE_OUT "$line\n";
	}
	print MODEL_FILE_OUT "$line\n";

	close(MODEL_FILE_IN);
	close(MODEL_FILE_OUT);

#	close (TEMPFILE);

	$command = "rm $py_file";
	doit($command);


	$command = "rm $model_file_in";
	doit($command);

}
else {
# Use a file which is either the one entered at a prompt, or one read from the source list file
	$model_file_out = $model_file;
}




# The file just created is the "master" combined file. 
# Rerun gtselect on it later to just use portions of the data file.

# get TSTART and TSTOP of master photon file

$command = "fkeypar \"$eventfile0\" TSTART";
    doit($command);

chomp($ph_tstart = `pget fkeypar value`);
qprint ("photon start time = $ph_tstart\n");

$command = "fkeypar \"$eventfile0\" TSTOP";
    doit($command);

chomp($ph_tstop = `pget fkeypar value`);
qprint ("photon stop time = $ph_tstop\n");

$time1 = max($ph_tstart, $start_time);
$time2 = min($ph_tstop, $stop_time);

########################################################
# Start of phase loop 
###########################


#$time1 = max($ph_tstart, $start_time);
#$finish_time = min($ph_tstop, $stop_time);

#$time1_mjd = met2mjd($time1);
#$finish_time_mjd = met2mjd($finish_time);


#$phase1 = 0.0;

$phase1 = $start_phase;


print "start phase = $start_phase, phase1 = $phase1\n";
print "stop phase = $stop_phase\n";
print "bin size = $bin_size\n";


#print "Using light curve start and end times of MET $time1 and $finish_time, = MJD $time1_mjd, $finish_time_mjd \n";

#while (($phase1 + $bin_size) <= 1.0) {
while (($phase1 + $bin_size) <= ($stop_phase + $bin_size/1000.0)) {

$BIN_COUNTER++;

$phase2 = $phase1 + $bin_size;


qprint ("-------------------------------------------------------------\n");
qprint ("Calculating for phase interval $phase1, $phase2\n");

# gtselect - potentially restrict time range if command line arguments were given to specific start/stop time. 
# Call gtselect differently depending on whether we use Pass 6, 7 or 8.
# (Since we don't restrict by time or phase here, maybe this doesn't do anything anymore...)


if($irf_code == 6 || $irf_code == 7 || $irf_code == 8  || $irf_code == 9){
$command = "gtselect chatter=$chatter zmax=$zenith_limit emin=$emin_to_use emax=$emax_to_use infile=\"$eventfile0\" outfile=\'$eventfile1\' ra=$ra dec=$dec rad=$roi evclass=$event_class evtype=$event_type tmin=$time1 tmax=$time2";
    doit($command);
}
elsif($irf_code >= 3 && $irf_code <= 5){
$command = "gtselect chatter=$chatter zmax=$zenith_limit emin=$emin_to_use emax=$emax_to_use infile=\"$eventfile0\" outfile=\'$eventfile1\' ra=$ra dec=$dec rad=$roi evclass=$event_class evclsmin=$event_class_min evclsmax=$event_class_max tmin=$time1 tmax=$time2";
    doit($command);
}
else{
$command = "gtselect chatter=$chatter zmax=$zenith_limit emin=$emin_to_use emax=$emax_to_use infile=\"$eventfile0\" outfile=\'$eventfile1\' ra=$ra dec=$dec rad=$roi evclsmin=$event_class_min evclsmax=$event_class_max tmin=$time1 tmax=$time2";
    doit($command);
}




# Loop through phase steps to make additional GTI selections. 
# Calculate time intervals



$t0_phase = $t0 + $phase1 * $period;
$p1s = ($time1 - $t0_phase)/$period;
$ncyc = POSIX::floor($p1s);
$p1s = $ncyc * $period + $t0_phase;
$p1e = $p1s + $bin_size * $period;



qprint ("phase1 = $phase1\n");
qprint ("time1 = $time1; t0 = $t0; t0_phase = $t0_phase; p1s = $p1s; ncyc = $ncyc; p1e = $p1e\n");
qprint ("(ph_start = $ph_tstart; start_time = $start_time; ph_tstop = $ph_tstop)\n");



# gtmktime - good time intervals - regular selections and first time selection
$command = "gtmktime chatter=$chatter scfile=\"$ft2_file\" filter=\"(STOP>$p1s) && (DATA_QUAL>0) && ABS(ROCK_ANGLE)<$rock_limit  && (LAT_CONFIG==1) && (angsep(RA_ZENITH,DEC_ZENITH,$ra,$dec)\+$roi<$zenith_limit) && (angsep($ra,$dec,RA_SUN,DEC_SUN)>$sun_minimum\+$roi) && (angsep($ra,$dec,RA_SCZ,DEC_SCZ)<$bore_limit)\" roicut=n evfile=\"$eventfile1\" outfile=\"$eventfile2\"";
    doit($command);

#$command = "gtmktime chatter=$chatter scfile=\"$ft2_file\" filter=\"(STOP>$p1s)\" roicut=n evfile=\"$eventfile2\" outfile=\"$eventfile3\"";
#$command = "gtmktime chatter=$chatter scfile=\"$ft2_file\"  roicut=n evfile=\"$eventfile2\" outfile=\"$eventfile3\"";
#    doit($command);
#$command = "mv $eventfile3 $eventfile2";
#	doit($command);


# How many binary cycles are there in the time range?
# (not needed?)
# $nbinary = ($time2 - $time1) / $period;


#while ($p1s < ($time2 + $period)){
while ($p1s < $time2){

$filter_time1 = $p1s;
$filter_time2 = $p1s + $bin_size * $period;
# assign good time intervals with gtmktime
$filter_string = "(STOP<$filter_time1) || ";
my $i = 0;
while (($i < $nslots) && ($filter_time1 < $time2)) {
	$filter_string = $filter_string . "((START>$filter_time1) && (STOP<$filter_time2))";
	$filter_time1 += $period;
	$filter_time2 += $period;
	$filter_string = $filter_string . " || ";
	$i++;
}
	$filter_string = $filter_string . "(START>$filter_time1)";

# Need to add a STOP < filter_time1 (or is it filter_time2) at start of expression


$p1s = $filter_time1;


	$command = "gtmktime chatter=$chatter scfile=\"$ft2_file\" filter=\"$filter_string\" roicut=n evfile=\"$eventfile2\" outfile=\"$eventfile3\"";
#		print "command = $command\n";

   		doit($command);

#exit;


#$ttest = ($p1e - $period);

#	$command = "gtmktime chatter=$chatter scfile=\"$ft2_file\" filter=\"(START>$p1s) || (STOP<($p1e - $period))\" roicut=n evfile=\"$eventfile2\" outfile=\"$eventfile3\"";
#	$command = "gtmktime chatter=$chatter scfile=\"$ft2_file\" filter=\"(START>$p1s) || (STOP<$ttest)\" roicut=n evfile=\"$eventfile2\" outfile=\"$eventfile3\"";

#	$p1s += $period;
#	$p1e += $period;

	$command = "mv $eventfile3 $eventfile2";
		doit($command);

#exit;

}

#exit;


# Check for existence of output file
# Can sometimes get:
# Caught St13runtime_error at the top level: Zero rows returned from FT2 file for this filter:

if(!check_file_exists("$eventfile2")){
	tidy_files(@tidy_file_list);
	$phase1 = $phase2;
	next;
}


if($use_binned){

########################
# Binned analysis
#
# 2. Make a counts map from the event data
# (not actually needed for fitting)
# $command = "gtbin chatter=$chatter evfile=$eventfile2 scfile=NONE outfile=$counts_map algorithm=CMAP 

# 3. Create a 3-D (binned) counts map

# hardcoded for 0.2 degree pixels
$binsz = 0.2;

$nxpix = $roi/$binsz;
$nypix = $roi/$binsz;

$command = "gtbin chatter=$chatter evfile=$eventfile2 scfile=NONE outfile=$counts_cube algorithm=CCUBE binsz=$binsz coordsys=CEL xref=$ra yref=$dec nxpix=$nxpix nypix=$nypix axisrot=0.0 proj=AIT ebinalg=LOG enumbins=37 emin=$emin_to_use emax=$emax_to_use";
	doit($command);

# 6. Compute livetimes and exposure

$command = "gtltcube chatter=$chatter evfile=$eventfile2 scfile=$ft2_file outfile=$ltcube dcostheta=0.025 binsz=1";
	doit($command);


# 7. Compute exposure map

# hardcoded for 0.2 degree pixels
$binsz = 0.2;

$nxpix = ($roi+30.0)/$binsz;
$nypix = ($roi+30.0)/$binsz;

$command = "gtexpcube2 chatter=$chatter infile=$ltcube outfile=$expcube irfs=$irfs_specific cmap=none binsz=$binsz coordsys=CEL xref=$ra yref=$dec nxpix=$nxpix nypix=$nypix axisrot=0.0 proj=AIT ebinalg=LOG enumbins=37 emin=$emin_to_use emax=$emax_to_use";
	doit($command);


# 8. Compute source map

$command = "gtsrcmaps chatter=$chatter scfile=$ft2_file expcube=$ltcube cmap=$counts_cube srcmdl=$model_file_out bexpmap=$expcube outfile=$source_map irfs=CALDB";
	doit($command);

# gtlike - do the actual likelihood analysis
$command = "gtlike results=$results_file specfile=$specfile chatter=$chatter irfs=$irfs expcube=$ltcube srcmdl=$model_file_out statistic=BINNED optimizer=MINUIT evfile=$eventfile2 scfile=$ft2_file expmap=$expmap cmap=$source_map bexpmap=$expcube";
	doit($command);


#exit;
######################
}

else{


########################
# Unbinned analysis
# gtltcube - live time cube
$command = "gtltcube chatter=$chatter evfile=$eventfile2 scfile=$ft2_file outfile=$expcube dcostheta=0.025 binsz=1";
    doit($command);

# gtexpmap - exposure map
$command = "gtexpmap chatter=$chatter evfile=$eventfile2 scfile=$ft2_file expcube=$expcube outfile=$expmap irfs=$irfs srcrad=$srcrad nlong=120 nlat=120 nenergies=20";
    doit($command);

# gtdiffrsp - calculate diffuse response columns. (If using other than DIFFUSE, which is pre-computed, this will take a long time!)
$command = "gtdiffrsp chatter=$chatter evfile=\'$eventfile2\' scfile=$ft2_file srcmdl=$model_file_out irfs=$irfs";
	doit($command);

# gtlike - do the actual likelihood analysis
$command = "gtlike results=$results_file specfile=$specfile chatter=$chatter irfs=$irfs expcube=$expcube srcmdl=$model_file_out statistic=UNBINNED optimizer=MINUIT evfile=$eventfile2 scfile=$ft2_file expmap=$expmap cmap=none bexpmap=none";
	doit($command);
#    $doit_result = doit($command);
########################

}




# Extract flux from results of gtlike
# skip to next time interval if no results file available
unless(open(RESULTS_FILE, "<$results_file")){ 
        print "Can't open gtlike output file: $results_file\n"; 
	print "(Probably gtlike failed. Possible solution is to simplify model file.)\n";
	print "Going to try to continue anyway with next time step.\n";
	tidy_files(@tidy_file_list);
	$phase1 = $phase2;
	next;
}


# Find the line containing the source flux and write to the output light curve file
$MY_SOURCE = $FALSE;
$FOUND_SOURCE = $FALSE;

qprint ("Looking for: $source_find ");


##########################



while(<RESULTS_FILE>){
	chomp($line = $_);

	if($line =~ /$source_find/){
			qprint ("and found the source.\n");
			$MY_SOURCE = $TRUE;
			$FOUND_SOURCE = $TRUE;
		}
		else {
			$MY_SOURCE = $FALSE;
		}

# Extract spectral information and read into array
	$spec_index = 0;
	if($MY_SOURCE){
		until($line =~ /\}/){

#		  print "$line\n";
		  $var = $line;

# remove anything before colon
		  $var =~ s/.*://;
# remove quotes and commas
		  $var =~ s/,//g;
		  $var =~ s/'//g;
# Now remove +/- sign
		$var =~ s/\+\/\-//g;

#	print "var = $var\n";

		if($line =~ /TS value/){
		   $ts_value = $var;
#		   print "got a TS value of: $ts_value\n";			
		}

		if($line =~ /Flux/){
			$out = $var;
		}

		  $spec_pars[$spec_index] = $var;
	      	  $spec_index++;
		  chomp($line = <RESULTS_FILE>);	
		}

#		print "@spec_pars\n";
# write output.

		$phase = ($phase1 + $phase2)/2;
		$phase_err = ($phase2 - $phase1)/2;
		print "Center of current phase bin = $phase\n";

# Exclude measurements with low TS values		
		if($ts_value > $ts_threshold){
				print OUTFILE "$phase $out $phase_err $ts_value @spec_pars\n";
		}
		else{
				print OUTFILE "$phase $out $phase_err $ts_value @spec_pars\n";
#			print "TS value ($ts_value) less than threshold ($ts_threshold)\n";
		}
	}
}




# tidy up files

tidy_files(@tidy_file_list);


if(!$FOUND_SOURCE){
	$command = "rm $eventfile0";
		doit($command);
	die
"ERROR: Couldn't find source $src in the results from gtlike.
The source specified in slist.dat is (probably) not present in the model file.
Check for typos in slist.dat and/or the model file.\n";}




###########################
# End of time loop for one source
###########################

$phase1 = $phase2;

}

print "End of analysis for $source\n";

close(OUTFILE);

# End of source loop

}

$command = "rm $eventfile0";
	doit($command);

print "End of analysis for all sources\n";


# Remove temporary PFILES directory and contents
 $command = "rm -rf $pfiles_tmp";
	doit($command);
#	print "command: $command\n";





############################
# Subroutines
##########






##########################################
# convert MJD to MET
# start and stop times are specified on command line in MJD units
# but science tools use MET

sub mjd2met{
my $value = shift;

# This makes an assumption about UT vs. TT....
$result = ($value - (51910.0 + 7.428703703703703E-4)) * 86400.0;

return($result)
}

##########################################
# convert MET to MJD
# start and stop times are specified on command line in MJD units
# but science tools use MET

sub met2mjd{
my $value = shift;

# This makes an assumption about UT vs. TT....
$result = ($value / 86400.0) + (51910.0 + 7.428703703703703E-4);

return($result)
}


#################################################
# System call, and write to screen what it's doing
# Usage: doit("string to execute");

sub doit{
my $doit_result;
    qprint ("\n$_[0]\n");
    $doit_result = system($_[0]);
#    print "system result = $doit_result\n";
    return $doit_result;
}

###########################
# regular print, if not batch mode
sub qprint{
    if(!$BATCH || $DEBUG){
		print "$_[0]";
	}
}

####################################################
# make a temporary filename to use that includes the
# process ID to (hopefully) make it unique
# Usage: $file_name = make_temp("file_name_base")

sub make_temp{

# The line below created output files on /tmp
# However, many systems only have limited space there.
#return "/tmp/".$$.$_[00];
	return "tmp_".$$.$_[00];
}

#######################################
# Prompt for and get a file name.
# Check that file exists before accepting its name.
# Usage: $file_name = get_file("String to prompt for name"[, "default name"]);
# The default file name is not required.

sub get_file{

my $prompt = $_[0];
my $default = $_[1];

my $waiting = 1;


if($BATCH){
	$file_name = $default;
	if (open(CHECK, $file_name)){ 
		close(CHECK);
    	}
	else {
        	die "Can't open default file: $file_name\n";
   	}

	return $default;
}


while($waiting){

if ($default){
	print "$prompt [default = $default]\: ";
	}
else{
	print "$prompt\: ";
}

chomp($file_name = <STDIN>);

if($file_name eq "" && $default)
       {$file_name = $default; print "(Using default: \"$default\")\n";}
if (open(CHECK, $file_name)){ 
#	print "opened $file_name\n";
	$waiting = 0;
	close(CHECK);
    }
else {
               print "Can't open file: $file_name\n";
   }
}
	return $file_name;
}





#######################################
# Prompt for and get a file name, then open it.
# This should probably be combined with get_file()
# Check that file exists before accepting its name.
# Usage: $file_name = get_file("String to prompt for name", FILE_HANDLE, [, "default name"]);
# The default file name is not required.

sub get_file_open{

my $prompt = $_[0];
my $check = $_[1];
my $default = $_[2];

if($BATCH){
        $file_name = $default;
        open($check, $file_name)||
               die "Can't open file: $file_name\n";
        return $file_name;
}

my $waiting = 1;


while($waiting){

if ($default){
	print "$prompt [default = $default]\: ";
	}
else{
	print "$prompt\: ";
}

chomp($file_name = <STDIN>);

if($file_name eq "" && $default)
       {$file_name = $default; print "(Using default: \"$default\")\n";}
if (open($check, $file_name)){ 
#	print "opened $file_name\n";
	$waiting = 0;
#	close(CHECK);
    }
else {
               print "Can't open file: $file_name\n";
   }
}
return $file_name;
}




############################################################
# Prompt for a Y/N response and return true or false (1 or 0)
# default value is determined by argument
# Usage: $response = get_yn("Prompt", $default);
# $default should be 1, 0, or "INDEF"

sub get_yn{
my $prompt = $_[0];
my $default = $_[1];

if($BATCH){
        return $default;
}


if($default eq 1){$default_yn = "Y";}
elsif($default eq 0){$default_yn = "N";}

my $INDEF = "INDEF";

my $waiting = 1;



while($waiting){

if($default eq $INDEF){
	print "$prompt [y/n]\: ";
}
else{
	if ($default){
		print "$prompt [y/n; default = y]\: ";
	}
	else{
		print "$prompt [y/n; default = n]\: ";
	}
}

chomp($response = substr(uc(<STDIN>),0,1));

#print "response = $response\n";

#if($response eq "" && ($default != $INDEF))
if($response eq "")
       {$value = $default; print "(Using default: \"$default_yn\")\n"; $waiting = 0;}
# Check response format OK
elsif($response eq "Y"){
#	print "got y\n";
	$value = 1; $waiting = 0;
}
elsif($response eq "N"){
#	print "got n\n";
	$value = 0; $waiting = 0;
}
else {
           print "Response ($response) is not \"y\" or \"n\"\n";
   }
}
return $value;


}


######################################
# Get a numerical value within specified range
# and with possible default value
sub get_value{

my $prompt = $_[0];
my $default = $_[1];
my $minimum = $_[2];
my $maximum = $_[3];

#print "default = $default, min = $minimum, max = $maximum\n";

if($BATCH){
	return $default;
}

my $INDEF = "INDEF";

my $waiting = 1;


while($waiting){


if ($default eq $INDEF){
	print "$prompt\: ";
	$use_default = 0;
	}
else{
	print "$prompt [default = $default]\: ";
	$use_default = 1;
}


chomp($value = <STDIN>);

if($value eq "" && ($use_default))
#if($value eq "")
       {$value = $default; print "(Using default: \"$default\")\n";}
# Check data range is OK

if(check_value($value, $minimum, $maximum)){
	 	$waiting = 0;
	}
#if (($value <= $maximum || $maximum eq $INDEF) && 
#	(  $value > $minimum )  || $minimum eq $INDEF     ){ 
#	$waiting = 0;
#    }
else {
               print "Value ($value) is not in the range $minimum to $maximum\n";
   }
}
return $value;
}




############################################################
# Prompt for and get a string.
# No range check.
# Usage: $value = get_string("String to prompt for value", default)
# If a default is not to be specified, call with value set to 
# the string "INDEF"

sub get_string{

$prompt = $_[0];
$default = $_[1];
#$minimum = $_[2];
#$maximum = $_[3];


if($BATCH){
	return $default;
}

$INDEF = "INDEF";

$waiting = 1;

while($waiting){

if ($default){
	print "$prompt [default = $default]\: ";
	}
else{
	print "$prompt\: ";
}

chomp($value = <STDIN>);
#print "value = $value\n";
if($value eq "" && ($value ne $INDEF))
       {$value = $default; print "(Using default: \"$default\")\n";
	$waiting = 0;}
if($value ne ""){$waiting = 0;}
# Check data range is OK
#if (($value <= $maximum || $maximum eq $INDEF) && 
#	(  $value > $minimum )  || $minimum eq $INDEF     ){ 
#	$waiting = 0;
#    }
#else {
#               print "Value ($value) is not in the range $minimum to $maximum\n";
#   }
}
return $value;
}


##########################################
# 
sub check_value{
my $value = shift;
my $minimum = shift;
my $maximum = shift;

my $result = $TRUE;



if($value !~ /[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?/){
	print "$value is not a number!\n";
	$result = $FALSE;
	return $result;
	}
#[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?

if($maximum ne "INDEF"){
	if($value > $maximum){$result = $FALSE;}
	}
if($minimum ne "INDEF"){
	if($value < $minimum){$result = $FALSE;}
	}

return $result;
}




###################################################
# 
# check for the existence of a specified file
# return TRUE if it exists, otherwise FALSE

sub check_file_exists{

my $file_name = shift;

if (open(CHECK, $file_name)){ 
#	print "opened $file_name\n";
	close(CHECK);
	return $TRUE;
    }
else{
	print "File $file_name not found\n";
	return $FALSE;
    }
}

################################################
# tidy_files
# rename or delete temporary files
sub tidy_files{

# Probably want to tidy names changing files
# to remove PID as well as "tmp"

$save_flag = $save_products;

my $command;

print "save_flag = $save_flag\n";

while(@_){
	my $file_name = shift(@_);
	if (-e $file_name){

		if($save_flag){
print "save\n";

			$save_name = $BIN_COUNTER . $prefix . $file_name;
			$save_name =~ s/tmp_//;
			$command = "mv $file_name $save_name" ;
		}
	     	else{
print "delete\n";
			$command = "rm $file_name";
		}
#		print "$command\n";
		doit ($command);
	}
}

}

##########################################

################################################
# tidy_files
# copy a bunch of intermediate files
sub tidy_files2{
my $save_flag = shift;
my @file_list = shift;

print "tidy_list = \n @file_list\n";

if($save_flag){
		save_files(@file_list);
	}
else{
		rm_files(@file_list);
	}
}


##########################################


################################################
# save_files
# copy a bunch of intermediate files
sub save_files{
while(@_){
	my $file_name = shift(@_);
	if (-e $file_name){
		$save_name = $BIN_COUNTER . $prefix . $file_name;
		my $command = "mv $file_name $save_name" ;
		print "SAVE_FILES: $command\n";
		doit ($command);
		}
	}
$BIN_COUNTER++;

}


##########################################

################################################
# rm_files
# delete a bunch of intermediate files
sub rm_files{
while(@_){
	my $file_name = shift(@_);
	if (-e $file_name){
		my $command = "rm $file_name";
		doit ($command);
		}
	}
}


##########################################
#
# tweak out spectral parameter and its error
# model independent version
sub new_extract_spectral_parameter{
my $source = shift;
my $in_line = shift;
my $name = shift;
my $result = shift;


#print "-------\nin new_extract_spectral_parameter\n";

#print "$name\n$in_line\n";

# Extract Index and error value
	if($source && ($in_line =~ /\'$name\'\:/)){
#		print "found the string!\n";
# parse index and error (yes, this could be rewritten more compactly!)
# but perhaps has more robustiness??
		$result = $in_line;
		$result =~ s/\'$name\'://;
		$result =~ s/\+\/\-//;
		$result =~ s/\'//g;
		$result =~ s/,//;
	}

#print "result = $result\n";
if($result ne "INDEF"){
# print "Got value of $result\n";
}
return $result;
}




##########################################
#
# tweak out spectral parameter and its error
sub extract_spectral_parameter{
my $source = shift;
my $in_line = shift;
my $name = shift;
my $result = shift;


#print "-------\nin extract_spectral_parameter\n";

#print "$name\n$in_line\n";

# Extract Index and error value
	if($source && ($in_line =~ /\'$name\'\:/)){
#		print "found the string!\n";
# parse index and error (yes, this could be rewritten more compactly!)
# but perhaps has more robustiness??
		$result = $in_line;
		$result =~ s/\'$name\'://;
		$result =~ s/\+\/\-//;
		$result =~ s/\'//g;
		$result =~ s/,//;
	}

#print "result = $result\n";
if($result ne "INDEF"){
# print "Got value of $result\n";
}
return $result;
}

##########################################
# Change default values if specified from a file
# Likely to be most useful if running in "batch" mode
sub set_defaults{

open(PARAMETERS, $parameter_filename) ||
     die "Can't open input parameter file: $parameter_filename\n";

#while(chomp($par_line = <PARAMETERS>)){

while($par_line = <PARAMETERS>){
	chomp($par_line);

# ignore comment lines
#print "first character:", substr($par_line, 0, 1), "\n";
	if(substr($par_line, 0, 1) ne "\#"){
		@par_fields = split(' ', $par_line);
# blank line OK
		if(@par_fields == 0){
		}		
		elsif((@par_fields == 1 || @par_fields > 2) && $par_fields[0] ne "prefix"){
			print "Invalid number of fields in parameter file, must have exactly 2\n";
			print "$par_line\n";
		}	

		if(@par_fields == 2){
# Yes, this could be done more elegantly!
# convert to upper case
			$par_fields[0] = lc($par_fields[0]);
			if($par_fields[0] eq "ft2") {$DEFAULT_FT2 = $par_fields[1];}
                        elsif($par_fields[0] eq "use_binned") {$DEFAULT_USE_BINNED = $par_fields[1];}
			elsif($par_fields[0] eq "bin_size") {$DEFAULT_BIN_SIZE = $par_fields[1];}
			elsif($par_fields[0] eq "period") {$DEFAULT_PERIOD = $par_fields[1];}
			elsif($par_fields[0] eq "t0") {$DEFAULT_T0 = $par_fields[1];}
			elsif($par_fields[0] eq "roi") {$DEFAULT_ROI = $par_fields[1];}
			elsif($par_fields[0] eq "zenith_limit") {$DEFAULT_ZENITH_LIMIT = $par_fields[1];}
			elsif($par_fields[0] eq "rock") {$DEFAULT_ROCK = $par_fields[1];}
			elsif($par_fields[0] eq "bore") {$DEFAULT_BORE = $par_fields[1];}
			elsif($par_fields[0] eq "event_class_min") {$DEFAULT_EVENT_CLASS_MIN = $par_fields[1];}
			elsif($par_fields[0] eq "source_list") {$DEFAULT_SOURCE_LIST = $par_fields[1];}
			elsif($par_fields[0] eq "ft1_list") {$DEFAULT_FT1_LIST = $par_fields[1];}
			elsif($par_fields[0] eq "catalog") {$DEFAULT_CATALOG = $par_fields[1];}
			elsif($par_fields[0] eq "irf_code") {$DEFAULT_IRF_CODE = $par_fields[1];}
			elsif($par_fields[0] eq "auto_model") {$DEFAULT_AUTO_MODEL = $par_fields[1];}
			elsif($par_fields[0] eq "model_file") {$DEFAULT_MODEL_FILE = $par_fields[1];}
			elsif($par_fields[0] eq "ts_threshold") {$DEFAULT_TS_THRESHOLD = $par_fields[1];}
			elsif($par_fields[0] eq "sun_minimum") {$DEFAULT_SUN_MINIMUM = $par_fields[1];}
			elsif($par_fields[0] eq "emin") {$DEFAULT_EMIN = $par_fields[1];}
			elsif($par_fields[0] eq "emax") {$DEFAULT_EMAX = $par_fields[1];}
			elsif($par_fields[0] eq "prefix") {$DEFAULT_PREFIX = $par_fields[1];}
			elsif($par_fields[0] eq "chatter") {$DEFAULT_CHATTER = $par_fields[1];}
			else{print "Unrecognized parameter: $par_fields[0]\n";}
		}
#		else{
#			print "Hmmm.... how did I get here?!\n";
#		}
	}	

}
close(PARAMETERS);


}

##########################################
# Write current file defaults to a parameter file
# This can then be edited (if desired) then used as
# future input.

sub write_defaults{

open(PARAMETERS, ">$output_parameter_filename") ||
     die "Can't open output parameter file: $output_parameter_filename\n";


print PARAMETERS
"\# Default parameters of the \"like_bphase\" binary phase resolved aperture photometry code
\# can be set with a file like this.
\# Use the syntax such as
\# \$ like_bphase -file like_bphase.par
\# This can be most useful when running in \"batch\" mode where parameters
\# are not prompted for.
\# e.g.
\# \$ like_bphase -batch -file like_bphase.par
use_binned $DEFAULT_USE_BINNED
ft2 $DEFAULT_FT2
bin_size $DEFAULT_BIN_SIZE
period $DEFAULT_PERIOD
t0 $DEFAULT_T0
roi $DEFAULT_ROI
zenith_limit $DEFAULT_ZENITH_LIMIT
rock $DEFAULT_ROCK
bore $DEFAULT_BORE
event_class_min $DEFAULT_EVENT_CLASS_MIN
source_list $DEFAULT_SOURCE_LIST
ft1_list $DEFAULT_FT1_LIST
catalog $DEFAULT_CATALOG
irf_code $DEFAULT_IRF_CODE
auto_model $DEFAULT_AUTO_MODEL
model_file $DEFAULT_MODEL_FILE
ts_threshold $DEFAULT_TS_THRESHOLD
sun_minimum $DEFAULT_SUN_MINIMUM
emin $DEFAULT_EMIN
emax $DEFAULT_EMAX
prefix $DEFAULT_PREFIX
chatter $DEFAULT_CHATTER\n";

close(PARAMETERS);

exit;

}