#!/usr/bin/env perl

# listcmf
#
# A program to read cmf files for skycells and produce text files that may be used 
# to read the data into opihi programs (mana, dvo) using the opihi command 'read'
#
# Depends on the HEASOFT program ftlist which must be in the path
#
# one line is produced per source (IPP_IDET)
#
# If "exended source analysis" (xsrc) and "extended source fits" (xfits) extensions exist
# in the cmf file selected parameters from those tables are included
# Sources that do not have extended source or models data have their corresponding columns set to NULL
# so when read in mana they get set to NAN
# 
# the option 
#   --opihi_read_file (mana filename) 
# is supplied a file that contains the command to read the vectors written to the specified filename
#
# The filter is appended to the names of the columns (except for IPP_IDET) so that list files
# made from the set of outputs from a multifilter staticskyRun may be read in sequence to yield
# a set of vectors with multi-filter information for each object
#
# For example
#   listcmf my.r.band.cmf r.list --opihi_read_file r.mana
#   listcmf my.i.band.cmf i.list --opihi-read-file i.mana
#
# To read the data into mana
#
#       datafile r.list
#       input r.mana
#
#   The vectors are populated and the variables $ZPT_r  and  $EXPTIME_r are set
#
#       datafile i.list
#       input i.mana
#
#       set i_minus_r = $ZPT_i + PSF_INST_MAG_i + 2.5 * log($EXPTIME_i) - ($ZPT_r + PSF_INST_MAG_r + 2.5 * log($EXPTIME_r))
#
# The model fits columns names are appended with a tag for the model name as well as the filter
# So to compare the de Vaucouleurs model magnitude in two filters
#
#       set dev_i_minus_r = $ZPT_i + EXT_INST_MAG_DEV_i + 2.5 * log($EXPTIME_i) - ($ZPT_r + EXT_INST_MAG_DEV_r + 2.5 * log($EXPTIME_r))
#
#
#       set sersic_i_minus_r = $ZPT_i + EXT_INST_MAG_SERSIC_i + 2.5 * log($EXPTIME_i) - ($ZPT_r + EXT_INST_MAG_SERSIC_r + 2.5 * log($EXPTIME_r))

use strict;
use warnings;

use Getopt::Long qw( GetOptions :config auto_help auto_version gnu_getopt );
use Pod::Usage qw( pod2usage );
use File::Basename;

#  make sure HEASOSFT is intialized
die "HEASOFT not initalized try running heainit\n" if !defined $ENV{HEADAS};


my ( $cmf, $out, $opihi_read_file, $append, $help, $verbose);

my ($exp_time, $filter, $zero_point);

GetOptions(
            'opihi-read-file|r=s'   => \$opihi_read_file,
            'append'                => \$append,
            'zero-point|z=s'        => \$zero_point,
            'exptime|e=s'           => \$exp_time,
            'filter|e=s'            => \$filter,
            'help|h'                => \$help,
            'verbose|v'             => \$verbose,
) or pod2usage( 2 );

pod2usage( -msg => "Unknown option: @ARGV", -exitval => 2 ) if @ARGV > 2;

$cmf = $ARGV[0];
$out = $ARGV[1];

# first read the psf extension
my $input = $cmf;

if ($help or !defined $cmf) {
    print STDERR "usage: listcmf <cmf_filename> [<listfile>] [--opihi-read-file <filename>]\n";
    die "usage: listcmf <cmf_filename> [<listfile>] [--append] [--opihi-read-file <filename>] [--zero-point <zp>] [--exptime <exptime>] [--filter <filter>]\n";
}

# output file handle
my $oh;
if (defined $out and $out ne "-") {
    my $to_open = $append ? ">>$out" : ">$out";
    open $oh, $to_open or die "failed to open output list file: $out\n";
} else {
    $oh = *STDOUT;
}


