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


line stmt bran cond sub pod time code
1             # --8<--8<--8<--8<--
2             #
3             # Copyright (C) 2008 Smithsonian Astrophysical Observatory
4             #
5             # This file is part of Astro::XSPEC::TableModel
6             #
7             # Astro::XSPEC::TableModel is free software: you can redistribute it
8             # and/or modify it under the terms of the GNU General Public License
9             # as published by the Free Software Foundation, either version 3 of
10             # the License, or (at your option) any later version.
11             #
12             # This program is distributed in the hope that it will be useful,
13             # but WITHOUT ANY WARRANTY; without even the implied warranty of
14             # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15             # GNU General Public License for more details.
16             #
17             # You should have received a copy of the GNU General Public License
18             # along with this program. If not, see .
19             #
20             # -->8-->8-->8-->8--
21              
22             package Astro::XSPEC::TableModel;
23              
24 2     2   517648 use 5.00800;
  2         9  
  2         69  
25 2     2   11 use strict;
  2         5  
  2         56  
26 2     2   18 use warnings;
  2         8  
  2         60  
27              
28 2     2   10 use base 'Exporter';
  2         2  
  2         253  
29             our @EXPORT_OK = qw( write_table );
30              
31 2     2   9 use Carp;
  2         6  
  2         124  
32              
33 2     2   10 use List::Util qw( max );
  2         2  
  2         193  
34              
35 2     2   1681 use Params::Validate qw( :all );
  2         21140  
  2         506  
36              
37 2     2   934 use Astro::FITS::CFITSIO qw( :constants );
  0            
  0            
