File Coverage

blib/lib/Astro/XSPEC/TableModel.pm

Criterion	Covered	Total	%
statement	22	24	91.6
branch			n/a
condition			n/a
subroutine	8	8	100.0
pod			n/a
total	30	32	93.7

# --8<--8<--8<--8<--
#
# Copyright (C) 2008 Smithsonian Astrophysical Observatory
#
# This file is part of Astro::XSPEC::TableModel
#
# Astro::XSPEC::TableModel is free software: you can redistribute it
# and/or modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see .
#
# -->8-->8-->8-->8--

package Astro::XSPEC::TableModel;

use 5.00800;
use strict;
use warnings;

use base 'Exporter';
our @EXPORT_OK = qw( write_table );

use Carp;

use List::Util qw( max );

use Params::Validate qw( :all );

use Astro::FITS::CFITSIO qw( :constants );
use Astro::FITS::CFITSIO::CheckStatus;
use Astro::FITS::Header::Item;
use Astro::FITS::Header::CFITSIO;

our $VERSION = '0.01';

# this rule is too broad
## no critic ( ProhibitAccessOfPrivateData )

sub write_table {

    my @fields;

    my %opt = validate_with(
        params => \@_,
        spec   => {
            output => { type => SCALAR },
            model  => 1,
            units  => 1,
            additive => { type => BOOLEAN,  default => '0' },
            redshift => { type => BOOLEAN,  default => '0' },
            ipars    => { type => ARRAYREF },
            apars    => { type => ARRAYREF, default => [] },
            energy   => { type => ARRAYREF },
            keywords => { type => ARRAYREF, default => [] },
        },
        normalize_keys => sub { lc $_[0] },
    );

    tie my $err, 'Astro::FITS::CFITSIO::CheckStatus';

    my $fits =
      Astro::FITS::CFITSIO::create_file( "!$opt{output}",
        $err = "Error creating $opt{output}:" );

    # -------------------------------------
    # create PRIMARY HDU
    $fits->create_img( 8, 0, 0, $err = 'Error creating primary HDU: ' );
    _write_cards(
        $fits,
        $opt{keywords},
        [
            STRING =>
              CONTENT => 'MODEL',
            'spectrum file contains time intervals and event'
        ],
        [ STRING => MODLNAME => $opt{model}, 'Model name' ],
        [ STRING => MODLUNIT => $opt{units}, 'Model units' ],
        [
            LOGICAL =>
              REDSHIFT => $opt{redshift},
            'If true then redshift will be included as a parameter'
        ],
        [
            LOGICAL =>
              ADDMODEL => $opt{additive},
            'If true then this is an additive table model'
        ],
        [ STRING => HDUCLASS => 'OGIP', 'format conforms to OGIP standard' ],
        [
            STRING =>
              HDUCLAS1 => 'XSPEC TABLE MODEL',
            'model spectra for XSPEC'
        ],
        [ STRING => HDUVERS1 => '1.0.0', 'version of format' ] );

    # -------------------------------------
    # create PARAMETERS HDU

    @fields = (
        { type => 'NAME',     form => '12A', unit => ' ' },
        { type => 'METHOD',   form => 'J',   unit => ' ' },
        { type => 'INITIAL',  form => 'E',   unit => ' ' },
        { type => 'DELTA',    form => 'E',   unit => ' ' },
        { type => 'MINIMUM',  form => 'E',   unit => ' ' },
        { type => 'BOTTOM',   form => 'E',   unit => ' ' },
        { type => 'TOP',      form => 'E',   unit => ' ' },
        { type => 'MAXIMUM',  form => 'E',   unit => ' ' },
        { type => 'NUMBVALS', form => 'J',   unit => ' ' },
        { type => 'VALUE',    form => 'E',   unit => ' ' },
    );
    my %field = ();
    for my $col ( 1 .. @fields )
    {
        my $f = $fields[ $col - 1 ];
        $f->{col} = $col;
        $field{ $f->{type} } = $f;
    }

    my @ipars = _validate_pars( \%field, $opt{ipars}, 'interpolated' );
    my @apars = _validate_pars( \%field, $opt{apars}, 'additional' );

    # update form for VALUE to equal maximum NUMBVALS
    my $nvals = max map { scalar @{ $_->{VALUE} } } @ipars;
    $field{VALUE}{form} = $nvals . $field{VALUE}{form};

    _create_tbl( $fits, 'PARAMETERS', @fields );

    _write_cards(
        $fits,
        [],
        [ STRING => HDUCLASS => 'OGIP', 'format conforms to OGIP standard' ],
        [
            STRING =>
              HDUCLAS1 => 'XSPEC TABLE MODEL',
            'model spectra for XSPEC'
        ],
        [
            STRING =>
              HDUCLAS2 => 'PARAMETERS',
            'extension containing parameter info'
        ],
        [ STRING => HDUVERS1 => '1.0.0', 'version of format' ],
        [
            INT =>
              NINTPARM => scalar @ipars,
            'Number of interpolation parameters'
        ],
        [ INT => NADDPARM => scalar @apars, 'Number of additional parameters' ],
    );

    # -------------------------------------
    # write the parameters

    my $row = 0;
    for my $par ( @ipars, @apars )
    {
        $row++;

        while ( my ( $key, $value ) = each %{$par} )
        {
            my $field = $field{$key};

            if ( $field->{type} eq 'VALUE' )
            {
                my $numbvals = @$value;
                $fits->write_col_dbl( $field->{col}, $row, 1, $numbvals, $value,
                    $err );
                $fits->write_col_dbl( $field{NUMBVALS}{col},
                    $row, 1, 1, $numbvals, $err );
            }
            else
            {
                my $datatype = $field->{form} =~ /A/ ? TSTRING : TDOUBLE;
                $fits->write_col( $datatype, $field->{col}, $row, 1, 1, $value, $err );
            }
        }
    }

    # -------------------------------------
    # create the ENERGIES HDU

    @fields = (
        { type => 'ENERG_LO', form => 'E', unit => ' ' },
        { type => 'ENERG_HI', form => 'E', unit => ' ' },
    );

    _create_tbl( $fits, 'ENERGIES', @fields );

    _write_cards(
        $fits,
        [],
        [ STRING => HDUCLASS => 'OGIP', 'format conforms to OGIP standard' ],
        [
            STRING =>
              HDUCLAS1 => 'XSPEC TABLE MODEL',
            'model spectra for XSPEC'
        ],
        [
            STRING =>
              HDUCLAS2 => 'ENERGIES',
            'extension containing energy bin info'
        ],
        [ STRING => HDUVERS1 => '1.0.0', 'version of format' ],
    );

    # ick.
    my @energy  = @{ $opt{energy} };
    my $nenergy = @energy - 1;
    $fits->write_col_dbl( 1, 1, 1, $nenergy, \@energy, $err );
    shift @energy;
    $fits->write_col_dbl( 2, 1, 1, $nenergy, \@energy, $err );

    # -------------------------------------
    # finally, create the SPECTRA HDU

    @fields = (
        { type => 'PARAMVAL', form => @ipars . 'E',   unit => ' ' },
        { type => 'INTPSPEC', form => $nenergy . 'E', unit => $opt{units} },
    );

    push @fields,
      { type => "ADDSP$_", form => $nenergy . 'E', unit => $opt{units} }
      for 1 .. @apars;

    _create_tbl( $fits, 'SPECTRA', @fields );

    _write_cards(
        $fits,
        [],
        [ STRING => HDUCLASS => 'OGIP', 'format conforms to OGIP standard' ],
        [
            STRING =>
              HDUCLAS1 => 'XSPEC TABLE MODEL',
            'model spectra for XSPEC'
        ],
        [
            STRING =>
              HDUCLAS2 => 'ENERGIES',
            'extension containing energy bin info'
        ],
        [ STRING => HDUVERS1 => '1.0.0', 'version of format' ],
    );

    return $fits;
}

