; $Id: sf.pro,v 1.4 1999/03/03 06:14:10 idl Exp $
;
;+
; NAME:
;       SF
;
; LAST CHANGED: ------------------ aph ------------------------ 02-may-04
;
; PURPOSE:
;       Function to calculate x-ray interaction properties of materials.
;
; CATEGORY:
;       X-ray optics
;
; CALLING SEQUENCE:
;       Result = SF(Compound [, Abscissa])
;
; INPUTS:
;       Compound:    A scalar string containing the compound molecular formula.
;                    The structure of the formula must be as follows:
;                      Name1 Amount1 Name2 Amount2 ...
;                    where Name1, Name2, etc. are the one or two letter
;                    standard abbreviations for elements 1 to 92, and
;                    Amount1, Amount2, etc. are the values describing the
;                    stoichiometry. The spaces are optional. Each element may
;                    appear only once in the formula.
;                    NOTES: 1. Formula names ARE case-sensitive
;                              (e.g., NI is NOT Ni).
;                           2. If an amount is not specified, it is assumed
;                              to be 1 (e.g., H2O is the same as H2 O1).
;
; OPTIONAL INPUTS:
;       Abscissa:    A vector of photon energies or wavelengths at which the
;                    result will be evaluated. The type and units of the
;                    abscissae are specified by the keyword parameters. The
;                    default is to assume the abscissae are photon energies
;                    in electron volts.
;
; KEYWORD PARAMETERS:
;       J:           If set, the abscissae, if supplied, are interpreted
;                    as x-ray photon energies in Joules.
;
;       keV:         If set, the abscissae, if supplied, are interpreted
;                    as x-ray photon energies in kiloelectron volts.
;
;       eV:          If set, the abscissae, if supplied, are interpreted
;                    as x-ray photon energies in electron volts.
;
;       m:           If set, the abscissae, if supplied, are interpreted
;                    as x-ray wavelengths in meters.
;
;       nm:          If set, the abscissae, if supplied, are interpreted
;                    as x-ray wavelengths in nanometers.
;
;       Angstrom:    If set, the abscissae, if supplied, are interpreted
;                    as x-ray wavelengths in Angstroms.
;
;       RESULT_TYPE: A string specifying what result to return. The following
;                    are the possible choices:
;                     "f1":      The compound f1 scattering factor for the
;                                molecule.
;
;                     "f2":      The compound f2 scattering factor for the
;                                molecule.
;
;                     "delta":   The real part of the complex refractive index.
;
;                     "beta":    The imaginary part of the complex refractive
;                                index.
;
;                     "n":       The complex refractive index.
;
;                     "mu":      The mass photoabsorption coefficient,
;                                in cm^2/g. This is the default result type.
;
;                     "mu_a":    The atomic photoabsorption coefficient,
;                                in barns/atom.
;
;                     "mu_l":    The linear photoabsorption coefficient,
;                                in cm^-1.
;
;                     "trans":   The transmissivity of the material.
;
;                     "reflect": The specular reflectivity of unpolarized light
;                                on an ideally smooth surface of the material.
;
;       DENSITY:     The density of the material, in g/cm^3. This is needed
;                    if the result type is "delta", "beta", "n", "mu_l",
;                    "trans", or "reflect". Otherwise, a default density
;                    of 1.0 g/cm^3 is used.
;
;       THICKNESS:   The thickness of the material in micrometers. This is
;                    needed if the result type is "trans". Otherwise, a
;                    default thickness of 1.0 micron is used.
;
;       ANGLE:       The angle of incidence in milliradians. This is needed
;                    if the result type is "reflect". Otherwise, a default
;                    angle of incidence of 10.0 mr is used.
;
;       ENERGY:      A variable to receive the x-ray photon energies (in eV)
;                    at which the result was evaluated.
;
;       MW:          A variable to receive the molecular weight (in g/mol) of
;                    the compound.
;
;       SF_DIR:      A string containing the name of the directory containing
;                    the SF data files.
;
; OUTPUTS:
;       Result:      The result specified by the RESULT_TYPE keyword.
;
; COMMON BLOCKS:
;       SF_COMMON
;
; SIDE EFFECTS:
;       SF data is stored into common block variables each time data is read
;       from the data files. This takes up memory to gain faster performance.
;       The memory can be freed using the SF_DESTROY routine.
;
; ROUTINES CONTAINED IN THIS MODULE:
;       SF_INIT           Procedure to initialize common block for module.
;       SF_DESTROY        Procedure to delete common block variables.
;       SF_GEN_EL_DATA    Procedure to generate data file for elements.
;       SF_READ_DATA      Procedure to read SF data for single elements.
;       SF_PARSE_COMPOUND Function to parse molecular formulas.
;       SF                Main function to compute materials x-ray properties.
;
; ROUTINES CALLED BY THIS MODULE:
;       IDL User Library routines:
;         UNIQ, INTERPOL
;       Bill Loo's IDL routines:
;         STRING2FILE, FILE2STRING, STR_COLUMNS, WVLEN2EN, EN2WVLEN
;
; DATA FILES:
;       This module requires the scattering factor data files compiled by
;       B.L. Henke, E.M. Gullikson, and J.C. Davis, available by anonymous
;       ftp at:
;         ftp://www-cxro.lbl.gov/pub/sf/sf.tar.Z
;
;       In addition, it requires a file containing the abbreviated element
;       names and atomic weights. Such a file can be generated manually
;       using the routine SF_GEN_EL_DATA. Just type "SF_GEN_EL_DATA" at the
;       IDL prompt, and enter the requested information.
;
; REFERENCE:
;       All formulas are published in:
;         B.L. Henke, E.M. Gullikson, and J.C. Davis, "X-ray interactions:
;         photoabsorption, scattering, transmission, and reflection at
;         E=50-30,000 eV, Z=1-92," Atomic Data and Nuclear Data Tables,
;         Vol. 54(2), pp. 181-342, July 1993.
;
; WRITTEN BY:
;       Billy W. Loo, Jr.
;       Bioengineering Graduate Group, UCB / UCSF
;       School of Medicine, UCD
;       Ernest Orlando Lawrence Berkeley National Lab
;       August, 1996
;
; MODIFICATION HISTORY:
;
; $Log: sf.pro,v $
; Revision 1.4  1999/03/03 06:14:10  idl
; * Changed function names for more consistent naming.
; * Minor changes to documentation.
;
; Revision 1.3  1999/02/28 00:26:10  idl
; * Return scalar when result has only one element.
;
; Revision 1.2  1997/12/16 01:55:20  loo
; * Major revision of SF program module.

