Index: trunk/tools/listcmf
===================================================================
--- trunk/tools/listcmf	(revision 33938)
+++ trunk/tools/listcmf	(revision 33938)
@@ -0,0 +1,367 @@
+#!/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
+#
+# (The mana command files may be re-used for different sets of inputs that have the same OUTPUT.FORMAT and
+# sets of extensions
+#
+# To read the data into mana
+#
+#       datafile r.list
+#       input r.mana
+#
+#       datafile i.list
+#       input i.mana
+#
+#       set i_minus_r = PSF_INST_MAG_i - PSF_INST_MAG_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 = EXT_INST_MAG_DEV_i - EXT_INST_MAG_DEV_r
+#
+# similarly
+#
+#       set sersic_i_minus_r = EXT_INST_MAG_SERSIC_i - EXT_INST_MAG_SERSIC_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, $save_temps, $verbose);
+
+GetOptions(
+            'opihi-read-file|r=s'   => \$opihi_read_file,
+            'append'                => \$append,
+	    'save-temps'            => \$save_temps,
+            '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;
+
+die "usage: listcmf <cmf_filename> [<listfile> <ext>]\n" unless defined($cmf);
+
+# 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 ($exp_time, $zero_point, $filter);
+my %extensions;
+# read the keywords and list of extensions from the header
+{
+    my $header_data = `ftlist $cmf kh include=EXPTIME,FPA.ZP,FPA.FILTERID`;
+
+    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;
+            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 eq 'hdr') {
+                # primary header
+                printf STDERR "HDU %2d %4s  %-12s\n", $hdu_num, $ext, "$type $cols" if $verbose;
+            }
+        } elsif ($line =~ /=/) {
+            # keyword line
+            my ($key, $value);
+            if ($line =~ /HIERARCH/) {
+                (undef, $key, undef, $value) = split " ", $line;
+                # print STDERR "$key $value\n" if $verbose;
+            } else {
+                ($key, undef, $value) = split " ", $line;
+                # print STDERR "$key $value\n" if $verbose;
+            }
+            if ($key eq 'EXPTIME') {
+                $exp_time = $value;
+            } elsif ($key eq 'FPA.FILTERID') {
+                $filter = substr $value, 1, 1;
+            } elsif ($key eq 'FPA.ZP') {
+                $zero_point = $value;
+            }
+        }
+    }
+    die "failed to find filter in header" unless defined $filter;
+    die "failed to find exposure time in header" unless defined $exp_time;
+    die "failed to find zero point in header" unless defined $zero_point;
+    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 1";
+    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";
+    }
+}
+
+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";
+    print OUT "read IPP_IDET 1";
+    my $n = 1;
+    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";
+        }
+    }
+}