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)

---

#!/usr/bin/perl
#
# Makes an emission spectrum plot from NIST data.
# Original version written by Christopher Thomas (2009).
#
# Usage: make_spectrum <infile|element> <outfile> [options]
#
# <infile> is a file containing a NIST data table, in text format.
# <element> is the name of an element to search NIST for.
# <outfile> 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 <element>.
#               Default is to read from <infile> 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 = <INFILE>;

      # 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";
<svg xmlns="http://www.w3.org/2000/svg"
  xmlns:xlink="http://www.w3.org/1999/xlink"
  width="$img_width" height="$img_height">

  <rect x="0" y="0" width="$img_width" height="$img_height"
    style="stroke:#000000; fill:#000000"/>

Endofblock

#  <line x1="0" y1="0" x2="100" y2="100" style="stroke:#ff00ff;"/>


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

<defs>
  <linearGradient id="specbar$lidx"
    x1="0%" y1="0%" x2="100%" y2="0%"
    spreadMethod="pad">
    <stop offset="0%" stop-color="#$hex1" stop-opacity="1"/>
    <stop offset="100%" stop-color="#$hex2" stop-opacity="1"/>
  </linearGradient>
</defs>

Endofblock

          # Rectangle.
          print OUTFILE '<rect x="' . ($xofs + $lambda)
            . '" y="' . ($yofs) . '" width="' . ($lambda2 - $lambda)
            . '" height="' . $img_spacing_small
            . '" style="fill:url(#specbar' . $lidx
            . '); stroke:url(#specbar' . $lidx . ');" />' . "\n";
        }
      }

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

      if (($vis_min <= $lambda_max) && ($vis_min >= $lambda_min))
      {
        print OUTFILE '<line x1="' . ($vis_min + $xofs)
          . '" y1="' . $yofs . '" x2="' . ($vis_min + $xofs)
          . '" y2="' . ($img_spacing_small + $yofs)
          . '" style="stroke:#000000;"/>' . "\n";
      }

      if (($vis_max <= $lambda_max) && ($vis_max >= $lambda_min))
      {
        print OUTFILE '<line x1="' . ($vis_max + $xofs)
          . '" y1="' . $yofs . '" x2="' . ($vis_max + $xofs)
          . '" y2="' . ($img_spacing_small + $yofs)
          . '" style="stroke:#000000;"/>' . "\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";

<text x="$xofs" y="$yofs"
  style="fill:#$hex1; font-size:$fsize;">UV</text>
Endofblock
      }

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

<text x="$xofs" y="$yofs"
  style="fill:#$hex2; font-size:$fsize;">IR</text>
Endofblock
      }

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

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

<text x="$xofs" y="$yofs"
  style="fill:#$hex1; font-size:$fsize;">UV</text>
Endofblock
      }

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

<text x="$xofs" y="$yofs"
  style="fill:#$hex2; font-size:$fsize;">IR</text>
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";

</svg>
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 <infile|element> <outfile> [options]

 <infile> is a file containing a NIST data table, in text format.
 <element> is the name of an element to search NIST for.
 <outfile> 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 <element>.
               Default is to read from <infile> 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.
#