; (02-may-04 aph) standardize header
;-

;=============================================================================
;=           Procedure to initialize common block for module                 =
;=============================================================================

PRO SF_INIT, sf_dir

  ;Declare common block for module.
  COMMON SF_COMMON, sf_dir_,   $ ;Directory containing SF data files.
                    sf_ext_,   $ ;File extension of SF data files.
                    el_names_, $ ;Array of element names.
                    at_wts_,   $ ;Array of element atomic weights.
                    energy_,   $ ;Array of pointers to energy for each element.
                    f1_,       $ ;Array of pointers to f1s for each element.
                    f2_          ;Array of pointers to f2s for each element.

  ;Initialize directory containing SF data files.
  IF (N_PARAMS() GT 0) THEN sf_dir_ = sf_dir ELSE sf_dir_ = ''

  ;Default file extension.
  sf_ext_ = '.nff'

  ;Try to load element data.
  IF (N_ELEMENTS(el_names_) EQ 0) THEN BEGIN
    ;Search for data file.
    fn = FINDFILE(sf_dir_ + 'elements.dat')

    IF (fn[0] NE '') THEN BEGIN
      data = FILE2STRING(fn[0])
      STR_COLUMNS, data

      el_names_ = data[0, *]
      at_wts_   = [[0], [FLOAT(data[1, 1:*])]]
    ENDIF $

    ELSE PRINT, 'Data file not found!' + STRING(7B)

    ;Initialize pointer arrays.
    energy_ = PTRARR(93)
    f1_     = PTRARR(93)
    f2_     = PTRARR(93)
  ENDIF

  RETURN
