Molly - PDL module for reading and writing molly format spectra.
use lib "/home/sousun/trm/code/molly/perl"; use Molly;
open(MOLLY, "test.mol);
if(cmol(*MOLLY)){
$spec1 = rmol(*MOLLY);
}
if(cmol(*MOLLY)){
$spec2 = rmol(*MOLLY);
}
close(MOLLY);
$flux = $spec->getflux; $wave = $spec->gettwav; line $wave, $flux;
rmol e.g. $spec = rmol(*FILEHANDLE); -- reads a spectrum.
wmol e.g. $spec->wmol(*FILEHANDLE); -- writes a spectrum.
smol e.g. smol(*FILEHANDLE); -- skips a spectrum.
hmol e.g. $hdr = hmol(*FILEHANDLE); -- reads headers.
cmol e.g. cmol(*FILEHANDLE); -- checks existence of data before
applying rmol, smol or hmol.
gettwav e.g. $wave = $spec->gettwav; -- gets wavelengths.
gethwav e.g. $wave = $spec->gethwav; -- gets wavelengths, correcting to
heliocentric scale.
getflux e.g. $flux = getflux $spec; -- gets fluxes.
getferr e.g. $ferr = getferr $spec; -- gets flux uncertainties.
getcnts e.g. $cnts = getcnts $spec; -- gets counts.
getcerr e.g. $cerr = getcerr $spec; -- gets errors on counts.
achar e.g. $spec->achar('Object','M31') -- adds a character header item.
areal e.g. $spec->areal('Time','100') -- adds a real header item.
Reads a molly spectrum. Requires a file to be open and must be used after a call to 'cmol' e.g.
open(MOLLY, 'spectrum.mol')
cmol(*MOLLY);
$spec1 = rmol(*MOLLY);
cmol(*MOLLY);
$spec2 = rmol(*MOLLY);
cmol(*MOLLY);
$spec3 = rmol(*MOLLY);
close(MOLLY);
$spec1 etc are then npix by 3 dimesional piddles representing the
first 3 molly spectra. They contain the standard counts, errors and
flux arrays of molly. 'cmol' is needed to skip over first 4 bytes
(and can be tested for success if number of spectra not known). This
is to work around a problem with testing pdls because it was otherwise
difficult to test for the end of the file.
NB: At the end of each read, one is left at the start of the next spectrum.
Much more information is stored in headers as a hash which can be
obtained via:
$hdr = $spec->gethdr;
Then
%{$hdr} are the headers (e.g. $$hdr{'Object'} is the object name)
$$hdr{'First record'} is a reference to an array containing the
first record i.e. fcode, units, npix, narc, nchar, ndoub,
nintr, nreal
$$hdr{'Arc'} is a reference to the arc coeffcients array
$$hdr{'Character names'} is a reference to an array containing the
names of character header items.
$$hdr{'Double names'} is a reference to an array containing the
names of double header items.
$$hdr{'Integer names'} is a reference to an array containing the
names of integer header items.
$$hdr{'Real names'} is a reference to an array containing the
names of real header items.
The latter 4 are required in order to be able to write molly spectra.
Writes a molly spectrum to a file opened for output on filehandle.
e.g.
open(OUTPUT,">out.mol");
wmol($spec,*OUTPUT);
close(OUTPUT);
or
open(OUTPUT,">out.mol");
$spec->wmol(*OUTPUT);
close(OUTPUT);
Skips over a molly spectrum more efficiently than using rmol to read it.
Like rmol you need first to use cmol to skip first 4 bytes if possible.
smol reads first record to work out how many bytes to skip and then skips
them.
Reads molly headers, then skips to end of spectrum. Like 'rmol'
and 'smol', must be used after a call to 'cmol'. Slightly more
efficient than using 'rmol' to read everything.
e.g.
$hdr = hmol(*MOLLY);
print "Object name = $$hdr{'Object'}\n";
=cut
sub hmol {PDL->hmol(@_)}
sub PDL::hmol{ barf 'Usage: $a =
hmol(FILEHANDLE)' if $#_ != 1; my ($class, $filehandle) = @_;
# Setup header hash
my $header = {};
my ($fcode,$units,$npix,$narc,$nchar,$ndoub,$nintr,$nreal);
my ($buf,$template,$nhead,$nbytes,$ncoeff,@header_names,@header_values);
# Read and unpack the first record
read($filehandle, $buf, 48) == 48
or barf "Couldn't read first record";
$$header{'First record'} = [unpack("ia16i6x4", $buf)];
$fcode = $$header{'First record'}[0];
$fcode == 3 or
barf "fcode = $fcode. Can only read fcode=3 molly data\n";
$npix = $$header{'First record'}[2];
$narc = $$header{'First record'}[3];
$nchar = $$header{'First record'}[4];
$ndoub = $$header{'First record'}[5];
$nintr = $$header{'First record'}[6];
$nreal = $$header{'First record'}[7];
# Read second record of header item names
$nhead = $nchar + $ndoub + $nintr + $nreal;
$nbytes = 16*$nhead+8;
read($filehandle, $buf, $nbytes) == $nbytes
or barf "Couldn't read second record";
$template = "x4";
for (1..$nhead) {
$template .= "A16";
}
@header_names = unpack($template, $buf);
# Remove trailing blanks.
foreach(@header_names){
s/ *$//;
}
# Read and interpret header item values
$nbytes = 32*$nchar + 8*$ndoub + 4*$nintr + 4*$nreal + 8;
read($filehandle, $buf, $nbytes) == $nbytes
or barf "Couldn't read third record";
$template = "x4";
for (1..$nchar) {
$template .= "A32";
}
$template .= "d$ndoub i$nintr f$nreal";
@header_values = unpack($template, $buf);
# Load the hash
for (0..$#header_names) {
$$header{$header_names[$_]} = $header_values[$_];
}
# Read arc coefficients, store as an array reference in the hash.
$nbytes = $narc < 0 ? -8*$narc + 8: 8*$narc + 8;
$ncoeff = $narc < 0 ? -$narc : $narc;
read($filehandle, $buf, $nbytes) == $nbytes
or barf "Couldn't read arc coefficients record";
$template = "x4d$ncoeff";
$$header{'Arc'} = [unpack($template, $buf)];
# Now read counts, errors, fluxes into the pdl
$nbytes = 12*$npix+8;
read($filehandle, $buf, $nbytes) == $nbytes
or barf "Failed to skip data";
return $header; }
Skips over first 4 bytes of a molly spectrum (which contains ignorable data) and returns 0 or 1 depending on success.
cmol should be used before every call to 'rmol', 'smol' and 'hmol' in order to be get the correct position correctly in the file. It is also needed when a file has an unknown number of spectra:
open(INPUT, 'file')
while(cmol(*INPUT)){
$spec = rmol(*INPUT);
.
.
Process spectra
.
.
}
close(INPUT);
returns a pdl of wavelengths generated from the arc coefficients. No correction to a heliocentric scale is applied (marked by the second 't') although of course the arc coefficients may already represent a heliocentric scale.
e.g.
$wave = $spec->gettwav; $wave = gettwav $spec;
e.g. $flux = $spec->getflux; $flux = getflux $spec;
returns a pdl of flux uncertainties.
e.g. $ferr = $spec->getferr; $ferr = getferr $spec;
returns a pdl of counts. This is a simple case of returning the first row of the pdl.
e.g. $cnts = $spec->getcnts; =cut
*getcnts = \&PDL::getcnts;
sub PDL::getcnts{ barf 'Usage: $flux =
getcnts($spec)' if $#_ > 1; my $spec = shift;
my($cnts); $cnts =
$spec->slice(':,(0)')->copy;
return $cnts; }
returns a pdl of errors on counts. This is a simple case of returning the second row of the pdl.
e.g. $cerr = $spec->getcerr;
The item has name $name and value $item.
adds a real header item to a molly piddle.
The item has name $name and value $item.