#! /usr/bin/env perl
#
#
# Create a response to a MOPS_DETECTABILITY_QUERY
#
#

use strict;
use warnings;

use Sys::Hostname;
use Carp;
use File::Basename;
use File::Copy;
use IPC::Cmd 0.36 qw( can_run run );
use Getopt::Long qw( GetOptions );
use Pod::Usage qw( pod2usage );
use File::Temp qw( tempfile tempdir);

use PS::IPP::PStamp::RequestFile qw( :standard );
use PS::IPP::PStamp::Job qw( :standard );
use PS::IPP::Config qw($PS_EXIT_SUCCESS
                       $PS_EXIT_UNKNOWN_ERROR
                       $PS_EXIT_SYS_ERROR
                       $PS_EXIT_CONFIG_ERROR
                       $PS_EXIT_PROG_ERROR
                       $PS_EXIT_DATA_ERROR
                       $PS_EXIT_TIMEOUT_ERROR
                       metadataLookupStr
                       metadataLookupBool
                       caturi
                       );
use Astro::FITS::CFITSIO qw( :constants );
Astro::FITS::CFITSIO::PerlyUnpacking(1);

#
# Set up
#
my $host = hostname();

my $EXTVER = 1.0;
my $EXTNAME = 'MOPS_DETECTABILITY_RESPONSE';
my ($req_id,$req_name,$product,$need_magic,$missing_tools,$project);
my ($request_file,$output,$workdir,$dbname,$dbserver,$verbose,$save_temps);
GetOptions(
    'input=s'         =>      \$request_file,
    'output=s'        =>      \$output,
    'workdir=s'       =>      \$workdir,
    'dbname=s'        =>      \$dbname,
    'dbserver=s'      =>      \$dbserver,
    'verbose'         =>      \$verbose,
    'save-temps'      =>      \$save_temps,
    ) or pod2usage(2);

pod2usage( -msg => "Unknown option: @ARGV", -exitval => 2 ) if @ARGV;
pod2usage( -msg => "Required options: --input --output --dbname",
	   -exitval => 3,
    ) unless
    defined $request_file and defined $output and defined $workdir and defined $dbname;

my $detect_query_read = can_run('detect_query_read') or (warn "Can't find detect_query_read" and $missing_tools = 1);
my $psphotForced      = can_run('psphotForced') or (warn "Can't find psphotForced" and $missing_tools = 1);
my $dquery_finish     = can_run('dquery_finish.pl') or (warn "Can't find dquery_finish.pl" and $missing_tools = 1);
my $ppCoord           = can_run('ppCoord') or (warn "Can't find ppCoord" and $missing_tools = 1);
if ($missing_tools) {
    warn("Can't find required tools.");
    exit($PS_EXIT_CONFIG_ERROR);
}

unless (defined $verbose) {
    $verbose = 0;
}

my $ipprc = PS::IPP::Config->new();

if (!$dbserver) {
    $dbserver = metadataLookupStr($ipprc->{_siteConfig}, 'PS_DBSERVER');
}
# This is hardcoded in for the moment.
$project = resolve_project($ipprc,"gpc1",$dbname,$dbserver);
my $imagedb = $project->{dbname};
if (!$imagedb) {
    carp("failed to find imagedb for project: $project");
}

#
# Parse input request file using detect_query_read (as it's already written).
#
my $dqr_command = "$detect_query_read --input $request_file";
my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
    run(command => $dqr_command, verbose => $verbose);
