Using Perl with molly

Perl can allow you to do things with molly that are otherwise rather difficult. One can send commands from Perl directly to molly, interact with the output files and then send revised commands to molly. This is especially useful for do-loop type operations. Perl is fairly easy to use so do not be put off.

Finally for plotting with pgperl its useful to be able to read in molly files directly and subroutines to do this are listed in this file.

A fairly common problem is that you want to run molly, and then take an action dependent upon the results. There is apparently a way of doing this in one go but I believe that it requires rewriting molly code which I would prefer not to do. Instead one can apply the following procedure: (1) run molly once, (2) get out of it process results with perl, (3) run molly again with actions based upon results of earlier run. You can see this in some of the scripts below.

Here are some example scripts:

  • reorder -- reorder spectra by JD.
  • strip -- strip trailing 's from object names.
  • slit -- slit correct spectra.
  • molly -- Perl subroutines to read molly data.
  • ascii -- read a set of ascii files into molly.

    • reorder

    Here is an example in which Perl is used to reorder a molly file, in this case so that the spectra are in order of time. If the file was called 'reorder' it would be invoked as 'reorder input output'. 
    
#!/usr/local/bin/perl

# Perl routine to reorder a molly file to make the spectra ascend
# in JD. This works for a maximum of $NMAX spectra and requires
# molly to be able to fit in 2*$NMAX spectra. Change the line below
# if you need to reorder more.
#
# Method: read in spectra, dump out list of RJDs, leave molly. Sort RJDS, 
# read in spectra a second time, reorder spectra dump to new file.
#
# Note that a file called junk is written.

$NMAX = 100;

# Arguments:  

# start up molly, read in spectra dump list of RJDs to a file
# close down.

open(MOLLY, "|molly") || die "Failed to start molly\n";
print MOLLY "load $ARGV[0] 1 $NMAX 1\n";
print MOLLY "flis 1 $NMAX junk 1 RJD\n";
print MOLLY "q\ny\n";
close(MOLLY);

# read RJDs in, make an array which associates the spectrum 
# number equivalent to each RJD

open(RJD, "junk") || die "Failed to open junk\n";
while(){
    ($number,$rjd) = split(' ');
    $which{$rjd} = $number;
}
close(RJD);

# sort the RJDs into an ascending order array

@order = sort byrjd keys(%which);

# now start molly again, load spectra reorder them
# and dump out

open(MOLLY, "|molly") || die "Failed to start molly\n";
print MOLLY "load $ARGV[0] 1 $NMAX 1\n";

$i = 0;
foreach(@order){
    $from = $which{$_};
    $to   = $NMAX + (++$i);
    print MOLLY "copy $from $from $to 1\n";
}
$start = $NMAX + 1;

print MOLLY "write $ARGV[1] $start $to N\n";
print MOLLY "q\ny\n";
close(MOLLY);

# comparison subroutine needed by sort command