my %extensions;
# read the keywords and list of extensions from the header
{
    my $header_data = `ftlist $cmf kh include=EXPTIME,FPA.ZP,FPA.FILTERID,FILTER,EXPOTIME`;

    die "failed to extract header information form $cmf\n" if !$header_data;

    my @lines = split "\n", $header_data;
    # require at least a header and psf extension and the 3 keywords
    die "failed to extract header information from $cmf\n" if scalar @lines < 11;

    foreach my $line (@lines) {
        next if !$line;
        next if $line =~ /------/;
        if ($line =~ /HDU/) {
            my (undef, $hdu_num, $hdu_name, $type, $cols, undef, undef, $rows) = split " ", $line;

            if (! ($hdu_name =~ /SkyChip/) ) {
#                print STDERR "unexpected extension name found: $hdu_name\n";
#                die "expected SkyChip.*\n";
            }
            my (undef, $ext) = split '\.', $hdu_name;
            $ext = $hdu_name if !defined $ext;
            if ($type eq 'BinTable') {
                printf STDERR "HDU %2d %4s  %-12s %4d cols %7d rows\n", $hdu_num, $ext, $type, $cols, $rows if $verbose;
                # save extension offset (number - 1) in extensions hash
                $extensions{$ext} = $hdu_num - 1;
            } elsif ($ext and $ext eq 'hdr') {
                # primary header
                printf STDERR "HDU %2d %4s  %-12s\n", $hdu_num, $ext, "$type $cols" if $verbose;
            }
        } elsif ($line =~ /=/) {
            # keyword line
            $line =~ s/=//; 
            my ($key, $value);
            if ($line =~ /HIERARCH/) {
                (undef, $key, $value) = split " ", $line;
                # print STDERR "$key $value\n" if $verbose;
            } else {
                ($key, $value) = split " ", $line;
                # print STDERR "$key $value\n" if $verbose;
            }
            if ($key eq 'EXPTIME') {
                $exp_time = $value unless ((lc($value) =~ /nan/) or $exp_time);
            } elsif ($key eq 'EXPOTIME') {
                $exp_time = $value unless $exp_time;
            } elsif ($key eq 'FPA.FILTERID') {
                $filter = substr $value, 1, 1 unless $filter;
            } elsif ($key eq 'FILTER') {
                $filter = substr $value, 1, 1 unless $filter;
                if (!$filter) {
                    $filter = undef;
#                    print STDERR "Null value for filter setting to r\n";
#                    $filter = "r";
                }
            } elsif ($key eq 'FPA.ZP') {
                $zero_point = $value unless $zero_point;
            }
        }
    }
    die "failed to find filter in header" unless defined $filter;
    die "failed to find exposure time in header" unless defined $exp_time;
    if (!defined $zero_point) {
        print STDERR "failed to find zero point in header using 25\n";
        $zero_point = 25.0;
    }
    printf STDERR "$cmf: %2d extensions filter: $filter zpt: $zero_point exposure_time: $exp_time\n", scalar (keys %extensions) if $verbose;
}

my $psf_ext = $extensions{psf};
die "no psf extension found in $cmf\n" unless $psf_ext;

print STDERR "psf extension offset is $psf_ext\n" if $verbose;

# get the psf extensions columns
$input = "$cmf+$psf_ext";
my @ids;
my @lines;
my @psf_cols;

read_extension($input, 'psf', \@ids, \@lines, \@psf_cols, undef, "PADDING");

printf STDERR "read %d cols %d rows from psf extension\n", scalar @psf_cols, scalar @lines if $verbose;

my @xsrc_cols;
my $xsrc_ext = $extensions{xsrc};
if ($xsrc_ext) {
    $input = "$cmf+$xsrc_ext";
    read_extension($input, 'xsrc', \@ids, \@lines, \@xsrc_cols, undef,
        'PROF_SB,PROF_FLUX,PROF_FILL');

    my $total_cols = scalar @psf_cols + scalar @xsrc_cols;
    my $n_psf = scalar @psf_cols;
    extend_incomplete_lines(\@lines, $total_cols, $n_psf);
} else {
    print STDERR "no xsrc extension found in $cmf\n" unless $xsrc_ext;
}

my @xfit_cols;
my $xfit_ext = $extensions{xfit};
if ($xfit_ext) {
    my $xfit_includes = 'IPP_IDET,X_EXT,Y_EXT,X_EXT_SIG,Y_EXT_SIG,EXT_INST_MAG,EXT_INST_MAG_SIG,NPARAMS,MODEL_TYPE,EXT_WIDTH_MAJ,EXT_WIDTH_MIN,EXT_THETA,EXT_PAR_07';

    $input = "$cmf+$xfit_ext";
    my @ext_cols;
    read_extension($input, 'xfit', \@ids, \@lines, \@ext_cols, $xfit_includes, 
        undef);

    my @model_types;
    for (my $i = 0; $i < scalar @lines; $i++) {
        next unless $lines[$i] =~ /PS_MODEL/;
        my @words = split " ", $lines[$i];
        for (my $j = 0; $j < scalar @words; $j++) {
            next unless $words[$j] =~ /PS_MODEL/;
            my (undef, undef, $type) = split /\_/, $words[$j];
            push @model_types, $type;
        }
        last;
    }
    die "no model types found\n" if !scalar @model_types;
    # create new column names one for each model type
    foreach my $model_type (@model_types) {
        foreach my $col (@ext_cols) {
            push @xfit_cols, "${col}_${model_type}";
        }
    }
    my $ncols_previous = scalar @psf_cols + scalar @xsrc_cols;
    my $total_cols = $ncols_previous + scalar @xfit_cols;
    extend_incomplete_lines(\@lines, $total_cols, $ncols_previous);
} else {
    print STDERR "no xfit extension found in $cmf\n" unless $xfit_ext;
}

my $base = basename($cmf) ;
{
    print $oh "# cmf: $base       FILTER: $filter       ZPT: $zero_point       EXPTIME: $exp_time\n";
    print $oh "#\n";
    print $oh "# ";

    print $oh "IPP_IDET ";
    for (my $i = 1; $i < scalar @psf_cols; $i++) {
        print $oh " $psf_cols[$i]_$filter";
    }
    for (my $i = 0; $i < scalar @xsrc_cols; $i++) {
        print $oh " $xsrc_cols[$i]_$filter"
    }
    for (my $i = 0; $i < scalar @xfit_cols; $i++) {
        print $oh " $xfit_cols[$i]_$filter";
    }
    print $oh "\n";
}

