User:Christopher Thomas/spectrum script v1

I wrote this Perl script to generate pictures of emission spectra, per this archived thread from WT:PHYS. This script may be freely copied and modified. Sample images from this version of the script are here.

Note: There are hooks for telling the script to fetch spectrum data from NIST, but this feature wasn't implemented. As-is, you have to fetch the table manually (in text mode, if I remember correctly), and store it as a local file for the script to read. The rest of the help screen should be accurate.

NIST's database query form is at http://physics.nist.gov/PhysRefData/ASD/lines_form.html --Christopher Thomas (talk) 17:36, 16 December 2009 (UTC)


 * 1) !/usr/bin/perl
 * 2) Makes an emission spectrum plot from NIST data.
 * 3) Original version written by Christopher Thomas (2009).
 * 4) Usage: make_spectrum  [options]
 * 5) is a file containing a NIST data table, in text format.
 * 6) is the name of an element to search NIST for.
 * 7) is the name of the SVG image file to write to.
 * 8) [options] are optional flags modifying behavior.
 * 9) Options are:
 * 10) --fetchnist   Attempts to fetch NIST data for.
 * 11)               Default is to read from instead.
 * 1) Options are:
 * 2) --fetchnist   Attempts to fetch NIST data for.
 * 3)               Default is to read from instead.
 * 1)               Default is to read from instead.


 * 1) Includes
 * 1) Includes

use strict;


 * 1) Constants
 * 1) Constants

my ($lambda_min_min, $lambda_max_max); $lambda_min_min = 200; $lambda_max_max = 2000;
 * 1) Ultimate limits to wavelength.

my ($vis_min, $vis_max); $vis_min = 380; $vis_max = 750;
 * 1) Visible spectrum limits.

my (%rgbmap_red, %rgbmap_green, %rgbmap_blue); %rgbmap_red = ( $lambda_max_max + 1 => 100,  750 => 100,  725 => 200,  675 => 250,  600 => 200,  525 => 0,  450 => 0,  425 => 100,  400 => 125,  375 => 100,  $lambda_min_min - 1 => 100 ); %rgbmap_green = ( $lambda_max_max + 1 => 50,  750 => 50,  725 => 0,  675 => 0,  600 => 200,  525 => 250,  500 => 200,  450 => 0,  400 => 0,  375 => 50,  $lambda_min_min - 1 => 50 ); %rgbmap_blue = ( $lambda_max_max + 1 => 0,  525 => 0,  500 => 200,  450 => 250,  400 => 200,  375 => 75,  $lambda_min_min - 1 => 75 );
 * 1) Piecewise-linear spectrum definition.
 * 2) Roughly matches FIXME:NAME 's mapping function.
 * 3) Tweaked to remove banding, and to put R/G/B at 675/525/450 nm.
 * 4) Also tweaked IR and UV values to be obviously-different.

my (@specbar_lambdas); @specbar_lambdas = ( $lambda_max_max,  750, 725, 675, 600, 525, 500, 450, 425, 400, 375,  $lambda_min_min );
 * 1) Spectrum bar landmark wavelengths.


 * 1) Configuration Variables
 * 1) Configuration Variables

my ($spacing_big_frac, $spacing_small_frac); $spacing_big_frac = 0.1; $spacing_small_frac = 0.03;
 * 1) Spacing. Make this a fraction of the wavelength range.

my ($lambda_min, $lambda_max); $lambda_min = 200; $lambda_max = 1000;
 * 1) Wavelength range to query data for.
 * 1) Test a smaller range.
 * 2) $lambda_min = 500;
 * 3) $lambda_max = 600;

my ($inten_bogus, $inten_default); $inten_bogus = 0; $inten_default = 1;
 * 1) Bogus and default line intensity values.

my ($img_width, $img_height); my ($img_spacing_big, $img_spacing_small); $img_width = 200; $img_height = 200; $img_spacing_big = 10; $img_spacing_small = 3;
 * 1) Actual spacing/size values. These will get reinitialized!

my ($inten_log_decades); $inten_log_decades = 2.5;
 * 1) Number of decades to map for log-scale colours.


 * 1) Functions
 * 1) Functions


 * 1) Performs a piecewise-linear mapping.
 * 2) Finds the nearest mapping vertices on either side of the input value,
 * 3) and interpolates the output value from these.
 * 4) Arg 1 is the value to be mapped.
 * 5) Arg 2 points to a hash containing mapping values.
 * 6) Returns the mapped value.

