#! /usr/bin/env perl

=head1 SYNOPSIS

removeStationaries [options] FILESET1 FILESET2

=head1 DESCRIPTION

Given two FITS files contain warp catalogs (presumably), remove
detections that lie within some specified distance from each other.
Write both TTI exposures to the datastore.

=cut


use strict;
use warnings;

use Carp;
use Getopt::Long;
use Pod::Usage;
use Cwd;

use FileHandle;
use Math::Trig;
use LWP::Simple;


use Astro::FITS::CFITSIO qw(:longnames :shortnames :constants);

#use PS::MOPS::FITS::SMF;
use PS::MOPS::Constants qw(:all);
use PS::MOPS::Lib qw(:all);
use PS::MOPS::DC::Instance;
use PS::MOPS::DC::Field;
use PS::MOPS::DC::Detection;


use subs qw(
    remove_stationaries
    get_smf
    open_fits
    read_chips
    grab_radecs
    remove_stationaries
    insert_smf
    _chkfs
    _stripquotes
    _comment2surveymode
    _m2tf
    _make_sigmas
);

our $PS1_MAGSIG2SNR = 2.5 / log(10); 
our $TINT32BIT = 41;    # FITS sucks

our @CHIP_NUMS = qw(
     67 57 47 37 27 17
  76 66 56 46 36 26 16 06
  75 65 55 45 35 25 15 05
  74 64 54 44 34 24 14 04
  73 63 53 43 33 23 13 03
  72 62 52 42 32 22 12 02
  71 61 51 41 31 21 11 01
     60 50 40 30 20 10
);

my $instance_name;
my $inst;
my $debug;
my $reject_flags = 0x3888;          # per IPP, see XXX HREF
my $smf1;
my $smf2;
my $sn1 = 4.5;                      # S/N cut for all reject
my $sn2 = 5.0;                      # S/N cut for stars at det limit
my $help;
GetOptions(
    'instance=s' => \$instance_name,
    'smf1=s' => \$smf1,
    'smf2=s' => \$smf2,
    debug => \$debug,
    help => \$help,
) or pod2usage(2);
pod2usage(-verbose => 3) if $help;
my ($fileset1, $fileset2) = @ARGV;
pod2usage(-msg => 'FILESET1 and FILESET2 must be specified') unless ($fileset1 and $fileset2) or ($smf1 and $smf2);


$inst = PS::MOPS::DC::Instance->new(DBNAME => $instance_name);
my $mops_logger = $inst->getLogger();
my $mops_config = $inst->getConfig();

# Some processing globals.
my $stationrad_minv_dpd = $mops_config->{ingest}->{stationrad_minv_dpd};
my $s2n_cut_sta = $mops_config->{ingest}->{s2n_cutoff_stationary} || 5.0;   # go this deep to find stationaries
my $s2n_cut_ing = $mops_config->{ingest}->{s2n_cutoff_ingest} || 6.0;       # ingest everything about this S/N
my $astro_floor_arcsec = $mops_config->{ingest}->{astrometry}->{floor_arcsec};
my $astro_max_arcsec = $mops_config->{ingest}->{astrometry}->{max_arcsec};
my $astro_force_arcsec = $mops_config->{ingest}->{astrometry}->{force_arcsec};


# Fetch stuff from datastore.  Need to get the fileset listing, then 
# find the SMF file in the listing (probably gzipped).
$smf1 = get_smf($fileset1) unless $smf1;
$smf2 = get_smf($fileset2) unless $smf2;


# Grab minimal info to do stationary removal.
my $dets = grab_radecs($smf1, $smf2, $reject_flags);
#my ($keep_map1, $keep_map2) = remove_stationaries($dets);
my $nonstationary_tbl = remove_stationaries($dets);

# Now we have our list of stationaries, so removem.
my $dbh = $inst->dbh();
my $field1;
my $field2;
$field1 = insert_smf($smf1, 'A', $nonstationary_tbl);
$field2 = insert_smf($smf2, 'B', $nonstationary_tbl);
exit;