foreach my $line (@lines) {
    print $oh "$line\n";
}
close $oh;

if ($opihi_read_file) {
    open OUT, ">$opihi_read_file" or die "failed to open $opihi_read_file\n";
    printf OUT '$EXPTIME_%s = %f' .  "\n", $filter, $exp_time;
    printf OUT '$ZPT_%s = %f' . "\n", $filter, $zero_point;
    printf OUT '$FILTER = %s' . "\n", $filter;
    print OUT "read IPP_IDET 1";
    my $n = 2;
    for (my $i = 1; $i < scalar @psf_cols; $i++) {
        print OUT " $psf_cols[$i]_$filter ",  $n++;
    }
    for (my $i = 0; $i < scalar @xsrc_cols; $i++) {
        print OUT " $xsrc_cols[$i]_$filter ",  $n++;
    }
    for (my $i = 0; $i < scalar @xfit_cols; $i++) {
        print OUT " $xfit_cols[$i]_$filter ",  $n++;
    }
    print OUT "\n";
    close OUT;
}


sub read_extension {
    my ($input, $ext, $r_ids, $r_lines, $r_cols, $include_list, $exclude_list) = @_;

    die "can't have both include list and exclude list\n" if defined $include_list and defined $exclude_list;

    my $col_list = `ftlist $input c`;
    die "failed to read column list\n" unless $col_list;

    my @lines = split "\n", $col_list;
    die "failed to read column list\n" unless scalar @lines;

    my $exclude_str = "";
    my %exclude;
    if ($exclude_list) {
        $exclude_str = '\'[col ';
        my @cols_to_skip = split '\,', $exclude_list;
        foreach my $col (@cols_to_skip) {
            $exclude{$col} = 1;
            $exclude_str .= "-$col;";
        }
        $exclude_str .= ']\'';
    }
    my %include;
    if ($include_list) {
        my @cols_to_keep = split '\,', $include_list;
        foreach my $col (@cols_to_keep) {
            $include{$col} = 1;
        }
    }

    # process the column list including and excluding columns based
    # on the arguments
    foreach my $line (@lines) {
        next if !$line;
        next if $line =~ /HDU/;
        my ($col, $name, $format) = split " ", $line;
        next if $col =~ /\D/;   # not a number
        if (defined($include_list) && ! $include{$name}) {
            print STDERR "excluding $name\n" if $verbose;
            next;
        } elsif ($exclude{$name}) {
            print STDERR "excluding $name\n" if $verbose;
            next;
        } elsif ($ext ne 'psf' and $name eq 'IPP_IDET') {
            print STDERR "excluding $name from $ext column list\n" if $verbose;
            next;
        } else {
            # it's a keeper
            # XXX TODO: check the format and insure that it is not a vector
            # because we currently don't support them
        }
        push @$r_cols, $name;
    }

    my $cmd = "ftlist $input$exclude_str t rownum=no colheader=no";
    $cmd .= " columns=$include_list" if $include_list;
    print STDERR "$cmd\n" if $verbose;

    @lines = `$cmd`;
    if ($ext eq 'psf') {
        # psf extension contains our list of ids
        # XXX: we are assuming here that the ids in a psf extension
        # are sequential
        foreach my $line (@lines) {
            # id is first column
            chomp $line;
            my ($ipp_idet) = split " ", $line;
            push @$r_ids, $ipp_idet;
            push @$r_lines, $line;
        }
    } else {
        my $num_objects = scalar @$r_lines;
        foreach my $line (@lines) {
            chomp $line;
            my ($id, $rest_of_line) = identify_line($line);
            if ($id > $num_objects) {
                die "TILT: id $id is greater than num objects\n";
            }
            # just append the line to this object's line
            $r_lines->[$id] .= "$rest_of_line";
        }
    }
}

# find the id for this line and return it and the rest of the string
# following the id keeping white space
sub identify_line {
    my $line = shift;

    my ($id) = split " ", $line;

    my $start_of_id = index $line, $id;
    my $rest_of_line = substr $line, ($start_of_id + length($id));

    return ($id, $rest_of_line);
}
sub extend_incomplete_lines {
    my $lines = shift;
    my $ntotal = shift;
    my $nprevious = shift;
    for (my $i = 0; $i < scalar @$lines; $i++) {
        my @words = split " ", $lines->[$i];
        my $nwords = scalar @words;
        if ($nwords == $ntotal) {
            # excellent
        } elsif ($nwords == $nprevious) {
            # need to add null words
            for (my $j = $nwords; $j < $ntotal; $j++) {
                $lines->[$i] .= " NULL";
            }
        } else {
            # we are confused
            die "unexpected number of words found $nwords on line $i expected $nprevious or $ntotal\n";
        }
    }
}