sub MapPWL { my ($inval, $map_p); my ($result); my ($lessidx, $moreidx, $lessval, $moreval); my ($delta_x, $delta_y); my ($thisidx);

# Initialize. $result = 0;

# Get args. $inval = $_[0]; $map_p = $_[1];

# Proceed if args are valid. if ((defined $inval) && (defined $map_p)) {   # Find the left and right bounding vertices.

$lessidx = undef; $moreidx = undef;

foreach $thisidx (sort keys %$map_p) {     if ($thisidx <= $inval) {       # Find the biggest index less than the input. if (!(defined $lessidx)) { $lessidx = $thisidx; } elsif ($thisidx > $lessidx) { $lessidx = $thisidx; } }     else {       # Find the smallest index greater than the input. if (!(defined $moreidx)) { $moreidx = $thisidx; } elsif ($thisidx < $moreidx) { $moreidx = $thisidx; } }   }

# Compute the interpolated value.

if ((defined $lessidx) && (defined $moreidx)) {     $lessval = $$map_p{$lessidx}; $moreval = $$map_p{$moreidx};

$delta_x = $moreidx - $lessidx; $delta_y = $moreval - $lessval;

# Sanity, even though PWL map should guarantee this. if ($delta_x > 1.0e-20) {       $result = $lessval; $result += $delta_y * ($inval - $lessidx) / $delta_x; }   }  }

# Return the mapped value. return $result; }


 * 1) Translates a wavelength (in nm) into an RGB value (range 0-255).
 * 2) Arg 1 is the wavelength.
 * 3) Returns a hash (by value), containing "red", "green", and "blue"
 * 4) entries.

sub WaveToRGB { my ($lambda); my (%rgb);

# Initialize. %rgb = ( 'red' => 0, 'green' => 0, 'blue' => 0);

# Perform mapping. if (defined ($lambda = $_[0])) {   $rgb{red} = int(MapPWL($lambda, \%rgbmap_red)); $rgb{green} = int(MapPWL($lambda, \%rgbmap_green)); $rgb{blue} = int(MapPWL($lambda, \%rgbmap_blue)); }

# Return the resulting triplet. return %rgb; }


 * 1) Performs a diagnostics sweep of the spectrum colours.
 * 2) Dumps output to STDERR.
 * 3) No arguments.
 * 4) No return value.

sub DebugRGBVals { my ($lambda); my (%rgb);

for ($lambda = $lambda_min_min;   $lambda <= $lambda_max_max;    $lambda += 10) {   %rgb = WaveToRGB($lambda); print STDERR "For $lambda nm, got (" . $rgb{red} . ',' . $rgb{green}     . ',' . $rgb{blue} . ").\n"; } }


 * 1) Retrieves text-format spectrum data from NIST.
 * 2) Arg 1 is the name of the element to fetch.
 * 3) Arg 2 points to an array to store raw data in.
 * 4) Returns 1 if successful and 0 if not.

sub GetDataNIST { my ($ename, $data_p); my ($is_ok);

# Initialize. $is_ok = 0;

# Get args. $ename = $_[0]; $data_p = $_[1];

# Proceed if we have args. if ((defined $ename) && (defined $data_p)) {   # Initialize. @$data_p = ;

# FIXME - NYI. print STDERR "### NIST fetch not yet implemented!\n"; }

# Return the resulting error flag. return $is_ok; }


 * 1) Retrieves text-format spectrum data from a file.
 * 2) Arg 1 is the name of the file to read.
 * 3) Arg 2 points to an array to store raw data in.
 * 4) Returns 1 if successful and 0 if not.

sub GetDataFile { my ($iname, $data_p); my ($is_ok);

# Initialize. $is_ok = 0;

# Get args. $iname = $_[0]; $data_p = $_[1];

# Proceed if we have args. if ((defined $iname) && (defined $data_p)) {   # Initialize. @$data_p = ;

# Try to open the data file. if (!open(INFILE, "<$iname")) {     print STDERR "### Unable to read from \"$iname\".\n"; }   else {     # Read file contents as-is. @$data_p = ;

# Close the input file. close(INFILE);

# Report success. $is_ok = 1; } }

# Return the resulting error flag. return $is_ok; }


 * 1) Extracts spectrum line data from raw NIST data tables.
 * 2) Arg 1 points to an array containing the raw NIST data table as text.
 * 3) Arg 2 points to a hash to store line and intensity data in.
 * 4) Returns 1 if successful, and 0 if not.