unless ($success) {
    # This is a problem, because I'm not sure how we handle a failure to read something.
    # We need to return a $PSTAMP_INVALID_REQUEST, I think, but if we can't read it, 
    # we can't send that response back.
    die("Unable to perform $dqr_command error code: $error_code");
}
my %query = ();
my $Nrows = 0;
{
    my @column_names = ();
    my $section = '';
    foreach my $entry (split /\n/, (join "", @$stdout_buf)) {
	if ($entry =~ /^#/) {
	    @column_names = split /\s+/, $entry;
	    shift(@column_names);  # Dump the hash sign.
	    if ($section eq 'HEADER') {
		$section = 'CONTENT';
	    }
	    else {
		$section = 'HEADER';
	    }
	}
	else {
	    # HEADER: 
	    # QUERY_ID FPA_ID MJD_OBS FILTER OBSCODE STAGE
	    # CONTENT:
	    # ROWNUM RA1_DEG DEC1_DEG RA2_DEG DEC2_DEG MAG
	    my @columns = split /\s+/, $entry;
	    for (my $i = 0; $i <= $#columns; $i++) {
		$query{$section}{$column_names[$i]}[$Nrows] = $columns[$i];
#		print "$section $column_names[$i] $Nrows $columns[$i]\n";
	    }
	    $Nrows++;
	}
    }
}

#
# Identify target images.  This should properly collate targets on a single imfile.
#
my %image_list_hash;
for (my $i = 1; $i < $Nrows; $i++) {
    # This could use the fact that locate_images now accepts position arrays, but
    # I'll save that for after I get the majority of things working.
    my $image_set_tmp  = find_image_set($query{HEADER}{FPA_ID}[0],$query{HEADER}{STAGE}[0],
					$query{HEADER}{MJD_OBS}[0],$query{HEADER}{FILTER}[0],
					$query{CONTENT}{RA1_DEG}[$i],$query{CONTENT}{DEC1_DEG}[$i],
					$query{CONTENT}{ROWNUM}[$i],$verbose);
    unless (%$image_set_tmp) {
	# No images were returned, so create a dummy entry that 
	$image_list_hash{'no_image'}{IMAGE}    = 'no_image';
	$image_list_hash{'no_image'}{PSF}      = 'no_psf';
	$image_list_hash{'no_image'}{MASK}     = 'no_mask';
	$image_list_hash{'no_image'}{WEIGHT}   = 'no_weight';
	$image_list_hash{'no_image'}{CATALOG}  = 'no_catalog';
	$image_list_hash{'no_image'}{CLASS_ID} = 'no_class';
	$image_list_hash{'no_image'}{ERROR}    = $PSTAMP_NO_IMAGE_MATCH;
	push @{ $image_list_hash{'no_image'}{SKY_COORDINATES} }, "$query{CONTENT}{RA1_DEG}[$i] $query{CONTENT}{DEC1_DEG}[$i]";
	push @{ $image_list_hash{'no_image'}{ROWNUM} }, $query{CONTENT}{ROWNUM}[$i];
	next;
    }
#     print "=== $image_set_tmp->{IMAGE}\n    $image_set_tmp->{PSF}\n";
#     print "    $image_set_tmp->{MASK}\n    $image_set_tmp->{WEIGHT}\n";
#     print "    $image_set_tmp->{SKY_COORDINATES}\n    $image_set_tmp->{ROWNUM}\n";

    # This indexes the results for identical images into the same hash.
    $image_list_hash{$image_set_tmp->{IMAGE}}{IMAGE}    = $image_set_tmp->{IMAGE};
    $image_list_hash{$image_set_tmp->{IMAGE}}{PSF}      = $image_set_tmp->{PSF};
    $image_list_hash{$image_set_tmp->{IMAGE}}{MASK}     = $image_set_tmp->{MASK};
    $image_list_hash{$image_set_tmp->{IMAGE}}{WEIGHT}   = $image_set_tmp->{WEIGHT};
    $image_list_hash{$image_set_tmp->{IMAGE}}{CATALOG}  = $image_set_tmp->{CATALOG};
    $image_list_hash{$image_set_tmp->{IMAGE}}{CLASS_ID} = $image_set_tmp->{CLASS_ID};
    $image_list_hash{$image_set_tmp->{IMAGE}}{ERROR}    = $image_set_tmp->{ERROR};
    push @{ $image_list_hash{$image_set_tmp->{IMAGE}}{SKY_COORDINATES} }, $image_set_tmp->{SKY_COORDINATES};
    push @{ $image_list_hash{$image_set_tmp->{IMAGE}}{ROWNUM} }, $image_set_tmp->{ROWNUM};
}