# --------------------------------
sub _create_tbl {

    my ( $fits, $extname, @fields ) = @_;

    tie my $err, 'Astro::FITS::CFITSIO::CheckStatus';

    $fits->create_tbl( BINARY_TBL,
        0,
        scalar @fields,
        [ map { $_->{type} } @fields ],
        [ map { $_->{form} } @fields ],
        [ map { $_->{unit} } @fields ],
        $extname,
        $err = "Error creating $extname HDU: "
    );

    return;
}

# --------------------------------
#  write FITS cards to the CHU.
#  Each card is an array with elements [ type, keyword, value, comment ]
#  where type is a string recognized by Astro::FITS::Header.

sub _write_cards {
    my ( $fptr, $items, @cards ) = @_;

    my $hdr = Astro::FITS::Header::CFITSIO->new( fitsID => $fptr );

    $hdr->insert( -1, $_ ) foreach @$items;

    $hdr->insert(
        -1,
        Astro::FITS::Header::Item->new(
            Type    => $_->[0],
            Keyword => $_->[1],
            Value   => $_->[2],
            Comment => $_->[3] ) ) foreach @cards;
    $hdr->writehdr( fitsID => $fptr );

    return;
}

# --------------------------------
#  validate parameter structures
sub _validate_pars {

    my ( $fields, $pars, $type ) = @_;

    # rewrite parameter hash with normalized keywords
    my @pars;

    my %keys = map { $_ => 1 } keys %{$fields};

    # exclude NUMBVALS,  as that's generated automatically
    my @delkeys = qw( NUMBVALS );

    # "additional" parameters don't need this
    push @delkeys, qw( METHOD VALUE )
      if $type eq 'additional';

    delete @keys{@delkeys};
    my @keys = keys %keys;

    my $row = 0;
    for my $par (@$pars)
    {
        $row++;

        # normalize keywords
        my %pars = map { uc $_ => $par->{$_} } keys %$par;

        my $name = $pars{NAME} || $row;

        # any missing?
        my @missing = grep { !defined $pars{$_} } @keys;
        croak( "$type parameter $name: missing attributes: ",
            join( ', ', @missing ), "\n" )
          if @missing;

        # any extra?  really meant to catch errors in "additional" parameters
        my @extra = grep { defined $pars{$_} } @delkeys;
        croak( "$type parameter $name: illegal extra attributes: ",
            join( ', ', @extra ), "\n" )
          if @extra;

        push @pars, \%pars;
    }

    return @pars;
}

1;

__END__