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.
#

Content Disclaimer

Informasi ini disarikan dari Wikipedia dan disajikan kembali untuk tujuan edukasi. Konten tersedia di bawah lisensi CC BY-SA 3.0. Kami tidak bertanggung jawab atas ketidakakuratan data yang bersumber dari kontribusi publik tersebut.

  1. The information displayed on this website is sourced in part or in whole from Wikipedia and has been adapted for the purpose of restating it. We strive to provide accurate and relevant information, however:
  2. There is no guarantee of absolute accuracy. Wikipedia is an open, collaborative project that can be edited by anyone, so information is subject to change.
  3. It is not intended to constitute professional advice. The content displayed is for informational and educational purposes only. For important decisions (e.g., medical, legal, or financial), please consult a professional.
  4. Content copyright. Wikipedia is licensed under the Creative Commons Attribution-ShareAlike License (CC BY-SA). This means that content may be reused with appropriate attribution and shared under a similar license.
  5. Responsible use. Any risk arising from the use of information from this website is entirely the responsibility of the user.