my $i = 0;
foreach my $k (keys %image_list_hash) {
    # If we errored out on finding an image, we need to not try to run psphot here.
    if ($image_list_hash{$k}{ERROR} != 0) {
	next;
    }
    # Write coordinates of the requested targets to a file.
    my ($coordfile,$coordname) = tempfile("${workdir}/detect.coords.$i.XXXX", 
					    UNLINK => !$save_temps);
    my ($targetfile,$targetname) = tempfile("${workdir}/detect.targets.$i.XXXX", 
					    UNLINK => !$save_temps);

    for (my $j = 0; $j <= $#{ $image_list_hash{$k}{SKY_COORDINATES} }; $j++) {
	print $coordfile "$image_list_hash{$k}{SKY_COORDINATES}[$j]\n";
    }
#    print "$k\n";
    # Convert the sky coordinates to image coordinates with ppCoord.
    my $command = "ppCoord -astrom $image_list_hash{$k}{CATALOG} -radec $coordname";
    my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
	run(command => $command, verbose => $verbose);
    unless ($success) {
	my_die("Unable to perform $command. Error_code: $error_code",
	       $query{HEADER}{QUERY_ID}[0],$query{HEADER}{FPA_ID}[0],
	       $query{HEADER}{MJD_OBS}[0],$query{HEADER}{FILTER}[0],$query{HEADER}{OBSCODE}[0],
	       $query{HEADER}{STAGE}[0],$error_code);
    }
    my @response = split /\n/, (join "", @$stdout_buf);
    foreach my $line (@response) {
	my ($r_ra,$r_dec,$trash,$r_x,$r_y,$r_chip) = split /\s+/, $line;
	print $targetfile "$r_x $r_y\n";
	$image_list_hash{$k}{EXTENSION_BASE} = $r_chip;
    }

#     print "psphot $image_list_hash{$k}{PSF}\n";
    # Run psphotForced on the target list.
    my $tmpdir  = tempdir("detect.$i.XXXX", DIR => "${workdir}/", CLEANUP => !$save_temps);
    $image_list_hash{$k}{OUTROOT} = "$tmpdir/detectability.$query{HEADER}{STAGE}[0].$query{HEADER}{FPA_ID}[0]";
    
    my $psphot_cmd = "$psphotForced -psf $image_list_hash{$k}{PSF} ";
    $psphot_cmd .= "-file $image_list_hash{$k}{IMAGE} ";
    $psphot_cmd .= "-mask $image_list_hash{$k}{MASK} ";
    $psphot_cmd .= "-variance $image_list_hash{$k}{WEIGHT} ";
    $psphot_cmd .= "-srctext $targetname $image_list_hash{$k}{OUTROOT}";

    ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
	run(command => $psphot_cmd, verbose => $verbose);
    unless ($success) {
	$image_list_hash{$k}{ERROR} = $PSTAMP_SYSTEM_ERROR;
    }
}

#
# Convert psphot output to response
#
my @rownums = ();
my @out_errors = ();
my @psphot_Npix = ();
my @psphot_Qfact= ();
my @psphot_flux = ();

foreach my $k (keys %image_list_hash) {
    if ($image_list_hash{$k}{ERROR} == 0) {
	my $cmf = "$image_list_hash{$k}{OUTROOT}.$image_list_hash{$k}{CLASS_ID}.cmf";
	
	my ($tmp_Npix,$tmp_Qfact,$tmp_flux) = read_cmf_file($cmf,$image_list_hash{$k}{EXTENSION_BASE});
	
	push @rownums,        @{ $image_list_hash{$k}{ROWNUM} };
	push @out_errors,     (map { $image_list_hash{$k}{ERROR} }  @{ $image_list_hash{$k}{ROWNUM} });
	push @psphot_Npix,    @{ $tmp_Npix };
	push @psphot_Qfact,   @{ $tmp_Qfact };
	push @psphot_flux,    @{ $tmp_flux };
    }
    else {
	push @rownums,        @{ $image_list_hash{$k}{ROWNUM} };
	push @out_errors,     (map { $image_list_hash{$k}{ERROR} }  @{ $image_list_hash{$k}{ROWNUM} });
	push @psphot_Npix,    (map { 0 }  @{ $image_list_hash{$k}{ROWNUM} });
	push @psphot_Qfact,   (map { 0.0 }  @{ $image_list_hash{$k}{ROWNUM} });
	push @psphot_flux,    (map { 0.0 }  @{ $image_list_hash{$k}{ROWNUM} });
    }	
}