sub ExtractLines { my ($data_p, $lines_p); my ($is_ok); my ($didx, $thisline); my ($lambda, $inten); my ($maxinten);

# Initialize. $is_ok = 0;

# Get args. $data_p = $_[0]; $lines_p = $_[1];

# If we have args, proceed. if ((defined $data_p) && (defined $lines_p)) {   # Scan the entire file, looking for valid-seeming lines. for ($didx = 0; defined ($thisline = $$data_p[$didx]); $didx++) {     # FIXME - Assuming a specific format! # Column 1 is ionization state/label, column 2 is measured # wavelength, column 3 is Ritz wavelength, column 4 is relative # intensity.

# Store defaults. $lambda = undef; $inten = $inten_bogus;

# First choice: Ritz wavelength exists, intensity exists. if ($thisline =~ m/^[^|]+\|[^|]+\|\s+([\d\.]+)[^|]+\|\s+(\d+)/) {       $lambda = $1; $inten = $2; }     # Second choice: Ritz wavelength exists, but no intensity. elsif ($thisline =~ m/^[^|]+\|[^|]+\|\s+([\d\.]+)/) {       $lambda = $1; }

# Ignore other cases. Ritz wavelength is always in the table, # so no need to check measured wavelength.

# If wavelength is defined, we have a valid data tuple. Store it. # If we found data on any line, report success. if (defined $lambda) {       # FIXME - Culling data that's out of range. # Otherwise we end up messing with normalization. if (($lambda >= $lambda_min) && ($lambda <= $lambda_max)) {         $$lines_p{$lambda} = $inten; $is_ok = 1; }     }    }

# If we're ok, normalize intensities to 1.0 max. if ($is_ok) {     $maxinten = 0; foreach $lambda (keys %$lines_p) {       $inten = $$lines_p{$lambda}; if ($maxinten < $inten) { $maxinten = $inten; } }

if (1.0e-20 < $maxinten) {       foreach $lambda (keys %$lines_p) {         $$lines_p{$lambda} /= $maxinten; }     }    }  }

# Return the resulting error flag. return $is_ok; }


 * 1) Performs a diagnostics readout of spectrum data.
 * 2) Arg 1 points to a hash containing spectral line and intensity data.
 * 3) No return value.

sub DebugLines { my ($lines_p); my ($lambda, $inten);

$lines_p = $_[0];

if (defined $lines_p) {   print STDERR "Spectrum data:\n";

foreach $lambda (sort keys %$lines_p) {     $inten = $$lines_p{$lambda};

if ($inten == $inten_bogus) { $inten = '???'; }

print STDERR " $lambda :   $inten\n"; }

print STDERR "End of spectrum data.\n"; } }


 * 1) Converts a (0..255) RGB integer tuple to a hex tuple.
 * 2) Scales by the supplied intensity, using either liner or log scale.
 * 3) Arg 1 points to a hash containing red, green, and blue components.
 * 4) Arg 2 is the intensity to scale by (0..1).
 * 5) Arg 3 is 'linear' or 'logarithmic'.
 * 6) Returns a 6-character hex string representing the colour.

sub RGBToHex { my ($rgb_p, $inten, $scalemode); my ($result); my ($rval, $gval, $bval); my ($is_log);

# Initialize. $result = "000000";

# Get args. $rgb_p = $_[0]; $inten = $_[1]; $scalemode = $_[2];

# Proceed if args are valid. if ((defined $rgb_p) && (defined $inten) && (defined $scalemode)) {   # Get derived inputs.

$rval = $$rgb_p{red}; $gval = $$rgb_p{green}; $bval = $$rgb_p{blue};

$is_log = 1; if ('linear' eq $scalemode) { $is_log = 0; }

# Clip intensity, and convert to log if appropriate.

if (1.0 < $inten) { $inten = 1.0; } elsif (0.0 > $inten) { $inten = 0.0; }

if ($is_log) {     if ($inten > 1.0e-20) {       $inten = log($inten); # Now -inf..0, base e.        $inten /= log(10); # Now -inf..0, base 10. $inten /= $inten_log_decades; # Now -inf..0, base 10^decades.

$inten += 1; if ($inten < 0.0) { $inten = 0.0; } # Now 0..1 again, but log-scale. }   }

# Scale colour values, and turn into hex.

$rval *= $inten; $gval *= $inten; $bval *= $inten;

$result = sprintf('%02x', $rval). sprintf('%02x', $gval) . sprintf('%02x', $bval); }

# Return the resulting colour. return $result; }


 * 1) Draws a SVG-format spectrum plot from supplied spectrum data.
 * 2) Arg 1 is the name of the file to write to.
 * 3) Arg 2 points to a hash containing line and intensity data.
 * 4) Returns 1 if successful, and 0 if not.