38             use Astro::FITS::CFITSIO::CheckStatus;
39             use Astro::FITS::Header::Item;
40             use Astro::FITS::Header::CFITSIO;
41              
42             our $VERSION = '0.01';
43              
44             # this rule is too broad
45             ## no critic ( ProhibitAccessOfPrivateData )
46              
47             sub write_table {
48              
49             my @fields;
50              
51             my %opt = validate_with(
52             params => \@_,
53             spec => {
54             output => { type => SCALAR },
55             model => 1,
56             units => 1,
57             additive => { type => BOOLEAN, default => '0' },
58             redshift => { type => BOOLEAN, default => '0' },
59             ipars => { type => ARRAYREF },
60             apars => { type => ARRAYREF, default => [] },
61             energy => { type => ARRAYREF },
62             keywords => { type => ARRAYREF, default => [] },
63             },
64             normalize_keys => sub { lc $_[0] },
65             );
66              
67             tie my $err, 'Astro::FITS::CFITSIO::CheckStatus';
68              
69             my $fits =
70             Astro::FITS::CFITSIO::create_file( "!$opt{output}",
71             $err = "Error creating $opt{output}:" );
72              
73             # -------------------------------------
74             # create PRIMARY HDU
75             $fits->create_img( 8, 0, 0, $err = 'Error creating primary HDU: ' );
76             _write_cards(
77             $fits,
78             $opt{keywords},
79             [
80             STRING =>
81             CONTENT => 'MODEL',
82             'spectrum file contains time intervals and event'
83             ],
84             [ STRING => MODLNAME => $opt{model}, 'Model name' ],
85             [ STRING => MODLUNIT => $opt{units}, 'Model units' ],
86             [
87             LOGICAL =>
88             REDSHIFT => $opt{redshift},
89             'If true then redshift will be included as a parameter'
90             ],
91             [
92             LOGICAL =>
93             ADDMODEL => $opt{additive},
94             'If true then this is an additive table model'
95             ],
96             [ STRING => HDUCLASS => 'OGIP', 'format conforms to OGIP standard' ],
97             [
98             STRING =>
99             HDUCLAS1 => 'XSPEC TABLE MODEL',
100             'model spectra for XSPEC'
101             ],
102             [ STRING => HDUVERS1 => '1.0.0', 'version of format' ] );
103              
104             # -------------------------------------
105             # create PARAMETERS HDU
106              
107             @fields = (
108             { type => 'NAME', form => '12A', unit => ' ' },
109             { type => 'METHOD', form => 'J', unit => ' ' },
110             { type => 'INITIAL', form => 'E', unit => ' ' },
111             { type => 'DELTA', form => 'E', unit => ' ' },
112             { type => 'MINIMUM', form => 'E', unit => ' ' },
113             { type => 'BOTTOM', form => 'E', unit => ' ' },
114             { type => 'TOP', form => 'E', unit => ' ' },
115             { type => 'MAXIMUM', form => 'E', unit => ' ' },
116             { type => 'NUMBVALS', form => 'J', unit => ' ' },
117             { type => 'VALUE', form => 'E', unit => ' ' },
118             );
119             my %field = ();
120             for my $col ( 1 .. @fields )
121             {
122             my $f = $fields[ $col - 1 ];
123             $f->{col} = $col;
124             $field{ $f->{type} } = $f;
125             }
126              
127             my @ipars = _validate_pars( \%field, $opt{ipars}, 'interpolated' );
128             my @apars = _validate_pars( \%field, $opt{apars}, 'additional' );
129              
130             # update form for VALUE to equal maximum NUMBVALS
131             my $nvals = max map { scalar @{ $_->{VALUE} } } @ipars;
132             $field{VALUE}{form} = $nvals . $field{VALUE}{form};
133              
134             _create_tbl( $fits, 'PARAMETERS', @fields );
135              
136             _write_cards(
137             $fits,
138             [],
139             [ STRING => HDUCLASS => 'OGIP', 'format conforms to OGIP standard' ],
140             [
141             STRING =>
142             HDUCLAS1 => 'XSPEC TABLE MODEL',
143             'model spectra for XSPEC'
144             ],
145             [
146             STRING =>
147             HDUCLAS2 => 'PARAMETERS',
148             'extension containing parameter info'
149             ],
150             [ STRING => HDUVERS1 => '1.0.0', 'version of format' ],
151             [
152             INT =>
153             NINTPARM => scalar @ipars,
154             'Number of interpolation parameters'
155             ],
156             [ INT => NADDPARM => scalar @apars, 'Number of additional parameters' ],
157             );
158              
159             # -------------------------------------
160             # write the parameters
161              
162             my $row = 0;
163             for my $par ( @ipars, @apars )
164             {
165             $row++;
166              
167             while ( my ( $key, $value ) = each %{$par} )
168             {
169             my $field = $field{$key};
170              
171             if ( $field->{type} eq 'VALUE' )
172             {
173             my $numbvals = @$value;
174             $fits->write_col_dbl( $field->{col}, $row, 1, $numbvals, $value,
175             $err );
176             $fits->write_col_dbl( $field{NUMBVALS}{col},
177             $row, 1, 1, $numbvals, $err );
178             }
179             else
180             {
181             my $datatype = $field->{form} =~ /A/ ? TSTRING : TDOUBLE;
182             $fits->write_col( $datatype, $field->{col}, $row, 1, 1, $value, $err );
183             }
184             }
185             }
186              
187             # -------------------------------------
188             # create the ENERGIES HDU
189              
190             @fields = (
191             { type => 'ENERG_LO', form => 'E', unit => ' ' },
192             { type => 'ENERG_HI', form => 'E', unit => ' ' },
193             );
194              
195             _create_tbl( $fits, 'ENERGIES', @fields );
196              
197             _write_cards(
198             $fits,
199             [],
200             [ STRING => HDUCLASS => 'OGIP', 'format conforms to OGIP standard' ],
201             [
202             STRING =>
203             HDUCLAS1 => 'XSPEC TABLE MODEL',
204             'model spectra for XSPEC'
205             ],
206             [
207             STRING =>
208             HDUCLAS2 => 'ENERGIES',
209             'extension containing energy bin info'
210             ],
211             [ STRING => HDUVERS1 => '1.0.0', 'version of format' ],
212             );
213              
214             # ick.
215             my @energy = @{ $opt{energy} };
216             my $nenergy = @energy - 1;
217             $fits->write_col_dbl( 1, 1, 1, $nenergy, \@energy, $err );
218             shift @energy;
219             $fits->write_col_dbl( 2, 1, 1, $nenergy, \@energy, $err );
220              
221             # -------------------------------------
222             # finally, create the SPECTRA HDU
223              
224             @fields = (
225             { type => 'PARAMVAL', form => @ipars . 'E', unit => ' ' },
226             { type => 'INTPSPEC', form => $nenergy . 'E', unit => $opt{units} },
227             );
228              
229             push @fields,
230             { type => "ADDSP$_", form => $nenergy . 'E', unit => $opt{units} }
231             for 1 .. @apars;
232              
233             _create_tbl( $fits, 'SPECTRA', @fields );
234              
235             _write_cards(
236             $fits,
237             [],
238             [ STRING => HDUCLASS => 'OGIP', 'format conforms to OGIP standard' ],
239             [
240             STRING =>
241             HDUCLAS1 => 'XSPEC TABLE MODEL',
242             'model spectra for XSPEC'
243             ],
244             [
245             STRING =>
246             HDUCLAS2 => 'ENERGIES',
247             'extension containing energy bin info'
248             ],
249             [ STRING => HDUVERS1 => '1.0.0', 'version of format' ],
250             );
251              
252             return $fits;
253             }
254              
255             # --------------------------------
256             sub _create_tbl {
257              
258             my ( $fits, $extname, @fields ) = @_;
259              
260             tie my $err, 'Astro::FITS::CFITSIO::CheckStatus';
261              
262             $fits->create_tbl( BINARY_TBL,
263             0,
264             scalar @fields,
265             [ map { $_->{type} } @fields ],
266             [ map { $_->{form} } @fields ],
267             [ map { $_->{unit} } @fields ],
268             $extname,
269             $err = "Error creating $extname HDU: "
270             );
271              
272             return;
273             }
274              
275             # --------------------------------
276             # write FITS cards to the CHU.
277             # Each card is an array with elements [ type, keyword, value, comment ]
278             # where type is a string recognized by Astro::FITS::Header.
279              
280             sub _write_cards {
281             my ( $fptr, $items, @cards ) = @_;
282              
283             my $hdr = Astro::FITS::Header::CFITSIO->new( fitsID => $fptr );
284              
285             $hdr->insert( -1, $_ ) foreach @$items;
286              
287             $hdr->insert(
288             -1,
289             Astro::FITS::Header::Item->new(
290             Type => $_->[0],
291             Keyword => $_->[1],
292             Value => $_->[2],
293             Comment => $_->[3] ) ) foreach @cards;
294             $hdr->writehdr( fitsID => $fptr );
295              
296             return;
297             }
298              
299             # --------------------------------
300             # validate parameter structures
301             sub _validate_pars {
302              
303             my ( $fields, $pars, $type ) = @_;
304              
305             # rewrite parameter hash with normalized keywords
306             my @pars;
307              
308             my %keys = map { $_ => 1 } keys %{$fields};
309              
310             # exclude NUMBVALS, as that's generated automatically
311             my @delkeys = qw( NUMBVALS );
312              
313             # "additional" parameters don't need this
314             push @delkeys, qw( METHOD VALUE )
315             if $type eq 'additional';
316              
317             delete @keys{@delkeys};
318             my @keys = keys %keys;
319              
320             my $row = 0;
321             for my $par (@$pars)
322             {
323             $row++;
324              
325             # normalize keywords
326             my %pars = map { uc $_ => $par->{$_} } keys %$par;
327              
328             my $name = $pars{NAME} || $row;
329              
330             # any missing?
331             my @missing = grep { !defined $pars{$_} } @keys;
332             croak( "$type parameter $name: missing attributes: ",
333             join( ', ', @missing ), "\n" )
334             if @missing;
335              
336             # any extra? really meant to catch errors in "additional" parameters
337             my @extra = grep { defined $pars{$_} } @delkeys;
338             croak( "$type parameter $name: illegal extra attributes: ",
339             join( ', ', @extra ), "\n" )
340             if @extra;
341              
342             push @pars, \%pars;
343             }
344              
345             return @pars;
346             }
347              
348             1;
349              
350             __END__