END ;sf_init

;=============================================================================
;=         Procedure to delete common block variables for module             =
;=============================================================================

PRO SF_DESTROY
  ;Declare common block for module.
  COMMON SF_COMMON

  ;Delete all common block variables for module.
  tmp = TEMPORARY(sf_dir_)
  tmp = TEMPORARY(sf_ext_)
  tmp = TEMPORARY(el_names_)
  tmp = TEMPORARY(at_wts_)

  PTR_FREE, energy_
  PTR_FREE, f1_
  PTR_FREE, f2_

  tmp = TEMPORARY(energy_)
  tmp = TEMPORARY(f1_)
  tmp = TEMPORARY(f2_)

  RETURN
END ;SF_DESTROY

;=============================================================================
;=Procedure to generate data file containing element names and atomic weights=
;=============================================================================

PRO SF_GEN_EL_DATA
  ;Declare common block for module.
  COMMON SF_COMMON

  ;Initialize storage.
  result = STRARR(2, 93)
  result[*, 0] = ['Element_name:', 'Atomic_weight:']

  ;User fills in names of elements.
  PRINT, 'Enter element names and atomic weights according to atomic number.'

  FOR i = 1, 92 DO BEGIN
    temp = ''
    READ, temp, PROMPT='Element #' + STRTRIM(i, 2) + ': '
    result[0, i] = temp

    READ, temp, PROMPT='Atomic weight: '
    result[1, i] = temp
  ENDFOR

  el_names_ = result[0, *]
  at_wts_   = [[0.0], [FLOAT(result[1, 1:*])]]

  result = result[0, *] + STRING(9B) + result[1, *]

  STRING2FILE, result, sf_dir_ + 'elements.dat'

  RETURN
END ;sf_gen_el_data

;=============================================================================
;=      Procedure to read scattering factor data for a single element        =
;=============================================================================

PRO SF_READ_DATA, element, energy, f1, f2, SF_DIR=sf_dir
  ;Declare common block for module.
  COMMON SF_COMMON

  ;Check that common block has been initialized. Initialize if required.
  IF (N_ELEMENTS(el_names_) EQ 0) THEN BEGIN
    IF (N_ELEMENTS(sf_dir) NE 0) THEN SF_INIT, sf_dir $
    ELSE SF_INIT
  ENDIF

  ;-------------------- Parse arguments and set defaults --------------------
  ;Return to main level on error
  ;ON_ERROR, 1

  ;Determine which element is desired.
  s = SIZE(element)

  ;If input is string containing element name:
  IF (s[s[0]+1] EQ 7) THEN BEGIN
    ;Convert to lower case.
    e_name = STRLOWCASE(element)

    ;Make sure it's a valid element name
    IF ((WHERE(STRLOWCASE(el_names_) EQ e_name))[0] EQ -1) THEN BEGIN
      MESSAGE, 'Invalid element name.' + STRING(7B)
    ENDIF
  ENDIF $

  ;If input is atomic number of element:
  ELSE BEGIN
    IF (element LT 1) OR (element GT 92) or (element NE FIX(element)) THEN $
      MESSAGE, 'Invalid atomic number.' + STRING(7B)

    ;Determine element name.
    e_name = STRLOWCASE(el_names_[element])
  ENDELSE

  ;-------------------------- Main body of program --------------------------

  ;Read scattering factor file of interest.

  f_name = sf_dir_ + e_name + sf_ext_

  ;Try to find correct data file.
  fname = FINDFILE(f_name)
  IF (fname[0] EQ '') THEN MESSAGE, 'SF data file not found!' + STRING(7B)

  ;Read in data.
  data = FILE2STRING(fname[0])
  STR_COLUMNS, data

  ;Fill in results arrays. First line is header.
  energy = FLOAT(data[0, 1:*])
  f1     = FLOAT(data[1, 1:*])
  f2     = FLOAT(data[2, 1:*])

  RETURN