sub byrjd {$a <=> $b;}

    strip

    This is similar to reorder but the purpose is to clean up object names (in this case it just strips off trailing '). Perl has powerful regular expression utilities which could be used on tougher problems. 
    #!/usr/local/bin/perl

# Perl routine to strip ' off end of object names
# Arguments:  

$NMAX = 100;

# start up molly, read in spectra dump list of object names.

open(MOLLY, "|molly") || die "Failed to start molly\n";
print MOLLY "load $ARGV[0] 1 $NMAX 1\n";
print MOLLY "flis 1 $NMAX junk 1 Object\n";
print MOLLY "q\ny\n";
close(MOLLY);

# Read in object names, removing final ' if present

open(OBJECT_NAMES, "junk") || die "Failed to open junk\n";
$i = 0;
while(){
    if(/ *\d+ (.*)/){
	($object[$i++] = $1) =~ s/\'//;
    }
}
close(OBJECT_NAMES);

# Generate a file to read into the header editor

open(OBJECT_NAMES, ">junk") || die "Failed to open junk\n";
print OBJECT_NAMES "Object\nC\n";
for(@object){
    $name = substr($_,0,32);
    print OBJECT_NAMES "\"$name\"\n";
}
close(OBJECT_NAMES);

# finally read in again, modify names and dump

open(MOLLY, "|molly") || die "Failed to start molly\n";
print MOLLY "load $ARGV[0] 1 $NMAX 1\n";
print MOLLY "edit 1 $NMAX\n";
print MOLLY "file\njunk\nquit\n";
print MOLLY "write $ARGV[1] 1 $NMAX n\n";
print MOLLY "q\ny\n";
close(MOLLY);

    slit

    The following program slit corrects spectra by matching objects within files of narrow and wide slit spectra. It matches by name only but assumes no ordering of spectra in either input file. It stores a header item indicating which wide slit spectrum has been used (if any). A mask must have been set-up prior to using this program. Checks are also made for matching positions. See comments for more details. 
    #!/usr/local/bin/perl

#
# Perl routine to apply slit corrections with molly
#
# Arguments: 
#
#  narrow -- input file with narrow slit spectra
#  wide   -- input file with wide slit spectra. Should only have one/object
#            or else only the last one will be used
#  output -- output file
#
# There must be a file called "slit.msk" to mask the spectra
# prior to polynomial fitting.
#
# $NMAX  = maximum number of spectra in either input spectra
# $NPOLY = number of poly coefficients
#
# molly must be able to handle 4*NARROW+WIDE spectra for this to work
# but it is not checked.
#
# Operation:
#
#  1) molly loads input files, dumps object names to two files. The record
#     numbers of the wide slit spectra are also dumped.
#  2) Perl loads object names, stripping trailing '. It also loads the
#     record numbers of the wide slit spectra.
#  3) For each wide slit spectrum, narrow slit spectra with precisely
#     matching names are selected and copied into contiguous slots.
#     A check on the relative positions is also made. A warning is printed
#     if they are too discrepant and no slit correction is performed.
#  4) The wide slit spectrum is then copied a matching number of times
#  5) The wide slit spectrum is rebinned onto the same scale as the narrow
#     slit spectra
#  6) The results from 5 are divided into the narrow slit spectra
#  7) Poly fits are made to the output from 6
#  8) The poly fits are divided into the narrow slit spectra
#  9) The now corrected spectra are copied back to their original places.
#     The record number of the wide slit star used is inserted in the headers
#     as an integer item called "Wide slit star"
#  10) The spectra ae dumped to disk
#

# ensure not overwriting input files

(($ARGV[0] ne $ARGV[2]) && ($ARGV[1] ne $ARGV[2])) || 
    die "Will not overwrite inputs\n";

$NMAX  = 100;
$NPOLY = 1;                     # number of poly coefficients
$DMAX  = 1.;                    # maximum separation in arcminutes

print STDERR "Slit will perform slit corrections\n";
print STDERR "    Narrow slit file (input): $ARGV[0]\n";
print STDERR "      Wide slit file (input): $ARGV[1]\n";
print STDERR "Slit corrected file (output): $ARGV[2]\n";
print STDERR "  Maximum number of spectra set at: $NMAX\n";
print STDERR "Number of poly coefficients set at: $NPOLY\n";
print STDERR "Maximum separation (arcmin) set at: $DMAX\n\n";

# First run molly to list object names (and record numbers of
# wide slit spectra.

open(STDOUT, ">logfile") || die "Failed to open log file\n";

open(MOLLY, "|molly")   || die "Failed to start molly\n";
print MOLLY "confirm\n";
print MOLLY "load $ARGV[0] 1 $NMAX 1\n";
print MOLLY "flis 1 $NMAX junk1 3 Object RA Dec\n";
print MOLLY "clear 1 $NMAX sure\n";
print MOLLY "load $ARGV[1] 1 $NMAX 1\n";
print MOLLY "flis 1 $NMAX junk2 4 Object Record RA Dec\n";
print MOLLY "q\ny\n";
close(MOLLY);

# Read in narrow slit object names stripping trailing ' and
# storing RAs and Decs

open(INPUT, "junk1") || die "No junk1 around\n";
$i = 0;			       
while(){		     
    chop;
    if(/^ *\d* *(.*) *\' *(\S*) *(\S*) *$/){
        $narrow_ra[$i]  = $2;
        $narrow_dec[$i] = $3;
	($narrow[$i++] = $1) =~ s/ *$//;
    }
}
close(INPUT);
$num_narrow = $i;

# Read in wide slit object names stripping trailing ' and storing
# record numbers, RA and Dec

open(INPUT, "junk2") || die "No junk2 around\n";
$i = 0;			       
while(){
    chop;
    if(/^ *\d* *(.*) *\' *(\S*) *(\S*) *(\S*) *$/){
	$wide_rec[$i] = $2;
        $wide_ra[$i]  = $3;
        $wide_dec[$i] = $4;
	($wide[$i++]  = $1) =~ s/ *$//;
    }
}
close(INPUT);

$num_wide = $i;			

# Now apply slit corrections.

open(MOLLY, "|molly")   || die "Failed to start molly\n";
print MOLLY "confirm\n";

# read input files

$nstart = 1;
$nend   = $num_narrow;
print MOLLY "load $ARGV[0] $nstart $nend 1\n";
$wstart = $num_narrow + 1;
$wend   = $wstart + $num_wide - 1;
print MOLLY "load $ARGV[1] $wstart $wend 1\n";
$bstart = $num_narrow+$num_wide+1;

$first = 1;
$twopi = 8.*atan2(1.,1.);

for($nw = 0; $nw < $num_wide; $nw++){

# First check for matching narrow slit spectra

    $spec   = $wstart + $nw;
    $wname  = $wide[$nw];
    $wra    = $wide_ra[$nw];
    $wdec   = $wide_dec[$nw];
    $cosd   = cos($twopi*$wdec/360.);

    $record = $wide_rec[$nw];
    @list = ();
    $nlist = 0;
    for($i = 0; $i < $num_narrow; $i++){
	if($narrow[$i] eq $wname){
            $diff = sqrt((15.*60.*$cosd*($narrow_ra[$i]-$wra))**2 +
                         (60.*($narrow_dec[$i]-$wdec))**2);
            if($diff < $DMAX){
                $moveto = $bstart + $nlist;
                $list[$nlist++] = ($nmatch = $nstart+$i);
                print MOLLY "copy $nmatch $nmatch $moveto 1\n";
            }
            else{
                print STDERR "Although the names of wide slit spectrum $nw"
                    . " and narrow slit spectrum $i match, they differ\n"
                        ." in poisiotn by $diff arcminutes\n"
                            . "and no correction will be applied\n";
            }
	}
    }

# There are narrow slit spectra to be fixed

    if($nlist > 0){
	$n1 = $bstart;		# range of narrow slit
	$n2 = $n1 + $nlist - 1;	# spectra
	$n3 = $n2 + 1;		# range of multiply copied wide slit
	$n4 = $n3 + $nlist - 1;	# spectrum 
	$n5 = $n4 + 1;		# range of rebinned wide slit
	$n6 = $n5 + $nlist - 1;	# spectra

	print MOLLY "copy $spec $spec $n3 $nlist\n";
	print MOLLY "tbin $n3 $n4 $n5 $n1 $n2 Q\n";
	print MOLLY "div $n1 $n2 $n5 $n6 $n3\n";

# $n3 .. $n4 now have ratio spectra

	if($first){
	    print MOLLY "pfit $n3 $n4 $n5 $NPOLY -3 3 10 0 y\n";
	    print MOLLY "load\nslit.msk\nq\n";
	    $first = 0;
	}
	else{
	    print MOLLY "pfit $n3 $n4 $n5 $NPOLY -3 3 10 0 n\n";
	}

# $n5 .. $n6 now have poly fits

	print MOLLY "div $n1 $n2 $n5 $n6 $n1\n";

# $n1 .. $n2 have slit corrected spectra. Copy back
# edit headers to include record number of wide slit spectrum.

	for($j = 0; $j < $nlist; $j++){
	    $corr = $n1 + $j;
	    print MOLLY "edit $corr $corr\n";
	    print MOLLY "term\nWide slit star\nI\n$record\nquit\n";
	    print MOLLY "copy $corr $corr $list[$j] 1\n";
	}
    }
}				# next wide slit spectrum

print MOLLY "write $ARGV[2] 1 $num_narrow n\n";
print MOLLY "q\ny\n";
close(MOLLY);
close(STDOUT);

    Perl reads molly

    For those who use pgperl the obvious way to plot a molly spectrum is to dump it to an ASCII file and read it in to Perl in the usual way. However this has drawbacks that you need a separate file for each spectrum and if you change the molly spectra you have to remember to change the ASCII spectra too. An alternative is to read the data into perl directly which can be done if you know the byte order. I have done this for molly files dumped on my SPARCstation and the same program will probably work on most UNIX machines.

    There is one drawback to this. It is SLOW. For more than about 20 spectra it starts to become impractical. Karl Glazebrook (the author of pgperl) has now developed a way of handling binary arrays inside perl. I have not used it yet but you may want to if you want to follow this route.

    The subroutines are listed below, but first here is an example showing them being called (assumes Perl5): 

    
$file = 'spectra.mol';
open(DATA, $file) || die "Failed to open $file\n";
$i = 0;

# Read in data and store the fluxes in a series of arrays. The counts
# and errors are not required later and are overwritten each time
# the loop breaks when the end of the file is reached

while(&read_molly(DATA, \$npix, \$narc, \%header, \@arc, \@counts, \@errors, 
		  \@{$flux[$i]})){
 
# calculate the wavelengths from the arc coefficients

    @{$wave[$i]} = &calculate_wavelength(\@arc, $narc, $npix);

# remove the Earth velocity to get a heliocentric scale

    &remove_earth_velocity(\@{$wave[$i]}, \%header);

# Correct the fluxes for zero count pixels

    &correct_fluxes(\@counts,\@errors,\@{$flux[$i]},\@{$ferr[$i]});
    $i++;
}
close(DATA);
$nspec = $i;
print "Read $nspec spectra\n";

    The flux arrays can be used later. For example assuming the plot has been opened etc, here is some pgperl which plots the spectra with an offset from one to the next: 

    
$off = 10.;
$add = 0.;
for($n = 0; $n < $nspec; $n++){
    foreach(@{$flux[$n]}){
	$_ += $add;
    }
    $add += $off;

# Plot spectrum

    &pgbin(scalar(@{$flux[$n]}),\@{$wave[$n]},\@{$flux[$n]},1);
}

    Here are the subroutines used above: 

    #read_molly
#
# read_molly($filehandle, \$npix, \$narc, \%header, \@arc, \@counts, \@errors, 
#            \@flux)
#
# Subroutine for reading a molly spectrum, assuming currently positioned
# at start of headers. Returns associative header array, counts, errors
# and fluxes. Header parameters can be accessed as in 
#
# $object_name = $header{'Object'};
#
#read_molly

sub read_molly{

    local($filehandle, *npix, *narc, *header, *arc,
	  *counts, *errors, *flux) = @_;
    local($fcode,$units,$nchar,$ndoub,$nintr,$nreal);
    local($buffer, @data);
    local($nhead,$template,$nbytes,$ndata,$ncoeff);
    local(@header_values,@header_names);
    
    @arc           = ();
    @data          = ();
    @counts        = ();
    @errors        = ();
    @flux          = ();
    @header_values = ();
    @header_names  = ();
    %header        = ();
#
# Read and interpret first line of header
#

    if(read($filehandle, $buffer, 52)) {
	
	($fcode,$units,$npix,$narc,$nchar,$ndoub,$nintr,$nreal) =
	    unpack("x4ia16i7", $buffer);

#
# Read and interpret header item names
#

	$nhead = $nchar + $ndoub + $nintr + $nreal;
	read($filehandle, $buffer, 16*$nhead+8);
	$template = "x4";
	for (1..$nhead) {
	    $template .= "A16";
	}
	@header_names = unpack($template, $buffer);
	foreach(@header_names){
	    s/ *$//;
	}

#
# Read and interpret header item values
#

	$nbytes = 32*$nchar + 8*$ndoub + 4*$nintr + 4*$nreal + 8;
	read($filehandle, $buffer, $nbytes);
	$template = "x4";
	for (1..$nchar) {
	    $template .= "A32";
	}
	$template .= "d$ndoub i$nintr f$nreal";
	@header_values = unpack($template, $buffer);
 
#
# Convert header to an associative array
#
	
	for (0..$#header_names) {
	    $header{$header_names[$_]} = $header_values[$_];
	}

#
# Read arc coefficients
#

	$nbytes = $narc < 0 ? -8*$narc + 8 : 8*$narc + 8;
	$ncoeff = $narc < 0 ? -$narc : $narc;
	
	read($filehandle, $buffer, $nbytes);
	$template = "x4d$ncoeff";
	@arc = unpack($template, $buffer);

#
# Read data
#

	$nbytes = 12*$npix + 8;
	$ndata  =  4*$npix;
	read($filehandle, $buffer, $nbytes);
	$template = "x4f$ndata";
	@data = unpack($template, $buffer);
	@counts = @data[0..$npix-1];
	@errors = @data[$npix..2*$npix-1];
	@flux   = @data[2*$npix..3*$npix-1];
    }
}				

#skip_molly
#
# skip_molly($filehandle, $nskip)
#
# Subroutine for skipping a molly spectrum, assuming currently positioned
# at start of headers. 
#
#skip_molly

sub skip_molly{			

    local($filehandle, $nskip) = @_;
    local($fcode,$units,$npix,$narc,$nchar,$ndoub,$nintr,$nreal);
    local($nbytes, $buffer);

    for (1..$nskip) {
#
# Read and interpret first line of header
#

	if(read($filehandle, $buffer, 52)) {
	    
	    ($fcode,$units,$npix,$narc,$nchar,$ndoub,$nintr,
	     $nreal) = unpack("x4ia16i7", $buffer);
			
#
# Compute number of bytes to skip
#
	
	    $nbytes  = 16*($nchar + $ndoub + $nintr + $nreal)+8;
	    $nbytes += 32*$nchar + 8*$ndoub + 4*$nintr + 4*$nreal + 8;
	    $nbytes += $narc < 0 ? -8*$narc + 8 : 8*$narc + 8;
	    $nbytes += 12*$npix + 8;

	    read($filehandle, $buffer, $nbytes);
	}
    }
}

#calculate_wavelength
#
# @wave = calculate_wavelength(\@arc, $narc, $npix)
#
# Calculates wavelength scale equivalent to molly arc coefficients. No
# correction for the Earth's velocity is applied.
#
#calculate_wavelength

sub calculate_wavelength{

    local(*arc,$narc,$npix) = @_;
    local($i,$j,$fac,$z);
    local(@wave) = ();
	
    for $i (0..$npix-1) {
	$wave[$i] = 0.;
	$z = ($i+1)/$npix;
	if($narc > 0) {
	    $fac = 1.;
	    for $j (0..$narc-1) {
		$wave[$i] += $arc[$j]*$fac;
		$fac  *= $z;
	    }
	}
	elsif($narc < 0) {		
	    $fac = 1.;
	    for $j (0..-$narc-1) {
		$wave[$i] += $arc[$j]*$fac;
		$fac  *= $z;
	    }
	    $wave[$i] = exp($wave[$i]);
	}
	else {
	    $wave[$i] = $i + 1;
	}
	
    }
    @wave;
}				

#remove_earth_velocity
#
# remove_earth_velocity(\@wave, \%header)
#
# Looks up Earth velocity and corrects wavelengths for it.
#
#remove_earth_velocity

sub remove_earth_velocity{

    local(*wave, *header) = @_;
    local($i,$j,$fac,$z);

    $corr_factor = 1. - $header{'Vearth'}/2.997925e5;
    
    for (@wave) {
	$_ *= $corr_factor;
    }
}

#correct_fluxes
#
# correct_fluxes(\@counts, \@errors, \@flux, \@ferr)
#
# Checks for 0 count pixels and corrects fluxes. Also generates
# array of flux errors.
#
#correct_fluxes

sub correct_fluxes{

    local(*counts, *errors, *flux, *ferr) = @_;
    local($i);
    
    @ferr = ();
    for $i (0..$#flux) {
	if($counts[$i] == 0.) {
	    $ferr[$i] = $errors[$i]/$flux[$i];
	    $flux[$i] = 0.;
	}
	else {
	    $ferr[$i] = $flux[$i]*$errors[$i]/$counts[$i];
	}
	$ferr[$i] = $ferr[$i]*$errors[$i] > 0 ? $ferr[$i] : -$ferr[$i];
    }
}

#wave_to_kms
#
# @velocity = wave_to_kms(\@wave, $wzero)
#
# Converts from wavelength to velocity in km/s given a central
# wavelength corresponding to zero velocity which must have the
# same units as the wavelength array.
#
#wave_to_kms

sub wave_to_kms{

    local(*wave, $wzero) = @_;
    local($i,$vlight);
    local(@velocity) = ();
    
    $vlight = 2.997925e5;

    for $i (0..$#wave) {
	$velocity[$i] = $vlight*($wave[$i]-$wzero)/$wzero;
    }
    @velocity;
}


1;

    ascii

    You may be trying to use molly on data from another system. At the moment you are stuck with FIGARO or ASCII for converting data formats. The ASCII reader 'lasc' only handles them one at a time. It could be made to run on lists but it does not. This is a classic case for using Perl to run molly. Here is a program written by Emilios Harlaftis to read a set of ASCII files and dump them into a molly file. The spectra are assumed to be in Angstroms and mJy. 
    
#!/usr/local/bin/perl
#
# Perl routine to read in ascii files into molly and dump out
# molly data file. It assumes that the list of ascii files is held
# in a file called 'files.lis'
# 
# Arguments: 
#
 
open(STDOUT, ">logfile") || die "Failed to open log file\n";
open(MOLLY, "| molly")   || die "Failed to start molly\n"; 
 
# load ascii spectra

open(FILE, "files.lis") || die "failed to open files.lis\n";

$nspec = 0;
while(){
   chop;
   $nspec++;
   print MOLLY "lasc $_ $nspec 1 2 -1 A MJY 0.01\n";
}
close(FILE);

# write out file, close down.

print MOLLY "write $ARGV[0] 1 $nspec N\n";
 
print MOLLY "q\ny\n";
close(MOLLY);
    Tom Marsh, Southampton, trm@astro.soton.ac.uk.