sub get_smf {
    my ($fileset) = @_;
    my $url;
    my $tmpdir;
    my ($exp_name) = ($fileset =~ /^(o\w+)/);
    $mops_logger->logdie("can't get exp_name from $fileset") unless $exp_name;

    print STDERR "Fetching $fileset...";

    $url = $mops_config->{ingest}->{index_url};
    $url =~ s|/$||;             # strip trailing /
    my $fileset_url = "$url/$fileset/index.txt";

    # We have seen that sometimes IPP times out on requests or returns an empty result.
    my $content;
    my $retries = 0;
    while (1) {
        $content = get($fileset_url);
        last if $content;       # got stuff, YAY
        $mops_logger->logdie("timeout fetching $fileset_url") if $retries > 1;
        $retries++;
        sleep(30);
    }

    my @lines = split /\n/, get($fileset_url) or $mops_logger->logdie("got empty content for fileset $fileset_url");
    my ($smf_entry) = grep /^$exp_name/, @lines;
    my ($smf_file, @dummy) = split /\|/, $smf_entry;
    die "strange smf_file: $smf_entry" unless $smf_file and $smf_file =~ /^$exp_name/;

    # Given a datastore URL, download the fileset index, find the
    # URL to the SMF file, then download the SMF file somewhere.
    # Return a path to the SMF file.  The dir will be auto_cleaned 
    # when the program exits.


    if ($smf_file =~ /gz$/) {
        system("wget --tries=3 -q -O - $url/$fileset/$smf_file | tar zxf -") == 0 or $mops_logger->logdie("wget failed for $url/$fileset/$smf_file");
        my (@smfs) = grep /\.smf$/, `/bin/ls $exp_name.*`;      # XXX LAME
        chomp @smfs;
        $mops_logger->logdie("can't find an SMF file in $smf_file") unless @smfs > 0;
        $mops_logger->logdie("found multiple SMF files in $smf_file") if @smfs > 1;
        $smf_file = $smfs[0];
    }
    elsif ($smf_file =~ /smf$/) {
        system("wget --tries=3 -q -O - $url/$fileset/$smf_file") == 0 or $mops_logger->logdie("wget failed for $url/$fileset/$smf_file");
        my (@smfs) = grep /\.smf$/, `/bin/ls $exp_name.*`;      # XXX LAME
        chomp @smfs;
        $mops_logger->logdie("can't find an SMF file in $smf_file") unless @smfs > 0;
        $mops_logger->logdie("found multiple SMF files in $smf_file") if @smfs > 1;
        $smf_file = $smfs[0];
    }

    print STDERR "done.\n";
    return $smf_file;
}


sub open_fits {
    my ($filename) = @_;
    my $fits_status = 0;
    my $fh = Astro::FITS::CFITSIO::open_file($filename, READONLY(), $fits_status);
    if ($fits_status) {
        fits_report_error('STDERR', $fits_status);
        die;
    }
    return $fh;
}