END ;sf_read_data

;=============================================================================
;=    Function to parse a compound name into elements and stoichiometry      =
;=============================================================================
; Usage: SF_PARSE_COMPOUND, Input
;
; Parameters:
;   Input: a string variable containing the molecular formula of the compound
;          The structure of the formula must be as follows:
;            Name1 Amount1 Name2 Amount2 ...
;          where Name1, Name2, etc. are the one or two letter standard
;          abbreviations for elements 1 to 92, and Amount1, Amount2, etc.
;          are the values describing the stoichiometry. The spaces are
;          optional. Each element may appear only once in the formula.
;          NOTES: 1. Formula names ARE case-sensitive (e.g., NI is NOT Ni)
;                 2. If an amount is not specified, it is assumed to be 1
;                    (e.g., H2O is the same as H2 O1)
;
; Returns: a 2 dimensional array with the following information:
;          Each row is as long as the number of elements in the formula.
;          Row(0) contains the atomic numbers of the elements in the compound.
;          Row(1) contains the corresponding amount of each element.
;=============================================================================

FUNCTION SF_PARSE_COMPOUND, input, SF_DIR=sf_dir
  ;Declare common block for module.
  COMMON SF_COMMON

  ;Check that common block has been initialized. Initialize if required.
  IF (N_ELEMENTS(el_names_) EQ 0) THEN BEGIN
    IF (N_ELEMENTS(sf_dir) NE 0) THEN SF_INIT, sf_dir $
    ELSE SF_INIT
  ENDIF

  ;Return to main level on error.
  ;ON_ERROR, 1

  ;Determine which elements have names with one character
  ;and which have two.
  lengths = STRLEN(el_names_)

  ones = WHERE(lengths EQ 1, count1)
  twos = WHERE(lengths EQ 2, count2)

  ;Remove all whitespace from compound name
  compound = STRCOMPRESS(input, /REMOVE_ALL)

  ;Create variables to store element number and position in the string.
  el_num  = 0
  str_pos = 0

  ;First search for presence of two letter elements.
  FOR i = 0, count2-1 DO BEGIN
    pos = STRPOS(compound, el_names_[twos[i]])

    ;If the element is in the compound...
    IF (pos NE -1) THEN BEGIN
      ;Record element number and position.
      el_num  = [el_num, twos[i]]
      str_pos = [str_pos, pos]

      ;Replace corresponding position in compound with two white spaces.
      STRPUT, compound, '  ', pos
    ENDIF
  ENDFOR

  ;Now search for presence of one letter elements.
  FOR i = 0, count1-1 DO BEGIN
    pos = STRPOS(compound, el_names_[ones[i]])

    ;If the element is in the compound...
    IF (pos NE -1) THEN BEGIN
      ;Record element number and position.
      el_num  = [el_num, ones[i]]
      str_pos = [str_pos, pos]

      ;Replace corresponding position in compound with one white space.
      STRPUT, compound, ' ', pos
    ENDIF
  ENDFOR

  ;If no compounds are found, abort.
  length  = N_ELEMENTS(el_num)
  IF (length EQ 1) THEN MESSAGE, 'Invalid compound.' + STRING(7B)

  ;Strip of leading dummy value.
  el_num  = el_num[1:length-1]
  str_pos = str_pos[1:length-1]

  ;If any non-numeric characters are left, abort.
  ;Convert to byte for testing.
  test     = BYTE(' .0123456789') ;Space and numerics only.
  compound = BYTE(compound)

  FOR i = 0, N_ELEMENTS(compound)-1 DO BEGIN
    pos = WHERE(test EQ compound[i])

    IF (pos[0] EQ -1) THEN MESSAGE, 'Invalid compound.' + STRING(7B)
  ENDFOR

  ;Now find numbers indicating stoichiometry.
  ;Note: at this point, all numbers must be separated by spaces.

  ;Variable to store stoichiometry values (as strings).
  ;Note: if unspecified, stoichiometry for an element is 1.
  values = REPLICATE('1.0', length-1)

  ;Find all non-space characters in string (space is ASCII character number 32).
  pos = WHERE(compound NE 32, count)

  ;Convert back to string.
  compound = STRING(compound)

  IF (pos[0] NE -1) THEN BEGIN
    i = 0

    ;Step through characters.
    WHILE (i LT count) DO BEGIN
      ;Number of characters in stoichiometry value.
      length = 1

      ;Find element associated with the stoichiometry value.
      el_pos = WHERE(str_pos LT pos[i])

      ;Check validity.
      IF (el_pos[0] EQ -1) THEN MESSAGE, 'Invalid compound.' + STRING(7B)

      el_pos = MAX(str_pos[el_pos])
      el_pos = WHERE(str_pos EQ el_pos)

      start = pos[i]

      ;Skip over consecutive positions to reach next number.
      WHILE (pos[(i+1) < (count-1)] EQ (pos[i]+1)) DO BEGIN
        length = length + 1
        i      = i + 1
      ENDWHILE

      values[el_pos] = STRMID(compound, start, length)

      i = i + 1
    ENDWHILE
  ENDIF

  ;Convert stoichiometry string values to numbers.
  stoich = FLTARR(N_ELEMENTS(el_num))

  ON_IOERROR, Invalid
  success = 0
  READS, values, stoich
  success = 1

  Invalid: IF NOT(success) THEN MESSAGE, 'Invalid compound.' + STRING(7B)

  ;Resort arrays according to ascending element number.
  st = SORT(el_num)

  RETURN, [[FLOAT(el_num[st])], [stoich[st]]]