sub MakeSpectrumSVG { my ($oname, $lines_p); my ($is_ok); my ($lambda, $inten); my ($xofs, $yofs); my (%rgb); my ($lidx, $lambda2); my ($hex1, $hex2); my ($fsize);

# Initialize. $is_ok = 0;

# Get args. $oname = $_[0]; $lines_p = $_[1];

# If we have args, proceed. if ((defined $oname) && (defined $lines_p)) {   # Try to open the SVG file. if (!open(OUTFILE, ">$oname")) {     print STDERR "### Unable to write to \"$oname\".\n"; }   else {     #      # Header.

print OUTFILE << "Endofblock"; 



Endofblock


 * 1)  

#     # Write spectrum bar.

$xofs = $img_spacing_big - $lambda_min; $yofs = 3 * $img_spacing_big + $img_spacing_small;

for ($lidx = 0;       defined ($lambda2 = $specbar_lambdas[$lidx + 1]);        $lidx++) {       $lambda = $specbar_lambdas[$lidx];

# Force lambda < lambda2. # Which is defaul depends on the bar list's order. if ($lambda > $lambda2) {         my ($scratch); $scratch = $lambda; $lambda = $lambda2; $lambda2 = $scratch; }

# Clip to range, or cull, if necessary. if (($lambda <= $lambda_max) && ($lambda2 >= $lambda_min)) {         # We haven't been culled. May still have to clip. if ($lambda < $lambda_min) { $lambda = $lambda_min; } if ($lambda2 > $lambda_max) { $lambda2 = $lambda_max; }

# Clipping done. Continue.

# Get colours. %rgb = WaveToRGB($lambda); $hex1 = RGBToHex(\%rgb, 1.0, 'linear'); %rgb = WaveToRGB($lambda2); $hex2 = RGBToHex(\%rgb, 1.0, 'linear');

# Gradient definition. # FIXME - Might not like multiple defs tags? print OUTFILE << "Endofblock";

   

Endofblock

# Rectangle. print OUTFILE ''. "\n"; }     }

# Indicate the boundaries of the visible spectrum, if in range.

if (($vis_min <= $lambda_max) && ($vis_min >= $lambda_min)) {       print OUTFILE ''. "\n"; }

if (($vis_max <= $lambda_max) && ($vis_max >= $lambda_min)) {       print OUTFILE ''. "\n"; }

#     # Write annotating text.

# FIXME - NYI.

# "UV" and "IR" indicators. # These can be on either side, or absent, depending on range!

%rgb = WaveToRGB($lambda_min_min); $hex1 = RGBToHex(\%rgb, 1.0, 'linear'); %rgb = WaveToRGB($lambda_max_max); $hex2 = RGBToHex(\%rgb, 1.0, 'linear');

$xofs = $img_spacing_small; $yofs = 3 * $img_spacing_big + 2 * $img_spacing_small;

$fsize = $img_spacing_small; $fsize .= "pt";

if ($lambda_min < $vis_min) {       print OUTFILE << "Endofblock";

UV Endofblock }

if ($lambda_min > $vis_max) {       print OUTFILE << "Endofblock";

IR Endofblock }

$xofs = $img_spacing_big + ($lambda_max - $lambda_min) + 0.5 * $img_spacing_small;

if ($lambda_max < $vis_min) {       print OUTFILE << "Endofblock";

UV Endofblock }

if ($lambda_max > $vis_max) {       print OUTFILE << "Endofblock";

IR Endofblock }

#     # Write spectral lines.

foreach $lambda (sort keys %$lines_p) {       $inten = $$lines_p{$lambda};

# Translate bogus intensities to minimum, if asked to do so. # FIXME - NYI.

# Avoid bogus intensities. # Also range-check. if (($inten != $inten_bogus)         && ($lambda >= $lambda_min) && ($lambda <= $lambda_max)) {         # Use lines, by default. # FIXME: Draw thick boxes if asked to do so.

%rgb = WaveToRGB($lambda); $xofs = $img_spacing_big - $lambda_min;

# Linear copy. $yofs = $img_spacing_big; print OUTFILE '<line x1="' . ($lambda + $xofs)           . '" y1="' . $yofs . '" x2="' . ($lambda + $xofs)            . '" y2="' . ($yofs + 2 * $img_spacing_big)            . '" style="stroke:#' . RGBToHex(\%rgb, $inten, 'linear')            . ';"/>'. "\n";

# Log copy. $yofs = 3 * $img_spacing_big + 3 * $img_spacing_small; print OUTFILE '<line x1="' . ($lambda + $xofs)           . '" y1="' . $yofs . '" x2="' . ($lambda + $xofs)            . '" y2="' . ($yofs + 2 * $img_spacing_big)            . '" style="stroke:#' . RGBToHex(\%rgb, $inten, 'logarithmic')            . ';"/>'. "\n"; }     }

#     # Footer.

print OUTFILE << "Endofblock";

Endofblock

# Close the SVG file. close(OUTFILE); } }

# Return the resulting error flag. return $is_ok; }


 * 1) Updates image dimensions/spacing to reflect flags.
 * 2) No arguments.
 * 3) No return value.