sub read_chips {
    my ($detsref, $file, $suffix) = @_;
    my $fs = 0;
    my $extname;
    my $exttype;
    my $extver;
    my $comment;
    my $chip;
    my $hdutype;
    my $fits = open_fits($file);

    my $epoch_mjd;
    my $hdunum = 1;
    $fits->movabs_hdu($hdunum, undef, $fs); _chkfs($fs);
    $fits->read_keyword('MJD-OBS', $epoch_mjd, $comment, $fs); _chkfs($fs);

    # Now just loop through HDUs and fetch detections. We need to be able to ID each detection
    # so there are no duplicates, so we will give them
    # $suffix.$extname.$detnum
    printf STDERR $file;
    foreach my $chipnum (@CHIP_NUMS) {
        $extname = "XY$chipnum.psf";
        $extver = 0;
        #$hdutype = ANY_HDU();
        $hdutype = BINARY_TBL();
        $fs = 0;
        $fits->movnam_hdu($hdutype, $extname, $extver, $fs);
        next if $fs == BAD_HDU_NUM();   # not found
        $fits->get_hdu_type($hdutype, $fs);
        next if $hdutype != BINARY_TBL();
        $chip = "XY$chipnum";

        $fits->read_keyword('EXTTYPE', $exttype, $comment, $fs);
        if ($exttype =~ /PS1_V2/) {
            # Got a chip, read it in.
            my $rows;
            $fits->get_num_rows($rows, $fs); _chkfs($fs);
            print STDERR '.';

            my $i;
            my $anynul;

            my @ipp_idet = ();
            my @ra_deg = ();
            my @dec_deg = ();
            my @mags = ();
            my @mag_sigma = ();
            my @flags = ();
            $fits->read_col($TINT32BIT, 1, 1, 1, $rows, 0, \@ipp_idet, $anynul, $fs);       # IPP_IDET
            $fits->read_col(TDOUBLE, 9, 1, 1, $rows, 0, \@mag_sigma, $anynul, $fs);         # PSF_INST_MAG_SIGMA
            $fits->read_col(TFLOAT, 13, 1, 1, $rows, 0, \@mags, $anynul, $fs);              # CAL_PSF_MAG
            $fits->read_col(TDOUBLE, 15, 1, 1, $rows, 0, \@ra_deg, $anynul, $fs);           # RA_PSF
            $fits->read_col(TDOUBLE, 16, 1, 1, $rows, 0, \@dec_deg, $anynul, $fs);          # DEC_PSF
            $fits->read_col($TINT32BIT, 31, 1, 1, $rows, 0, \@flags, $anynul, $fs);         # FLAGS

            # Now save our stuff.
            $i = 0;
            my $s2n;
            my $id;
            while ($i < $rows) {
                $id = "$suffix.$chip." . $ipp_idet[$i];
                if (!(($flags[$i] & 0x3888) or $mag_sigma[$i] <= 0)) {
                    $s2n = $PS1_MAGSIG2SNR / $mag_sigma[$i];
                    if ($s2n >= $s2n_cut_sta) {
                        push @{$detsref}, [ $id, $ra_deg[$i], $dec_deg[$i], $mags[$i] ];
                    }
                }
                $i++;
            }
        }   # if $exttype
    }
    $fits->close_file($fs);
    printf STDERR "\n";

    return $epoch_mjd;
}


sub grab_radecs {
    my ($smf1, $smf2) = @_;

    # Given two FITS file handles for two exposure catalogs, fetch all
    # chips from each device into a single list for each exposure, and
    # extract RA, DEC, MAGS and SIGMAS.
    my @all_dets;
    my $suffix_A = 'A';
    my $suffix_B = 'B';
    my $epochA_mjd = read_chips(\@all_dets, $smf1, $suffix_A);
    my $epochB_mjd = read_chips(\@all_dets, $smf2, $suffix_B);
    my $diff_sec = abs($epochA_mjd - $epochB_mjd) * 86400;

    printf STDERR "TTI is %.1f seconds (%.2f min).\n", $diff_sec, $diff_sec / 60;

    return {
        tti_min => $diff_sec / 60,
        dets => \@all_dets,
    };
}