write_response_file($output,
		    $query{HEADER}{QUERY_ID}[0],$query{HEADER}{FPA_ID}[0],
		    $query{HEADER}{MJD_OBS}[0],$query{HEADER}{filter}[0],
		    $query{HEADER}{obscode}[0],
		    \@rownums, \@out_errors, \@psphot_Npix, \@psphot_Qfact, \@psphot_flux);
# print "Wrote response file $output\n";
#
# Add to datastore
#
# Files are added to the datastore by dquery_finish.pl
#
# Cleanup
#
# Since everything is written to temporary files, there should be nothing to cleanup.

#
# Utilities
#
sub find_image_set {
    my $FPA_ID = shift;
    my $stage  = lc(shift);
    my $mjd    = shift;
    my $filter = shift;
    my $ra     = shift;
    my $dec    = shift;
    my $index  = shift;
    my $verbose = shift;

    # This is the set of things that we need in order to run psphotForced
    my $option_mask |= 1;
    $option_mask |= $PSTAMP_SELECT_IMAGE;
    $option_mask |= $PSTAMP_SELECT_MASK;
    $option_mask |= $PSTAMP_SELECT_VARIANCE;
    $option_mask |= $PSTAMP_SELECT_PSF;
    my $need_magic = 1;
    my $mjd_min = $mjd;
    my $mjd_max = $mjd + 1;

    # Construct a row list. 
    my @rowList;
    $rowList[0]->{CENTER_X} = $ra;
    $rowList[0]->{CENTER_Y} = $dec;
    $rowList[0]->{ID} = 1;
    $rowList[0]->{STAGE} = $stage;
    $rowList[0]->{COORD_MASK} = 0;

    #    print "$stage\n";
    # Call the PStamp code to find the images that contain the target on the given MJD in the specified filter.
    my @images = locate_images($ipprc,$imagedb,
			        \@rowList,
			       "bycoord",$stage,
			       undef,undef,undef,
			       $option_mask,$need_magic,
			       # $ra,$dec,
			       $mjd_min,$mjd_max,$filter . ".00000",undef,$verbose);  

    my %image_info  = ();
    foreach my $i (@images) {           # Scan over each result
	foreach my $j (@{ $i }) {       # Scan over each image in the result
	    # We only care about an image if it matches the FPA_ID in the request
	    if ($stage eq 'diff') {
		# Diffs match if either exposure name is defined and matches the FPA_ID
		unless ((defined(${ $j }{exp_name_1}) && (${ $j }{exp_name_1} eq $FPA_ID))||
			(defined(${ $j }{exp_name_2}) && (${ $j }{exp_name_2} eq $FPA_ID))||
			(${ $j }{exp_id_1} eq $FPA_ID)||(${ $j }{exp_id_2} eq $FPA_ID)) {
		    next;
		}
	    }
	    elsif ($stage eq 'stack') {
		# Stacks hide the exposure name very well, so we can only match against stage_id
		if (${ $j }{stage_id} ne $FPA_ID) {
		    next;
		}
	    }
	    else {
		# For all the other stages (warp and chip are the ones I've tested), we can simply 
		# directly match the exposure name to the FPA_ID
		unless ((${ $j }{exp_name} eq $FPA_ID)||(${ $j }{exp_id} eq $FPA_ID)) {
		    next;
		}
	    }
	    # Debug prints of all the components of this image
#   	    foreach my $k (keys %{ $j }) {
# 		if ($k eq 'row_index') {
# 		    print "$i $j $k @{${ $j }{$k} }\n";
# 		}
#  		print "$i $j $k ${ $j }{$k}\n";
#  	    }

	    # Check for existance of the images. Drawn mostly from pstampparse.pl
	    my $run_state = ${ $j }{state};
	    my $data_state = ${ $j }{data_state};
	    $data_state = $run_state if $stage eq 'stack';
	    my $fault = 0;
	    if (($run_state eq 'goto_purged') or ($data_state eq 'purged') or
		($run_state eq 'drop') or 
		($run_state eq 'error_cleaned') or 
		($run_state eq 'goto_scrubbed') or ($data_state eq 'scrubbed')) {
		# image is gone and it's not coming back
		$fault = $PSTAMP_GONE;
	    }
	    elsif  (($data_state ne 'full') or ($need_magic and (${ $j }{magicked} < 0))) {
		if (($stage eq 'stack')||($stage eq 'diff')) {
		    # updating stacks and diffs isn't implemented
		    $fault = $PSTAMP_NOT_IMPLEMENTED;
		}
		if ($stage eq 'chip') {
		    my $burntool_state = ${ $j }{burntool_state};
		    if ($burntool_state and (abs($burntool_state) < 14)) {
			$fault = $PSTAMP_NOT_AVAILABLE;
		    }
		}
		
		if ($fault == 0) {
		    # This bombs us out to dqueryparse, which will then flag a job for this run to be updated.
		    my_die_for_update($data_state,$query{HEADER}{STAGE}[0],
				      ${ $j }{stage_id},${ $j }{class_id} || ${ $j }{skycell_id},
				      $need_magic,$imagedb,$PSTAMP_NOT_AVAILABLE);
		}
	    }

	    # This image matches, so we want to save the information into our output structure
	    $image_info{ROWNUM} = $index;
	    $image_info{IMAGE}  = ${ $j }{image};
	    $image_info{PSF}    = ${ $j }{psf};
	    $image_info{MASK}   = ${ $j }{mask};
	    $image_info{WEIGHT} = ${ $j }{weight};
	    $image_info{ERROR}  = $fault;
	    $image_info{SKY_COORDINATES} = "$ra $dec";
	    # To do sky->image coordinate transformations, we need to use the cmf/smf file. If 
	    # an astrom reference (the camera stage smf file) exists, then use that, as we're dealing with
	    # with the chip stage. Otherwise, use the stage-dependent cmf (and set the class_id to fpa).
	    # The EXTENSION_BASE stores the basename of the extension that will be generated by psphotForced.
	    if (exists(${ $j }{astrom})) {
		$image_info{CATALOG} = ${ $j }{astrom};
		$image_info{CLASS_ID} = ${ $j }{class_id};

	    }
	    else {
		$image_info{CATALOG} = ${ $j }{cmf};
		$image_info{CLASS_ID} = 'fpa';

	    }

	}
    }
    return(\%image_info);
}

