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 http://physics.nist.gov/PhysRefData/ASD/lines_form.html] --Christopher Thomas (talk) 17:36, 16 December 2009 (UTC)

----

!/usr/bin/perl

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

Usage: make_spectrum   [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.




Includes


use strict;


Constants



Ultimate limits to wavelength.

my ($lambda_min_min, $lambda_max_max);
$lambda_min_min = 200;
$lambda_max_max = 2000;

Visible spectrum limits.

my ($vis_min, $vis_max);
$vis_min = 380;
$vis_max = 750;

Piecewise-linear spectrum definition.
Roughly matches FIXME:NAME 's mapping function.
Tweaked to remove banding, and to put R/G/B at 675/525/450 nm.
Also tweaked IR and UV values to be obviously-different.

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
);

Spectrum bar landmark wavelengths.

my (@specbar_lambdas);
@specbar_lambdas = (
$lambda_max_max,
750, 725, 675, 600, 525, 500, 450, 425, 400, 375,
$lambda_min_min
);


Configuration Variables



Spacing. Make this a fraction of the wavelength range.

my ($spacing_big_frac, $spacing_small_frac);
$spacing_big_frac = 0.1;
$spacing_small_frac = 0.03;

Wavelength range to query data for.

my ($lambda_min, $lambda_max);
$lambda_min = 200;
$lambda_max = 1000;

Test a smaller range.
$lambda_min = 500;
$lambda_max = 600;


Bogus and default line intensity values.

my ($inten_bogus, $inten_default);
$inten_bogus = 0;
$inten_default = 1;

Actual spacing/size values. These will get reinitialized!

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;

Number of decades to map for log-scale colours.

my ($inten_log_decades);
$inten_log_decades = 2.5;


Functions



Performs a piecewise-linear mapping.
Finds the nearest mapping vertices on either side of the input value,
and interpolates the output value from these.
Arg 1 is the value to be mapped.
Arg 2 points to a hash containing mapping values.
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;
}

Translates a wavelength (in nm) into an RGB value (range 0-255).
Arg 1 is the wavelength.
Returns a hash (by value), containing "red", "green", and "blue"
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;
}

Performs a diagnostics sweep of the spectrum colours.
Dumps output to STDERR.
No arguments.
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";
}
}

Retrieves text-format spectrum data from NIST.
Arg 1 is the name of the element to fetch.
Arg 2 points to an array to store raw data in.
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;
}

Retrieves text-format spectrum data from a file.
Arg 1 is the name of the file to read.
Arg 2 points to an array to store raw data in.
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;
}

Extracts spectrum line data from raw NIST data tables.
Arg 1 points to an array containing the raw NIST data table as text.
Arg 2 points to a hash to store line and intensity data in.
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;
}

Performs a diagnostics readout of spectrum data.
Arg 1 points to a hash containing spectral line and intensity data.
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";
}
}

Converts a (0..255) RGB integer tuple to a hex tuple.
Scales by the supplied intensity, using either liner or log scale.
Arg 1 points to a hash containing red, green, and blue components.
Arg 2 is the intensity to scale by (0..1).
Arg 3 is 'linear' or 'logarithmic'.
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;
}

Draws a SVG-format spectrum plot from supplied spectrum data.
Arg 1 is the name of the file to write to.
Arg 2 points to a hash containing line and intensity data.
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";
xmlns:xlink="http://www.w3.org/1999/xlink"
width="$img_width" height="$img_height">

style="stroke:#000000; fill:#000000"/>
Endofblock



#
# 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";


x1="0%" y1="0%" x2="100%" y2="0%"
spreadMethod="pad">




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";

style="fill:#$hex1; font-size:$fsize;">UV
Endofblock
}
if ($lambda_min > $vis_max)
{
print OUTFILE << "Endofblock";

style="fill:#$hex2; font-size:$fsize;">IR
Endofblock
}
$xofs = $img_spacing_big + ($lambda_max - $lambda_min)
+ 0.5 * $img_spacing_small;
if ($lambda_max < $vis_min)
{
print OUTFILE << "Endofblock";

style="fill:#$hex1; font-size:$fsize;">UV
Endofblock
}
if ($lambda_max > $vis_max)
{
print OUTFILE << "Endofblock";

style="fill:#$hex2; font-size:$fsize;">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 '' . "\n";
# Log copy.
$yofs = 3 * $img_spacing_big + 3 * $img_spacing_small;
print OUTFILE '' . "\n";
}
}
#
# Footer.
print OUTFILE << "Endofblock";

Endofblock
# Close the SVG file.
close(OUTFILE);
}
}
# Return the resulting error flag.
return $is_ok;
}

Updates image dimensions/spacing to reflect flags.
No arguments.
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;
}


Main Program


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

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;
}

Check args. Proceed if valid. Print a help screen if not.

if (!((defined $iname) && (defined $oname)))
{
print << "Endofblock";
Makes an emission spectrum plot from NIST data.
Original version written by Christopher Thomas (2009).
Usage: make_spectrum   [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)
{

FIXME - Test.
DebugLines(\%lines);

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

FIXME - Test.
DebugRGBVals();


Done.



This is the end of the file.