sub UpdateImageSizes { my ($delta_l);

# NOTE: These can be non-integer values!

# Get wavelength span. $delta_l = $lambda_max - $lambda_min; if ($delta_l < 1.0e-20) { $delta_l = 1; }

# Get spacing lengths, derived from span. $img_spacing_big = $spacing_big_frac * $delta_l; $img_spacing_small = $spacing_small_frac * $delta_l;

# Get image dimensions, derived from spacing and span. $img_width = $delta_l + 2 * $img_spacing_big; $img_height = 6 * $img_spacing_big + 3 * $img_spacing_small; }


 * 1) Main Program
 * 1) Main Program

my ($iname, $oname, %flags); my ($thisarg, $aidx); my (@rawdata, %lines); my ($is_ok);


 * 1) Get args.

$iname = $ARGV[0]; $oname = $ARGV[1];

for ($aidx = 2; defined ($thisarg = $ARGV[$aidx]); $aidx++) { # FIXME - Blithely assuming that all remaining args are valid. # FIXME - Blithely assuming that all flags have no arguments. $flags{$thisarg} = 1; }

if (!((defined $iname) && (defined $oname))) { print << "Endofblock";
 * 1) Check args. Proceed if valid. Print a help screen if not.

Makes an emission spectrum plot from NIST data. Original version written by Christopher Thomas (2009).

Usage: make_spectrum <infile|element> [options]

is a file containing a NIST data table, in text format. is the name of an element to search NIST for. is the name of the SVG image file to write to. [options] are optional flags modifying behavior.

Options are: --fetchnist  Attempts to fetch NIST data for. Default is to read from instead.

Endofblock } else { # Args look ok. Proceed.

# Process flags. # FIXME - NYI.

# Adjust image dimensions/spacing. UpdateImageSizes;

# Fetch raw input data.

@rawdata = ; $is_ok = 0;

if (defined ($flags{'--fetchnist'})) {   $is_ok = GetDataNIST($iname, \@rawdata); } else {   $is_ok = GetDataFile($iname, \@rawdata); }

# If successful, turn raw input data into a hash of lines and # intensities.

%lines = ;

if ($is_ok) {   $is_ok = ExtractLines(\@rawdata, \%lines); }

# If successful, draw the resulting spectrum chart. if ($is_ok) {
 * 1) FIXME - Test.
 * 2) DebugLines(\%lines);

# Ignore return value. MakeSpectrumSVG($oname, \%lines); } }


 * 1) FIXME - Test.
 * 2) DebugRGBVals;


 * 1) Done.


 * 1) This is the end of the file.
 * 1) This is the end of the file.