# Taken largely from detect_query_read
sub read_cmf_file {
    my $cmf_file = shift;
    my $class_id = shift;
    my $extname = $class_id . ".psf";

    my @tmp_Npix = ();
    my @tmp_Qfact = ();
    my @tmp_flux = ();

    my ($detect_F_col,$detect_N_col,$detect_flux_col,$detect_mag_col) = (-1, -1, -1, -1);

    
    my $status = 0;
    my $inFits = Astro::FITS::CFITSIO::open_file( $cmf_file, READONLY, $status ); # Open CMF
    check_fitsio($status);
    $inFits->movnam_hdu(BINARY_TBL, $extname, 0, $status) and check_fitsio($status);

    my $inHeader = $inFits->read_header();

    my $CMFversion = $inHeader->{EXTTYPE};
    $CMFversion =~ s/\'//g;
    $CMFversion =~ s/\s+//g;

    # This is the only data we actually care about
    my $column_defs;
    if ($CMFversion eq 'PS1_V2') {
	$column_defs = [
	    { name => 'PSF_QF',   type => '1E', writetype => TDOUBLE },
	    { name => 'PSF_NPIX', type => '1J', writetype => TLONG },
	    { name => 'PSF_INST_MAG', type => '1E', writetype => TDOUBLE },
	    ];
    }
    elsif ($CMFversion eq 'PS1_DV2') {
	$column_defs = [
	    { name => 'PSF_QF',   type => '1E', writetype => TDOUBLE },
	    { name => 'PSF_NPIX', type => '1J', writetype => TLONG },
	    { name => 'PSF_INST_FLUX', type => '1E', writetype => TDOUBLE },
	    ];
    }

    my @colNames; 
    my @colTypes;
    my @colWriteType;
    my %colData;
    foreach my $colSpec (@$column_defs) {
	push @colNames, $colSpec->{name};
	push @colTypes, $colSpec->{type};
	push @colWriteType, $colSpec->{writetype};
    }
    my $numRows;
    $inFits->get_num_rows($numRows, $status) and check_fitsio($status);

    foreach my $col (@$column_defs) {
	my ($col_num,$col_type,$col_data);
	$inFits->get_colnum(0, $col->{name}, $col_num, $status) and check_fitsio($status);
	if ($status != 0) {
	    $status = 0;
	    next;
	}
	$inFits->get_coltype($col_num, $col_type, undef, undef, $status) and check_fitsio($status);
	$inFits->read_col($col_type, $col_num, 1, 1, $numRows, 0, $col_data, undef, $status) and check_fitsio($status);
	
	if ($col->{name} eq 'PSF_QF') {
	    @tmp_Qfact = @{ $col_data };
	}
	elsif ($col->{name} eq 'PSF_NPIX') {
	    @tmp_Npix = @{ $col_data };
	}
	elsif ($col->{name} eq 'PSF_INST_MAG') {
	    @tmp_flux = map { $_ = 10**(-0.4 * $_) } @{ $col_data };
	}
	elsif ($col->{name} eq 'PSF_INST_FLUX') {
	    @tmp_flux = @{ $col_data };
	}
    }
    $inFits->close_file( $status ) and check_fitsio($status);    
    return(\@tmp_Npix, \@tmp_Qfact, \@tmp_flux);
}