sub remove_stationaries {
    # Process the current set of detections in $stuff, send it to the cleaners, and
    # remove anything that was scrubbed by astroclean.
    my ($stuff) = @_;
    my $tti_min = $stuff->{tti_min};
    my $dets = $stuff->{dets};
    my ($ac_infh, $ac_infilename);
    my ($id, $ra_deg, $dec_deg, $mag);

    $ac_infilename = 'ac_input';
    $ac_infh = new FileHandle ">$ac_infilename" or die "can't create astroClean input file";

    # Write out the detections from both files to a single file, using a common
    # epoch_mjd so that they appear to astroClean to be in the same field.  Rewrite the ID
    # so that we can back out which file the detections came from.
    my $fake_epoch_mjd = 54242.0;
    my $fake_obscode = 'F51';           # PS1, but doesn't matter
    my $det;

    foreach $det (@{$dets}) {
        ($id, $ra_deg, $dec_deg, $mag) = @{$det};
        print $ac_infh join(' ', $id, $fake_epoch_mjd, $ra_deg, $dec_deg, $mag, $fake_obscode), "\n";
    }
    $ac_infh->close();

    # Convert our min tracklet velocity to a stationary radius for processing.  If the min
    # velocity is not set, use 1 arcsec.
    my $stationrad_thresh_arcsec;
    if (defined($stationrad_minv_dpd)) {
        $stationrad_thresh_arcsec = $stationrad_minv_dpd * 3600 * $tti_min / 1440;
        printf STDERR "Using stationary radius %.1f arcsec (%.1f deg/day).\n", $stationrad_thresh_arcsec, $stationrad_minv_dpd;
    }
    else {
        $stationrad_thresh_arcsec = 1.0;
        printf STDERR "Using stationary radius %.1f arcsec.\n", $stationrad_thresh_arcsec;
    }
    my $stationrad_str = "stationrad " . $stationrad_thresh_arcsec / 3600;      # in degrees

    my $ac_outfilename;
    my $sfn = 'stationaries';
    my $sfh;
    $ac_outfilename = 'ac_outfile';
    #astroclean file $ac_infilename outtype 1 density 1000 Dradius .0010 proxrad $proxrad_thresh_deg clean_file $ac_outfilename
    my $cmd = <<"EOF";
astroclean file $ac_infilename outtype 1 $stationrad_str clean_file $ac_outfilename > ac.out && cut -d' ' -f1 $ac_outfilename > $sfn
EOF
    system($cmd) == 0 or die "command failed: $? : $cmd";

    # Now read the result.  Note that the output of astroclean is the cleaned data,
    # thus the transients -- stationaries are removed.  So our table will be a dict
    # of detections to preserve when we re-read the FITS file.
    $sfh = new FileHandle $sfn or die "can't open $sfn";
    my $line;
    my $nonstationary_tbl = { 'A' => {}, 'B' => {} };
    while ($line = <$sfh>) {
        if ($line =~ /^(A|B)\.(XY\d\d)\.(\d+)$/) {
            $nonstationary_tbl->{$1}->{$2}->{$3} = 1;      # mark A/B, chip, det number
        }
    }
    $sfh->close();                         # just need filename, not fh

    return $nonstationary_tbl;
}