END ;sf_parse_compound

;=============================================================================
;=        Main function to compute various SF data for compounds             =
;=============================================================================

FUNCTION SF, compound, abscissa, J=j, KEV=kev, EV=ev,            $
                                 M=m, NM=nm, ANGSTROM=angstrom,  $
                                 RESULT_TYPE=result_type,        $
                                 DENSITY=density, ANGLE=angle,   $
                                 THICKNESS=thickness,            $
                                 ENERGY=energy, MW=mw, SF_DIR=sf_dir
  ;Declare common block for module.
  COMMON SF_COMMON

  ;Check that common block has been initialized. Initialize if required.
  IF (N_ELEMENTS(el_names_) EQ 0) THEN BEGIN
    IF (N_ELEMENTS(sf_dir) NE 0) THEN SF_INIT, sf_dir $
    ELSE SF_INIT
  ENDIF

  ;Define some useful constants.
  Na   = 6.0220453E23  ;Avogadro's number.
  barn = 1E-24         ;(cm^2 / barn) conversion factor.
  r0   = 2.8179380E-15 ;(m) Classical electron radius.

  ;Default density is 1.0 g/cm^3.
  IF (N_ELEMENTS(density) EQ 0) THEN density = 1.0

  ;Default angle of incidence is 10 milliradians.
  IF (N_ELEMENTS(angle) EQ 0) THEN angle = 10.0

  ;Default sample thickness is 1.0 micrometers.
  IF (N_ELEMENTS(thickness) EQ 0) THEN thickness = 1.0

  ;Parse compound into elements and stoichiometries.
  comp = SF_PARSE_COMPOUND(compound)
  el  = comp[*, 0]
  amt = (TEMPORARY(comp))[*, 1]

  ;Number of elements in compound.
  num_els = N_ELEMENTS(el)

  ;Calculate molecular weight.
  MW = TOTAL(amt * at_wts_[el])

  ;Read SF data for each element, if not already in memory.
  FOR i = 0, num_els-1 DO BEGIN
    IF NOT(PTR_VALID(energy_[el[i]])) THEN BEGIN
      energy_[el[i]] = PTR_NEW(/ALLOCATE_HEAP)
      f1_[el[i]]     = PTR_NEW(/ALLOCATE_HEAP)
      f2_[el[i]]     = PTR_NEW(/ALLOCATE_HEAP)

      SF_READ_DATA, el[i], *energy_[el[i]], *f1_[el[i]], *f2_[el[i]]
    ENDIF
  ENDFOR

  ;Convert abscissa to energy in eV.
  IF (N_PARAMS() GT 1) THEN BEGIN
    IF KEYWORD_SET(m)        THEN energy = WVLEN2EN(abscissa, /m, /eV)
    IF KEYWORD_SET(nm)       THEN energy = WVLEN2EN(abscissa, /nm, /eV)
    IF KEYWORD_SET(Angstrom) THEN energy = WVLEN2EN(abscissa, /Angstrom, /eV)

    IF KEYWORD_SET(J)   THEN energy = abscissa / 1.6021892E-19 ;(J/eV)
    IF KEYWORD_SET(keV) THEN energy = abscissa * 1000.0        ;(eV/keV)
    IF KEYWORD_SET(eV)  THEN energy = abscissa

    ;Check if energy is already in eV.
    IF (N_ELEMENTS(energy) EQ 0) THEN energy = abscissa
  ENDIF $

  ;If no abscissae are specified, use whole range of tabulated energies.
  ELSE BEGIN
    energy = 0.0
    FOR i = 0, num_els-1 DO energy = [[energy], [*energy_[el[i]]]]

    energy = energy[1:*]
    energy = energy[UNIQ(energy, SORT(energy))]
  ENDELSE

  ;Convert to log values for interpolation.
  energy = ALOG(TEMPORARY(energy))

  ;Check result type.
  IF (N_ELEMENTS(result_type) EQ 0) THEN res_type = 'mu' $
  ELSE res_type = STRLOWCASE(result_type)

  ;Check whether f1 and f2 both need to be evaluated.
  IF (res_type EQ 'f1')      OR  $
     (res_type EQ 'delta')   OR  $
     (res_type EQ 'n')       OR  $
     (res_type EQ 'reflect') THEN f1_needed = 1
  IF (res_type NE 'f1')    AND $
     (res_type NE 'delta') THEN f2_needed = 1

  ;Extract compound f1 and f2 values for specified energies.
  IF KEYWORD_SET(f1_needed) THEN f1 = FLTARR(N_ELEMENTS(energy))
  IF KEYWORD_SET(f2_needed) THEN f2 = FLTARR(N_ELEMENTS(energy))

  FOR i = 0, num_els-1 DO BEGIN
    tmp_en = ALOG(*energy_[el[i]])
    IF KEYWORD_SET(f1_needed) THEN $
      f1 = TEMPORARY(f1) + amt[i] * INTERPOL(*f1_[el[i]], tmp_en, energy)
    IF KEYWORD_SET(f2_needed) THEN $
      f2 = TEMPORARY(f2) + amt[i] * INTERPOL(*f2_[el[i]], tmp_en, energy)
  ENDFOR

  ;Convert single element arrays to scalars.
  IF KEYWORD_SET(f1_needed) AND (N_ELEMENTS(f1) EQ 1) THEN f1 = f1[0]
  IF KEYWORD_SET(f2_needed) AND (N_ELEMENTS(f2) EQ 1) THEN f2 = f2[0]

  ;Convert energy back to eV.
  energy = EXP(TEMPORARY(energy))

  ;Just return compound f1 or f2, if specified.
  IF (res_type EQ 'f1') THEN RETURN, f1
  IF (res_type EQ 'f2') THEN RETURN, f2

  ;Get wavelength in meters.
  wvln = EN2WVLEN(energy, /eV, /m)

  ;Compute real part of refractive index, if required.
  IF (res_type EQ 'delta')   OR $
     (res_type EQ 'n')       OR $
     (res_type EQ 'reflect') THEN BEGIN
    ;(1 = 1E6 cm^3/m^3)
    delta = (r0 / (2.0 * !PI)) * (density * Na * 1E6 / MW) * wvln^2.0 * f1

    ;Return real part of refractive index, if specified.
    IF (res_type EQ 'delta') THEN RETURN, delta
  ENDIF

  ;Compute imaginary part of refractive index, if required.
  IF (res_type EQ 'beta')    OR $
     (res_type EQ 'n')       OR $
     (res_type EQ 'reflect') THEN BEGIN
    ;(1 = 1E6 cm^3/m^3)
    beta = (r0 / (2.0 * !PI)) * (density * Na * 1E6 / MW) * wvln^2.0 * f2

    ;Return real part of refractive index, if specified.
    IF (res_type EQ 'beta') THEN RETURN, beta
  ENDIF

  ;Compute and return refractive index, if specified.
  IF (res_type EQ 'n') THEN RETURN, COMPLEX(1.0 - delta, -beta)

  ;Compute and return specular reflectivity, if specified.
  IF (res_type EQ 'reflect') THEN BEGIN
    theta = angle * 1E-3 ;radian/milliradian

    rhosq = (1.0/2.0) *                   $
            (sin(theta)^2.0 - 2.0*delta + $
             SQRT((sin(theta)^2.0 - 2.0*delta)^2.0 + 4.0*beta^2.0))
    rho = SQRT(rhosq)

    Isig = (rhosq * (sin(theta) - rho)^2.0 + beta^2.0) / $
           (rhosq * (sin(theta) + rho)^2.0 + beta^2.0)

    Ipisig = (rhosq * (rho - cos(theta) / tan(theta))^2.0 + beta^2.0) / $
             (rhosq * (rho + cos(theta) / tan(theta))^2.0 + beta^2.0)

    RETURN, Isig * (1 + Ipisig) / 2.0
  ENDIF

  ;Compute atomic photoabsorption coefficient (barns/atom), if required.
  IF (res_type EQ 'mu_a')  OR $
     (res_type EQ 'mu_l')  OR $
     (res_type EQ 'mu')    OR $
     (res_type EQ 'trans') THEN BEGIN

    mu_a = f2 * wvln * (2.0 * r0 * 1E28) ;(1 = 1E28 barn/m^2)

    ;Return atomic photoabsorption coefficient, if specified.
    IF (res_type EQ 'mu_a') THEN RETURN, mu_a
  ENDIF

  ;Compute mass photoabsorption coefficient (cm^2/g), if required.
  IF (res_type EQ 'mu_l')  OR $
     (res_type EQ 'mu')    OR $
     (res_type EQ 'trans') THEN BEGIN
    mu = mu_a * (Na * 1e-24 / MW) ;(1 = 1E-24 cm^2/barn)

    ;Return mass photoabsorption coefficient, if specified.
    IF (res_type EQ 'mu') THEN RETURN, mu
  ENDIF

  ;Compute and linear photoabsorption coefficient (cm^-1), if required.
  IF (res_type EQ 'mu_l')  OR $
     (res_type EQ 'trans') THEN BEGIN
    mu_l = mu * density

    ;Return linear photoabsorption coefficient, if specified.
    IF (res_type EQ 'mu_l') THEN RETURN, mu * density
  ENDIF

  ;Compute and return transmissivity, if specified.
  ;(1 = 1E-4 cm/micron)
  IF (res_type EQ 'trans') THEN RETURN, EXP(-mu_l * thickness * 1E-4)

  ;If result has not been returned, result type is invalid.
  MESSAGE, 'Unrecognized result type.' + STRING(7B)
END ;sf
