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

###################################
# Robin Corbet (UMBC), 2010-08-09 #
###################################

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

$VERSION = "2.1";


$MODEL_MAKER = "make4FGLxml";
$MODEL_MAKER_FLAG = "";


$help = "
like_lc.pl

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

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

** Note, from version 2.1 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. 

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    4.0137901E+01   6.1228100E+01
	i.e. the 4FGL catalog name must be given 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, MJD),  flux (cts/cm^2/s), flux error, time bin half width (days), 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.
** The output format was changed in version 1.7 to enable any spectral model to be output.
** One difference is that generally the normalization and its error will be the first two parameters.
** Thus the columns are shifted by two entries compared to the output in previous versions.


Other command line options:

-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:
2.1     Use either binned or unbinned likelihood, with binned as default, change prompt order
2.02    Filter on DATA_QUAL > 0 rather than DATA_QUAL==1
2.01	Catalog update to V27
2.0     Catalog update to V22, updated IRFs to P8R3_SOURCE_V3_v1
1.91    Updated default catalog to gll_psc_v20.fit
1.9	Defaults to 4FGL models, and updated galactic diffuse 
1.8	P8R3 default data, FL8Y, and updated isotropic diffuse background (2018-02)
1.72	bug fixes for batch mode
1.71	Create and delete temporary PFILES directory to avoid collisions.
1.7	Output full set of parameters for any spectral model. Includes e.g. normalization which
	was missing before, along with Npred, ROI distance.
1.62	Added -prefix command line option
1.61	Updated help information, added -version option
1.6	Use latest version of make3FGLxml with backwards compatibility for XML files
1.5	Allow use of PASS8 data
1.47	Use 3FGL catalog as default
1.46	Allow reading from/writing to parameter file.
1.45	Include spectral output for log parabola and broken power law models. Make default maximum energy 300 GeV,
	allow optional prefix for output light curve name.
1.44	Command line options to specify start and stop times of light curve.
        Minor changes to avoid spurious error messages.
1.43	Support for updated background files for P7REP and use of \"CALDB\" for IRFs
1.42	Changed name of output xml files to decrease name collision possibility (2014 February)
1.41	Typo correction in prompts (October 2013).
1.4	Support for Pass 7 reprocessed data. (September 2013).
1.35	Change to checks on existence of environment variables.
1.34	Print time step at each loop iteration.
1.33	Changed (yet again) location of diffuse files. (2012 April).
1.32	Change to diffuse models used.
1.31	Bug fix for Pass 6 data.
1.3     Support for Pass 7 data (default) as well as Pass 6. Use 2FGL catalog and make2FGLxml.py.
1.2	Fix for bug causing mangled output if logParabola model was used
1.1	Support for selecting either P6_V3 or P6_V11 IRFs

";



##########################################
# 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 minimum start time
$start_time = 0;
# default maximum stop time
$stop_time = 9e99;

# minimal output from tools
$DEFAULT_CHATTER = 0;

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

$DEFAULT_BIN_SIZE = 5;
$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 DR2 + new IRF

$DEFAULT_SUN_MINIMUM = 5.0;

$DEFAULT_AUTO_MODEL = $TRUE;

$DEFAULT_USE_BINNED = $TRUE;


# 4FGL DR2
$DEFAULT_CATALOG = "gll_psc_v27.fit";


$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_lc.par";

$catalog = $DEFAULT_CATALOG;


$DEFAULT_FT2 = "lat_spacecraft_merged.fits";

print "like_lc.pl: V$VERSION (R. Corbet)\n";
# print "*** Output format changed in V1.7 and later! \"like_lc.pl -help\" for details ***\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] =~ /-18m/) {
#		print "Using 18 month catalog instead of 1FGL\n";
#		$catalog = "gll_psc18month_uw11b.fits";
#		}
	if($ARGV[$argnum] =~ /-help|-h|--help/) {
		print $help;
		exit;
	}
	elsif($ARGV[$argnum] =~ /-version/) {
		print "like_lc.pl: Version $VERSION\n";
		exit;
	}


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


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

	elsif($ARGV[$argnum] =~ /-file/) {
		if($argnum+1 < $numArgs){
			$parameter_filename = $ARGV[$argnum+1];
			print "parameter_filename = $parameter_filename\n";
			set_defaults();
		}
		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] =~ /-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] =~ /-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] =~ /-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] =~ /-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);
}


$bin_size = get_value("Give bin size (days)", $DEFAULT_BIN_SIZE, 0, 1e10);

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

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

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

$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_";

# assume we're in conda environment 
#         $pweight = "pweight_ore";


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







$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($ENV{'MY_FERMI_DIR'} !~/\S/ )
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 catalog FITS file (e.g. 4FGL) 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. xFGL 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);
}


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

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


$roi = get_value("Give data region 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);



$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_likelc_$$.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_lc.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");
$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);

# 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_lc_xmlmodel.xml");
	$model_file_out = "edit_like_lc_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 make2FGLxml 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";



	$MY_SOURCE = $FALSE;

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

	while(chomp($line = <MODEL_FILE_IN>)){

# Check if it's a source name
		if($line =~ /source/){
# 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){
#debug
#			print "$line\n";
			$line =~ s/parameter free=\"1\"/parameter free=\"0\"/g;
#			print "$line\n";
			}


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

	close(MODEL_FILE_IN);
	close(MODEL_FILE_OUT);

	$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");


########################################################
# Start of time loop 
###########################


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

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

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

while ($time1 < $finish_time) {

$time2 = $time1 + $bin_size;

$time1_mjd = met2mjd($time1);
$time2_mjd = met2mjd($time2);
qprint ("-------------------------------------------------------------\n");
qprint ("Calculating for time interval MET $time1, $time2; MJD $time1_mjd, $time2_mjd\n");

# gtselect - reselect time range. Call gtselect differently depending on whether we use Pass 6 or 7.


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



# gtmktime - good time intervals
$command = "gtmktime chatter=$chatter scfile=\"$ft2_file\" filter=\"(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);


# 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);
	$time1 = $time2;
	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);
	$time1 = $time2;
	next;
}

# The check on the existence of a results file seems to take
# care of crashes of gtlike
#if($doit_result){
#        print "Failed call to gtlike\n"; 
#	tidy_files(@tidy_file_list);
#	$time1 = $time2;
#	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;

		if($line =~ /TS value/){
		   $ts_value = $var;
		}

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

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

# write output.
# write output.
# Convert MET to MJD
		$time = ($time1 + $time2)/2;
		$time_met = $time;
		$time = ($time/86400) + $MJDREF;
		$terr = ($time2 - $time1)/(86400 * 2);
		qprint ("Center of current time bin = MJD $time (MET $time_met)\n");


# Exclude measurements with low TS values		
		if($ts_value > $ts_threshold){
				print OUTFILE "$time $out $terr $ts_value @spec_pars\n";
		}
		else{
			qprint ("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 $source 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
###########################

$time1 = $time2;

}

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){
#    if(!$BATCH){
		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
# delete a bunch if intermediate files
sub tidy_files{
while(@_){
	my $file_name = shift(@_);
	if (-e $file_name){
		my $command = "rm $file_name";
		doit ($command);
		}
	}
}


##########################################
#
# 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/)){
	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);

#print "par_line = $par_line\n";
# ignore comment lines
#print "first character:", substr($par_line, 0, 1), "\n";
	if(substr($par_line, 0, 1) ne "\#"){
#print "doing split\n";
		@par_fields = split(' ', $par_line);
#print "did split\n";
# 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 "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";
#		}
	}	

}
#print "about to close parameters\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
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;

}