sub insert_smf {
    my ($file, $suffix, $nonstationary_tbl) = @_;

    # Read FITS file, get
    #   epoch_mjd
    #   fpa_id (EXPNAME)
    #   exp_id
    #   boresight RA, DEC
    #   exposure time
    #   filter
    #   posang
    my $fits = open_fits($file);
    my $fs = 0;     # FITS status
    my $comment;
    my $hdutype;
    my $exttype;
    my $extname;
    my $extver;

    my ($epoch_mjd, $ra_deg, $dec_deg, $posang_deg, $exptime_d, $filter_id, $survey_mode, $fpa_id);
    my ($limiting_mag, $obscode);

    # These keywords are in the Primary FITS extension.
    $fits->read_keyword('MJD-OBS', $epoch_mjd, $comment, $fs); _chkfs($fs);
    $fits->read_keyword('RA', $ra_deg, $comment, $fs); _chkfs($fs);
    $fits->read_keyword('DEC', $dec_deg, $comment, $fs); _chkfs($fs);
    $fits->read_keyword('POSANGLE', $posang_deg, $comment, $fs); _chkfs($fs);
    $fits->read_keyword('EXPTIME', $exptime_d, $comment, $fs); _chkfs($fs);
    $exptime_d /= 86400;      # convert seconds to days
    $epoch_mjd += $exptime_d / 2;
    $fits->read_keyword('FILTERID', $filter_id, $comment, $fs); _chkfs($fs);
    $filter_id =~ s/^\W(\w).*/$1/;           

    # These keywords are a little harder to get, since there may be missing chip 
    # tables.  So we have to go hunting for them.
    foreach my $chip_num (@CHIP_NUMS) {
        $extname = "XY$chip_num.hdr";
        $extver = 0;
        $hdutype = IMAGE_HDU();
        $fs = 0;
        $fits->movnam_hdu($hdutype, $extname, $extver, $fs);
        next if $fs == BAD_HDU_NUM();   # not found
        $fits->get_hdu_type($hdutype, $fs);
        next if $hdutype != IMAGE_HDU();

        $fits->read_keyword('CMMTOBS', $survey_mode, $comment, $fs);
        $survey_mode = _comment2surveymode($survey_mode) if $survey_mode;

        $fits->read_keyword('FILENAME', $fpa_id, $comment, $fs);
        # We see some SMFs with 'UNKNOWN' in the filename fields.
        # If so, sigh, try to get from filename.
        if ($fpa_id =~ /UNKNOWN/) {
            ($fpa_id) = ($file =~ /^(o\w+)/);
            die "can't get an FPA_ID from $file" unless $fpa_id;
        }
        else {
            $fpa_id = substr(_stripquotes($fpa_id), 0, -2) if $fpa_id;     # take all of filename except last trailing digits
        }


        last if $survey_mode and $fpa_id;
    }

    # If we have reached this point with undefined $survey_mode and $fpa_id,
    # then there were no good chips in the entire FITS file.  So bail.
    if (!$survey_mode and !$fpa_id) {
        $mops_logger->logwarn("Could not get FILENAME from $file");
        $fits->close_file($fs);
        return;
    }

    $limiting_mag = $mops_config->{site}->{s2n_config}->{$filter_id}->{limiting_mag}
        or die "can't get limiting mag in site config for $filter_id";
    $obscode = $mops_config->{site}->{obscode};
 
    my $field = PS::MOPS::DC::Field->new($inst,
        epoch => $epoch_mjd,
        ra => $ra_deg,
        dec => $dec_deg,
        pa_deg => $posang_deg,
        surveyMode => $survey_mode,
        timeStart => $epoch_mjd - $exptime_d / 2,
        timeStop => $epoch_mjd + $exptime_d / 2,
        status => $FIELD_STATUS_NEW,
        filter => $filter_id,
        limitingMag => $limiting_mag,
        raSigma => 0.0,
        decSigma => 0.0,
        obscode => $obscode,                # XXX from CONFIG
        de => [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        fpaId => $fpa_id,
    );

    my $found = modcf_retrieve($inst, fpaId => $fpa_id);

    if (!$found) {
        my $num_added;
        my $dbh = $inst->dbh();
        print STDERR "Reading $file";
        eval {
            $dbh->begin_work();
            $field->insert();


            # Read detections here.
            my $chip;
            my $chipnum;
            my @dets;

            my $dets_seen;
            my $dets_used;

            foreach my $chipnum (@CHIP_NUMS) {
                $extname = "XY$chipnum.psf";
                $extver = 0;
                $hdutype = BINARY_TBL();
                $fs = 0;
                $fits->movnam_hdu($hdutype, $extname, $extver, $fs);
                next if $fs == BAD_HDU_NUM();
                $chip = "XY$chipnum";

                $fits->read_keyword('EXTTYPE', $exttype, $comment, $fs);
                if ($exttype =~ /PS1_V2/) {
                    $dets_seen = 0;
                    $dets_used = 0;

                    # Got a chip, read it in.
                    my $rows;
                    $fits->get_num_rows($rows, $fs); _chkfs($fs);

                    print STDERR '.' unless $debug;

                    my $i;
                    my $anynul;

                    # FITS column indexes.  Later probably need to get these programmatically
                    # so that we are more robust against FITS schema changes.
                    my $IPP_IDET_COL = 1;   # IPP_IDET
                    my $MAG_SIGMA_COL = 9;  # PSF_INST_MAG_SIGMA
                    my $MAG_COL = 13;       # CAL_PSF_MAG
                    my $RA_DEG_COL = 15;    # RA_PSF
                    my $DEC_DEG_COL = 16;   # DEC_PSF
                    my $FLAGS_COL = 31;
                    my $DET_POSANG_DEG_COL = 6;
                    my $PLTSCALE_COL = 7;
                    my $PSF_X_SIG_COL = 4;
                    my $PSF_Y_SIG_COL = 5;
                    my $MXX_COL = 28;
                    my $MYY_COL = 30;

                    # Need stub arrays for all our columns here.
                    my @ipp_idet = ();
                    my @ra_deg = ();
                    my @ra_sigma_deg = ();
                    my @dec_deg = ();
                    my @dec_sigma_deg = ();
                    my @mags = ();
                    my @mag_sigma = ();
                    my @flags = ();
                    my @det_posang_deg = ();
                    my @pltscale = ();
                    my @psf_x_sig = ();
                    my @psf_y_sig = ();
                    my @mxx = ();
                    my @myy = ();

                    # XXX Hard-coded col indices are for losers; please fix.
                    $fits->read_col($TINT32BIT, $IPP_IDET_COL, 1, 1, $rows, 0, \@ipp_idet, $anynul, $fs);
                    $fits->read_col(TFLOAT, 9, 1, 1, $rows, 0, \@mag_sigma, $anynul, $fs);
                    $fits->read_col(TFLOAT, 13, 1, 1, $rows, 0, \@mags, $anynul, $fs);
                    $fits->read_col(TFLOAT, $PSF_X_SIG_COL, 1, 1, $rows, 0, \@psf_x_sig, $anynul, $fs);
                    $fits->read_col(TFLOAT, $PSF_Y_SIG_COL, 1, 1, $rows, 0, \@psf_y_sig, $anynul, $fs);
                    $fits->read_col(TFLOAT, $DET_POSANG_DEG_COL, 1, 1, $rows, 0, \@det_posang_deg, $anynul, $fs);
                    $fits->read_col(TFLOAT, $PLTSCALE_COL, 1, 1, $rows, 0, \@pltscale, $anynul, $fs);
                    $fits->read_col(TDOUBLE, 15, 1, 1, $rows, 0, \@ra_deg, $anynul, $fs);
                    $fits->read_col(TDOUBLE, 16, 1, 1, $rows, 0, \@dec_deg, $anynul, $fs);
                    $fits->read_col($TINT32BIT, 31, 1, 1, $rows, 0, \@flags, $anynul, $fs);
                    $fits->read_col(TFLOAT, $MXX_COL, 1, 1, $rows, 0, \@mxx, $anynul, $fs);
                    $fits->read_col(TFLOAT, $MYY_COL, 1, 1, $rows, 0, \@myy, $anynul, $fs);


                    # Now save our stuff.
                    $i = 0;
                    my $s2n;
                    my $detnum;
                    my $idet;
                    my $det;

                    my $ra_sigma_deg;
                    my $dec_sigma_deg;
                    my $min_mag = $mops_config->{ingest}->{photometry}->{min_mag} || -99;
                    my $max_mag = $mops_config->{ingest}->{photometry}->{max_mag} || 999;
                    while ($i < $rows) {
                        $dets_seen++;
                        $idet = $ipp_idet[$i];
                        if ($nonstationary_tbl->{$suffix}->{$chip}->{$idet} && $mag_sigma[$i] > 0 and !($flags[$i] & 0x3888)) {
                            $s2n = $PS1_MAGSIG2SNR / $mag_sigma[$i];
                            my $mag = $mags[$i];
                            if ($mag > $min_mag and $mag < $max_mag and ($s2n > $s2n_cut_ing) and $mxx[$i] > 0 and $myy[$i] > 0) {
                                if ($psf_x_sig[$i] ne 'nan' and $psf_y_sig[$i] ne 'nan') {
                                    $dets_used++;

                                    # OK, we're here.  If it's in the keep list, keep it!
                                    ($ra_sigma_deg, $dec_sigma_deg) = _make_sigmas(
                                        $s2n,
                                        $mxx[$i],
                                        $myy[$i],
                                        $psf_x_sig[$i],
                                        $psf_y_sig[$i],
                                        $det_posang_deg[$i],
                                        $pltscale[$i],
                                    );

                                    $det = PS::MOPS::DC::Detection->new($inst,
                                        ra => $ra_deg[$i],
                                        raSigma => $ra_sigma_deg,
                                        dec => $dec_deg[$i],
                                        decSigma => $dec_sigma_deg,
                                        epoch => $epoch_mjd,
                                        mag => $mags[$i],
                                        magSigma => $mag_sigma[$i],
                                        filter => $filter_id,
                                        s2n => $s2n,
                                        isSynthetic => 0,
                                        orient_deg => 0,
                                        length_deg => 0,
                                        obscode => $obscode,
                                        objectName => $MOPS_NONSYNTHETIC_OBJECT_NAME,
                                        detNum => ($idet * 100 + $chipnum),
                                    );
                                    push @dets, $det;
                                }   # nan
                            }   # mag cut
                        }   # flags
                        $i++;
                    }

                    print STDERR "Used $dets_used/$dets_seen detections from $file/$chip.\n" if $debug;
                }
            }
         

            $fits->close_file($fs);
            $field->addDetections(@dets);
        };
        if ($@) {
            $dbh->rollback();
            $mops_logger->logdie($@);       # or just warn?
        }
        else {
            $dbh->commit();
            print STDERR "\nInserted $fpa_id.\n";
        }
    }
    else {
        $mops_logger->warn("already found field for $fpa_id; skipping");
    }
}


sub _make_sigmas {
    my ($s2n, $mxx, $myy, $psf_x_in, $psf_y_in, $posang_deg, $pltscale) = @_;
    my ($ra_sig_deg, $dec_sig_deg);

    if ($astro_force_arcsec) {
        $ra_sig_deg = $astro_force_arcsec / 3600;
        $dec_sig_deg = $astro_force_arcsec / 3600;
    }
    else {
        # Convert PSF X and Y sigmas to RA/DEC.  We have
        # PSF X and Y sigmas in pixels
        # local position angle (PA from north at detection)
        # plate scale, arcsec/pix
        my $psf_x = sqrt($mxx) / $s2n;
        my $psf_y = sqrt($myy) / $s2n;


        my ($det_ra_sig, $det_dec_sig);
        my $cos_posang = cos(deg2rad($posang_deg));
        my $sin_posang = sin(deg2rad($posang_deg));

        $det_ra_sig = $pltscale * sqrt(
            $psf_x * $psf_x * $cos_posang * $cos_posang + 
            $psf_y * $psf_y * $sin_posang * $sin_posang
        );
        $det_dec_sig = $pltscale * sqrt(
            $psf_x * $psf_x * $sin_posang * $sin_posang + 
            $psf_y * $psf_y * $cos_posang * $cos_posang
        );

        # Add floor and Poisson uncertainties in quadrature.
        $ra_sig_deg = sqrt($astro_floor_arcsec * $astro_floor_arcsec + $det_ra_sig * $det_ra_sig);
        $dec_sig_deg = sqrt($astro_floor_arcsec * $astro_floor_arcsec + $det_dec_sig * $det_dec_sig);

        $ra_sig_deg = $astro_max_arcsec if ($ra_sig_deg > $astro_max_arcsec);
        $dec_sig_deg = $astro_max_arcsec if ($dec_sig_deg > $astro_max_arcsec);

        # Correct RA for declination and convert to degrees.
        $ra_sig_deg *= 1 / 3600;
        $dec_sig_deg *= 1 / 3600;
    }

    return ($ra_sig_deg, $dec_sig_deg);
}


sub _chkfs {
    # Util routine to check FITS status.
    my ($status, $msg) = @_;
    if ($status != 0) {
        fits_report_error('STDERR', $status);
        if ($msg) {
            croak $msg;
        }
        else {
            croak;
        }
    }
}


sub _stripquotes {
    my ($str) = @_;
    $str =~ s/['"\s]+//g;
    return $str;
}


sub _comment2surveymode {
    # Convert IPP/OTIS comment string to MOPS survey mode.
    my ($cmt) = @_;
    if ($cmt =~ /3PI/) {
        return "3PI ($cmt)";
    }
    elsif ($cmt =~ /([EM]?SS)/) {
        return "$1 ($cmt)";
    }
    else {
        return 'NA';
    }
}


sub _m2tf {
    # Convert detection morphological params (moments) to detection
    # length and orientation.
    my ($mxx, $myy, $mxy, $pixel_size_deg) = @_;
    my ($g1, $g2, $g3);
    my ($s1, $s2, $ang);

    $g1 = $mxx + $myy;
    $g2 = $mxx - $myy;
    $g3 = sqrt($g2 * $g2 - 4 * $mxy * $mxy);
    $s1 = sqrt(($g1 + $g3) / 2);
    $s2 = sqrt(($g1 - $g3) / 2);

    $s1 = $pixel_size_deg * $s1;

    if ($2 * $mxy == 0 and $g2 == 0) {
        $ang = 0;
    }
    else {
        $ang = rad2deg(atan2(2 * $mxy, $g2));
    }

    return ($s1, $ang);     # length, orientation
}