# A lot of this pulled verbatim from Bill's detect_response_create file
sub write_response_file {
    my $outfile = shift;
    my $query_id = shift;
    my $FPA_ID = shift;
    my $MJD_OBS = shift;
    my $filter = shift;
    my $obscode = shift;
    my $rownum_ref = shift;
    my $out_err_ref = shift;
    my $psphot_Npix_ref = shift;
    my $psphot_Qfact_ref = shift;
    my $psphot_flux_ref = shift;
    my $status = 0;

    # Specification of columns to write
    my $columns = [
	# matching rownum from detectability original request
        { name => 'ROWNUM',   type => 'V', writetype => TULONG }, 
	# any errors that occurred during processing
        { name => 'ERROR_CODE',   type => 'V', writetype => TULONG }, 
        # number of pixels used in hypothetical PSF for the query detection
        { name => 'DETECT_N', type => 'V',   writetype => TULONG },
        # detectibility, indicating the fraction of PSF pixels detetable by IPP
        { name => 'DETECT_F', type => 'D',   writetype => TDOUBLE },
	# flux of the target source
	{ name => 'TARGET_FLUX', type => 'D', writetype => TDOUBLE },
	];
    
    # Header translation table
    my $headers = {
	'QUERY_ID' => { name => 'QUERY_ID', type => TSTRING, 
                        comment => 'MOPS Query ID for this batch query' },
	'FPA_ID' => { name => 'FPA_ID',   type => TSTRING, 
		      comment => 'original FPA_ID used at ingest' },
    };

    # Parse the list of columns
    my @colNames;                   # Names of columns
    my @colTypes;                   # Types of columns
    my %colData;                    # Data for each column
    my @colWriteType;                 # type to use to write
    foreach my $colSpec ( @$columns) {
	push @colNames, $colSpec->{name};
	push @colTypes, $colSpec->{type};
	push @colWriteType, $colSpec->{writetype};
	$colData{$colSpec->{name}} = [];
    }

    my $numRows = $#{ $rownum_ref } + 1;
    my $inHeader = { };

    # Hack to force the data to match the detect_response_create formats
    $inHeader->{QUERY_ID}->{value} = $query_id;
    $inHeader->{FPA_ID}->{value} = $FPA_ID;
    $inHeader->{MJD_OBS}->{value} = $MJD_OBS;
    $inHeader->{FILTER}->{value} = $filter;
    $inHeader->{OBSCODE}->{value} = $obscode;
    
    # Fill the table columns with the data, making sure the flux is defined
    for (my $i = 0; $i < $numRows; $i++) {
	push @{$colData{'ROWNUM'}},      ${ $rownum_ref }[$i];
	push @{$colData{'ERROR_CODE'}},  ${ $out_err_ref }[$i];
	push @{$colData{'DETECT_N'}},    ${ $psphot_Npix_ref }[$i];
	push @{$colData{'DETECT_F'}},    ${ $psphot_Qfact_ref }[$i];
	push @{$colData{'TARGET_FLUX'}}, ${ $psphot_flux_ref }[$i];
	unless (defined(${ $colData{'TARGET_FLUX'}}[-1])) {
	    $colData{'TARGET_FLUX'}[-1] = 0.0;
	}
    }

    # Back to detect_response_create:
    my $outFits = Astro::FITS::CFITSIO::create_file( $output, $status );
    check_fitsio( $status );
    $outFits->create_img( 16, 0, undef, $status );
    check_fitsio( $status );
    my $EXTNAME = 'MOPS_DETECTABILITY_RESPONSE';
    $outFits->create_tbl( BINARY_TBL(), $numRows, scalar @colNames, \@colNames, \@colTypes, undef, $EXTNAME, $status);
    check_fitsio( $status );

    # TODO: get the extension version number from somewhere common
    $outFits->write_key( TINT, 'EXTVER', $EXTVER, 'Version of this Extension', $status );
    check_fitsio( $status );

    # Write the Extension keywords
    foreach my $keyword ( keys %$headers ) {
	my $value = $inHeader->{$keyword}->{value}; # Header keyword value
	unless (defined $value) {
	    print "Can't find header keyword $keyword\n";
	    next;
	}
	$value =~ s/\'//g;
	my $name    = $headers->{$keyword}->{name}; # New name
	my $type    = $headers->{$keyword}->{type}; # Type
	my $comment = $headers->{$keyword}->{comment}; # Comment
	$outFits->write_key( $type, $name, $value, $comment, $status );
	check_fitsio( $status );
    }
    for (my $i = 0; $i < scalar @colNames; $i++) {
	my $colName = $colNames[$i];# Column name
	my $writeType = $colWriteType[$i];
	$outFits->write_col( $writeType, $i + 1, 1, 1, $numRows, $colData{$colName}, $status );
	check_fitsio( $status );
    }
    $outFits->close_file( $status );
    return($status);
}

# From Astro::FITS::CFITSIO demo
sub check_fitsio
{
    my $status = shift;         # Status of FITSIO calls

    if ($status != 0) {
        my $msg;                # Message to output
        Astro::FITS::CFITSIO::fits_get_errstatus( $status , $msg );
        die "CFITSIO error: $msg\n";
    }
}

sub my_die {
    my $message = shift;
    my $QUERY_ID = shift;
    my $FPA_ID = shift;
    my $MJD_OBS = shift;
    my $FILTER = shift;
    my $OBSCODE = shift;
    my $STAGE = shift;

    my $exit_code = shift;      # Exit code to add

    $exit_code = $PS_EXIT_PROG_ERROR unless defined $exit_code;

    carp("$message : $QUERY_ID $FPA_ID $MJD_OBS $FILTER $OBSCODE $STAGE");
    exit($exit_code);

}

sub my_die_for_update {
    my $state = shift;
    my $stage = shift;
    my $stage_id = shift;
    my $component = shift;
    my $need_magic = shift;
    my $imagedb = shift;
    my $exit_code = shift;

    print "$state $stage $stage_id $component $need_magic $imagedb\n";
    print STDERR "$state $stage $stage_id $component $need_magic $imagedb\n";
    exit($exit_code);
}
