Users guide, Part I: Optics programs

2. The Kramers-Kronig Toolbox

V. Ossenkopf, MPG AG “Dust in Star-Forming Regions” Jena

Abstract

The Kramers-Kronig toolbox (KKTOOL) is an interactive program for manipulating optical constants/dielectric functions of different materials. At any given step the Kramers-Kronig consistency of the materials can be checked or reproduced. Starting from an arbitrary set of constants (n, k, ϵ′, or ϵ′′) the other functions may be produced by a Kramers-Kronig analysis. A convenient extrapolation of the constants outside of the range covered by the material data has to be done by the user. Several standard functions for the extrapolations are available and the continuation can be visually controlled. The computed constants can be stored or used in further calculations. All computations can be controlled by eye via graphical plots.

1 Basic algorithms

Starting from an arbitrary set of optical constants from the pairs n and k or ϵ′ and ϵ′′ the second set of the pair may be reproduced from the Kramers-Kronig integrals

ϵ′(λ) − 1 = \frac{2}{π} P ∫_{0}^{∞} \frac{λ^{2}}{Λ} \frac{ϵ′′(Λ)}{λ^{2} − Λ^{2}} dΛ  (1)

and ϵ′′(λ) = −\frac{2}{π} P ∫_{0}^{∞} \frac{λ(ϵ′(Λ) − 1)}{λ^{2} − Λ^{2}}dΛ  (2)

The same equations written here for ϵ also hold for the constants n and k where one has simply to substitute ϵ′ by n and ϵ′′ by k.

Two severe problems limit the applicability of this formula: Optical data are never given at all wavelengths between 0 and ∞, and there is no numerical implementation of a principal value of an integral around a singularity with arbitrary accuracy.

The first problem asks for an extrapolation of the data outside of the given range. Here, the user has to give analytic functions. In case of power-law extrapolations with integer exponents, the contributions are computed analytically. Otherwise, they are computed numerically using the Simpson rule. As criterion for the dynamical stepwidth and the most outer interval computed, a relative accuracy of all contributions of 10^{−7} is required.

The singularity in the principal value integral may be circumvented by leaving out a small interval around the critical wavelength. This interval gives no contribution to the integral when the interval is symmetric to the singularity and the denominator function is constant within it. The interval omitted by the used integrator is determined by the condition that the optical data must be constant within 10^{−4} of its value in it.

In the range of given data values the set of optical data is approximated by cubic splines. Then, the Kramers-Kronig integral may be computed analytically for each interval between two wavelength points. This approach provides a much better accuracy than a direct cubic spline integration of the whole integrator function since it better fits the behaviour near the singularities.

The other pair of constants is reproduced from

n = \sqrt{1/2 (|ϵ| + ϵ′)}
k = \sqrt{1/2 (|ϵ| − ϵ′)}  (3)
ϵ′ = n^{2} − k^{2}
ϵ′′ = 2nk .  (4)

2 User interface

The program is designed for a completely interactive data processing. All parts contain menus which show the operations available in each status of the program. The menus follow a hierarchical structure. A step back to the next upper menu is always possible. Some points in a menu can only be reached after the definition of starting parameters. In case of operations for single materials the selection of a ”0.” material leads back to the last menu. All materials and menu functions are represented by integer numbers.

2.1 Menu structure

Main menu:

                           Kramers-Kronig analysis program
Material 1 :             0.000E+00 <lambda< 0.000E+00
Material 2 :             0.000E+00 <lambda< 0.000E+00
Material 3 :             0.000E+00 <lambda< 0.000E+00
Material 4 :             0.000E+00 <lambda< 0.000E+00
Material 5 :             0.000E+00 <lambda< 0.000E+00
Material 6 :             0.000E+00 <lambda< 0.000E+00

     (1) Read optical constants                       (2) Save optical constants
     (3) Calculation / Modification                   (4) Combination of media
     (5) Graphic                                      (6) Quit

From this menu one can call the basic routines for the file input and output of optical data. The combination of media allows to perform a weighted sum up of the optical constants of several materials. As there is no deeper physical meaning in this procedure one has to take care when to apply it. This menu point does not contain an effective media calculator. After combining the optical constants of different media, the Kramers-Kronig consistency of the resulting material has to be reproduced by menu point (3). The graphics routine called by menu point (5) shows one set of constants (either n, k, ϵ′, or ϵ′′) depending on the wavelength λ for up to 4 materials. The calculation/modification menu contains the preparation and execution of the Kramers-Kronig analysis. At the beginning the user has to decide which set of constants from which material is to be used for the integration. Then the menu opens up.

Parameter modification:

     (1) Short wavelength extrapolation (2) Long wavelength extrapolation
     (3) Select special wavelength modus(4) Calculation of new constants
     (5) Extrapolation grafic           (6) Main menu

Before the Kramers-Kronig integral can be computed, one has to define functions for the extrapolation of the set of optical constants outside of the range covered by the given data. When an appropriate extrapolation is found, it can also be used directly for other (of course only similar) materials. The extrapolation graphic allows an visual control whether smooth transitions between the given data and the extrapolating functions occur and the whole behaviour seems physical. Various functions may be chosen for the extrapolation at the two wavelength ends.

Extrapolation menus:

      Possible kinds of short wavelength extrapolation:
      (1) Lorentz/Drude in lambda->0 limit (2 parameters)
      (2) a*lambda^n+b behaviour           (3 parameters)
      (3) a*exp(-1/lambda)+b               (2 parameters)

      Possible kinds of long wavelength extrapolation:
      (1) Lorentz in lambda->inf.  limit (1 parameter)
      (2) a*lambda^(-n)+b behaviour      (3 parameters)
      (3) a*exp(-lambda)+b               (2 parameters)

The Lorentz behaviour is a special case of the power low extrapolation in the menu points (2). In the short wavelength range we have aλ^{2}+b for n and ϵ′ and aλ^{3} for k and ϵ′′. In both cases the contributions from wavelengths much shorter than the interval considered to the real constants n or ϵ′ has to be given. In the long wavelength extrapolation the Lorentz behaviour is according to aλ^{0} for n and ϵ′ and aλ^{−1} for k and ϵ′′. Chosing a power law behaviour for the extrapolation the exponent must be above 0 for the short wavelength range and below 1 or 2 for the real and imaginary constants, respectively, in the long wavelength range. The highest exponents have to be used with care, however, due to an unphysical static behaviour of the such materials (Purcell 1969). The number of parameters which have to be given by the user is less by one than the number of parameters required for the extrapolation formula since this parameter can be determined from the first or last value of the given constants.

When the extrapolation at both wavelength ends is defined, the new constants can be computed from the parameter modification menu. Three separated routines are called for the Kramers-Kronig integration, responsible for the two extrapolated regions and the range of given data.

Instead of the computation of the new set of optical constants at the same wavelength points as in the given set, one can chose another wavelength scaling in the wavelength modus menu:

Wavelengths modus:

     Use a given wavelengths scale (1)                 (2) Main menu
     Use a linear scaling :        (3) in wavelengths  (4) in wavenumbers
     Use a logarithmic scaling :   (5) in wavelengths  (6) in wavenumbers

Here, one can select the best wavelength scaling for the new material. In case of linearily or logarithmically equidistant points one can chose the starting point, the stepwidth and the number of steps. The selection of a special wavelength modus can be dangerous because the strict separation of data and extrapolation range is lost and one might forget that the extrapolated data are no real measured data.

2.2 Data files

All data files treated as input or output of the program have the same structure. They are ASCII files containing lines with three numbers separated by blanks. Every line represents one wavelength point. The files have to consist of at least 2 lines. Four different file types are recognized automatically according to their file extensions. The first letter of a file extension states whether the first column contains wavelengths (in cm) (letter “l”) or wavenumbers (in cm^{−1}) (letter “z”). The next two letters stand for the content of the second and third column. An “nk” states that they have to be considered as optical constants n and k and “ee” leads to an interpretation of the data as real and imaginary part of the dielectric function ϵ′ and ϵ′′. All files with other extensions are rejected by the program.

3 System requirements

The program is written in standard FORTRAN 77 so that it can be compiled on all machines providing an f77 compiler. Due to the lack of standardized graphics routines an adaption to the graphics libraries available on different machines is necessary. This adaption has to be done in the minimal graphics library of the program (GRAPHMIN) which provides all necessary graphics routines based on simple functions as line drawing or label setting. Furthermore, this file contains the system dependent definition of SMLPAUSE which provides a delay of 3 seconds used to underline different warning messages. There are no system dependent functions and subroutines outside of GRAPHMIN ! You might only want to change the path delimiter symbol in NAMECHEC for the DOS version from ‘/’ to ‘álthough this is not necessary.

Uptil now, there exist versions of GRAPHMIN using the IBM Professional Fortran graphics library for PCs, the Microsoft Fortran graphics library for PCs and the Mongo graphics subroutines for SUNs and VAXes. The programs based on the MONGO graphics functions may be executed within gterm (SunView) or within xterm (XWindows, OpenWindows). These terminal programs have to be able to create a Tektronix4101 graphics window and to close it after the graphics channel is closed.

The program does not use large data arrays so that the storage requirement is mainly dependent on the compiler efficiency and the graphics routines linked. The running size is between 170 and 200 kB for the PC versions and 680 kB for the SUN versions.

4 Local particularities - Jena

All PC versions available on the net are compiled to use a numerical coprocessor. On machines without physical coprocessor it may be emulated by Franke or other software coprocessors. On the SUNs the gterm should be called with the options gterm -G -Gopen show -Gclose blank kktool to enable a correct handling of the graphics window and to start the Kramers-Kronig toolbox directly. The xterm available on the Alpha machines does not close the graphics window after its release so that it has to be closed by hand.

5 Limitations and bugs

In this version the number of wavelengths to be treated for each material is restricted to 256. At most six materials, can be loaded into the program. The files containing the optical data of the different materials have to follow the conventions for names and structures discussed above. Different file structures may cause abnormal program termination.

The integration routine sometimes selects inappropriate cubic splines at the ends of the interval of given optical constants when the extrapolation has a strongly different slope. This may result in wrong values for the optical constants at the two most outer points of the given data interval. Inner data points are not influenced.

6 Description of the different parts of the program

6.1 Main program (KKTOOL, HEADMENU, INITDATA)

These routines provide the outer surface of the program. They show the most important parameters of all sets of optical constants loaded into the program, initialize the data at the start of the program and set up the main menu for the data handling.

HEADMENU is the routine showing the parameters of all materials loaded. It is called by all selection menus so that they contain always a survey on the most recent parameters. INITDATA resets the names and wavelength ranges of all materials at the program start. KKTOOL contains the main menu allowing the selection of the file operations, material combination, graphics functions, and the preparation of the Kramers-Kronig analysis.

6.2 File operations (READCNST, READFILE, SAVECNST, SAVEFILE, NAMECHEC)

The file operations treat always one material which may be given as optical constants or dielectric functions versus wavelength (in cm) or wavenumber (in cm^{−1}). The files have to consist of lines containing three ASCII numbers separated by blanks. The name of a data file is used to denote the material in the menus. The extension determines which kind of data shall be written to the file or how the data lines in the file have to be interpreted: lnk - wavelength n k znk - wavenumber n k lee - wavelength ϵ′ ϵ′′. zee - wavenumber ϵ′ ϵ′′. All files with other extensions are rejected.

The read and write of data files contains two steps. READCNST and SAVECNST are interface procedures that ask for the component to be saved, that select all constants of this component, that extract the minimum and maximum wavelength from the data, and that organize the data in an ascending order with wavelength. They produce the dielectric functions ϵ from the optical constants n and k. The routines READFILE and SAVEFILE use only one field of optical constants. They ask for a filename and transform the data depending on the NAMECHEC result into wavelength and n and k or the format required for storage. They are the procedures reading and writing the data on the disks. NAMECHEC analyses the filename given by the user, removes the path from it and extracts its extension which is translated into the NAMECHEC result.

6.3 Graphics (GRAPHIC, EXGRAPH , PLOT, GRAPHMIN)

The graphics routines allow to show the sets of optical constants of up to four materials. One set of constants (n, k, ϵ′, or ϵ′′) will be shown. According to experience with plots of optical data, the scales in one direction are chosen to be logarithimally when the data span a range of more than the factor 80. Otherwise, the program uses a linear scale.

The routine GRAPHIC is called from the Main menu. It asks for the materials to be shown and for the set of constants to be plotted. Then it selects those constants for the materials and calls the subroutine PLOT. EXGRAPH is the subroutine for the demonstration of the extrapolation. It is called from the parameter modification menu. It uses the set of constants chosen for the Kramers-Kronig analysis and produces a second data field which extends by a factor 10 to lower and higher wavelengths using the selected extrapolation. With these two data fields it calls PLOT to show the behaviour of th extrapolation graphically. PLOT sets up the graph of up to four functions which are passed to it as parameters. PLOT uses the basic graphics functions provided by GRAPHMIN to set up the plot. After an initialization of the graphics area and the set up of the coordinate system the data given to PLOT are shown as linegraphs. The routine waits for an ⟨Enter⟩ (⟨Return⟩) of the user to show the second graph or to return to the main menu. GRAPHMIN contains the graphics routines: GRINIT initializes the graphics window and sets up the coordinate system with the kind of axis scaling chosen by the calling program, GRDRAW plots linegraphs for single functions, GRAXIS produces the axis labels and uses GRPRINT, GRCLOSE closes the graphics window, WAIT stops program execution until the user presses ⟨Enter⟩ and SMLPAUSE produces a three seconds pause.

6.4 Combination of materials (ADDFILES, VALUEAT)

With these functions it is possible to change the optical constants of given materials by simple algebraic functions. The manipulations may be done either in the optical constants n and k or the dielectric functions ϵ. The second set is reproduced automatically. All functions may be multiplied by a fixed value and the functions of different materials may be summed up. It is, e.g., possible to obtain the arithmetic averages of the cosntants of different materials. The functions produced in this way are, however, in gerneral not Kramers-Kronig consistent so that the Kramers-Kronig integration has to be added.

The first material used in the combination procedure determines the wavelength points where the optical constants are computed. The addition of further materials may only remove some wavelength points which are not covered by the range where the optical constants of these materials are given.

The procedure VALUEAT interpolates the optical constants of the additional materials at the wavelength points given by the first one by a linear function. The procedure ADDFILES asks for the materials to be combined and where to store the new constants and carries out the arithmetic operations. Finally it reproduces the second set of constants.

6.5 Parameter modifications (PARMMENU, TAKEALL)

Before the parameter modifications menu is entered, one has to chose which material and which set of constants shall serve as starting point for the Kramers-Kronig analysis and to which material the results are to be assigned. This is necessary to enable an extrapolation of the function. From this menu the user can call the routines for the extrapolation of the set of optical constants at the short and long wavelength end, the extrapolation graphic, the selection of a special wavelength modus for the computation of the new constants and the Kramers-Kronig integration routine. When an appropriate extrapolation is found, it can also be used directly for other (of course only similar) materials.

PARMMENU is the menu to call these functions. TAKEALL selects those data necessary for the Kramers-Kronig integration from the data fields of the source material when PARMMENU is called. Furthermore, it adjusts the extrapolation parameters of the last material to the new one. This enables that different materials may be treated without calling the extrapolation routines for each of them.

6.6 Wavelength scale (WAVEMENU, WAVETAKE, WAVESCAL, ZSCAL)

Instead of the computation of the new set of optical constants at the same wavelength points as in the given set, these routines allow to chose another wavelength scaling. One can chose equidistant scales in wavelengths or wavenumbers, logarithmic scales and a system of wavelength points given by an arbitrary material loaded into the program. If the wavelength range covered by the new scale exceeds the range given by the source material, an extrapolation has to be chosen before the new wavelength scale can be initialized.

WAVEMENU is the selection menu to chose the five different treatments. WAVETAKE asks for the material which wavelength scale is to be used for the mixture and computes the set of constants at the new wavelength points from a linear interpolation between the given data points using VALUEAT or from the analytical extrapolation formulae. WAVESCAL asks the user whether to take a linear or logarithmic scale in wavelength, it asks for the starting point to be used, the stepwidth, and the number of steps (<256). For the sake of comparison it shows always the shortest and longest wavelength covered by the data set from the source material. Then, the array of wavelength points is filled with the new values and the optical constants are computed at the these data points. ZSCALE does the same except for wavenumbers instead of wavelengths.

6.7 Extrapolation routines (LEXTRA, REXTRA, LFUNCT, RFUNCT)

These routines provide the menus for the selection of an appropriate extrapolation function at the ends of the wavelength range where the optical constants of the source material are given. Until now, two general types of extrapolations are implemented: power-laws and exponential decays. As special case of a power low behaviour the Lorentz/Drude behaviour at very long and short wavelengths is predefined. This means that the exponent of the λ dependence of the real constants is 2 at short wavelengths and 0 at long wavelengths and the exponent of the λ dependence of the imaginary constants 3 at short wavelengths and -1 at long wavelengths.

To enable a convergence of the Kramers-Kronig integral the exponents have to be above 0 at short wavelengths and below 1 and 2 for the real and imaginary constants at long wavelengths. One has to be aware that chosing an non-integer exponent for the power law extrapolation results in a numerical integration outside of the wavelength range of given data which is much slower and also less accurate than the analytic integration which is possible in case of integer exponents.

Each extrapolation function may be characterized by up to 6 parameters. The existing functions use only at most 5 of them but this redundance enables a expansion of further user-defined extrapolation functions which may be simply added to the extrapolation routines. The first parameter is the number for the type of function, the other parameters may be used in the function. In the short wavelength extrapolation the real and imaginary constants do not behave symmetric because contributions from resonances at wavelengths short compared to the shortest wavelength of the given interval only affect the real constants. Independent from the set of constants used for the Kramers-Kronig integration the program asks, therefore, always for the short wavelength limit of the real constants. They are stored either in the fourth or fifth parameter. The other one is zero.

LEXTRA and REXTRA are the routines to input the extrapolation parameters. They check whether they are in the allowed ranges and compute the parameter given by the condition that the extrapolation has to meet the given data at the first and last data point. The functions LFUNCT and RFUNCT compute the values of the constants at any given wavelength in the extrapolated range.

6.8 The integration routines (KKCALC, LINT, CUBINT, RINT, SPLINE)

The Kramers-Kronig integration routines divide the total integration range for each wavelength point into three parts. The first part is a small interval around the singularity which gives no contribution to the principal value in the integral. This is fulfilled when the interval is symmetric to the singularity and the denominator function is constant within it. The interval omitted by KKCALC is determined by the condition that the function of optical data must be constant within 10^{−4} of its value. This criterion must not be weakened because the analysis may, otherwise, result in unphysical optical data (see e.g. Maron 1990). The second range is the interval where the data points are given. For each step between two wavelength points the integral is computed from the analytical formula of a cubic spline fit. The third range covers the wavelengths where the behaviour is only extrapolated. In case of power-law extrapolations with integer exponents, the contributions are computed analytically. Otherwise, they are computed numerically using the Simpson rule. As criterion for the dynamical stepwidth and the most outer interval computed, a relative accuracy of all contributions of 10^{−7} is required.

At first, KKCALC fits the behaviour of the given set of constants by cubic splines using the routine SPLINE which computes the cubic spline derivatives at all data points. SPLINE is a modified version of the SPLINE routine from the “Numerical Recipies”. Then KKCALC determines the small interval to be taken out of the integration. Then it calls CUBINT for each wavelength interval between two given data points and CUBINT computes the Kramers-Kronig integral for these intervals from the given cubic spline constants. To reduce the errors in the principal value of the integral always one interval on the left and on the right side of the wavelength point considered are treated and their contributions are sumed up. When reaching the last point of given data, the routines for the integration of the extrapolated functions LINT and RINT are called. They use the extrapolation parameters and decide whether an analytical integration is possible. In case of a numerical integration by the Simpson rule, the stepwidth is dynamically adjusted by computing the integral for each interval with two different step sizes and comparing their difference to the allowed relative deviation of 10^{−7}. Finally, KKCALC computes the other pair of optical constants from the newly determined data.

6.9 Used data structures and the data exchange between different parts

The main data fields for wavelength and the optical constants n and k and ϵ′ and ϵ′′ for all materials form the central COMMON block which is accessed by most routines. Moreover, there are COMMON blocks for all datas shown in the menus as material names, minimum and maximum wavelength, and the number of wavelength points. Their COMMON blocks are OPTCST, FNAMES, MINMAX, and ANZ. These quantities remain static during the run of the program unless the user changes them explicitiely. The same holds for the extrapolation parameters LPARM and RPARM which are not related to a special material but after their first definition available for all media which are chosen to be the result of the Kramers-Kronig analysis when entering the parameter modification menu. The data exchange concerning all other quantities is resricted to single materials or particular quantities and is performed via parameters in subroutine calls.

In the data exchange between the different parts of the programs all real data are defined as real*4. Only in the routines for the Kramers-Kronig integration all contributions to the integral are computed as real*8 to minimize the errors of the numerical integration.

7 Source files

The following files contain the sources for the different routines:

kktool.for: KKTOOL, HEADMENU, INITDATA

kkiofcts.for: READCNST, SAVECNST, READFILE, SAVEFILE, NAMECHEC

kkgraph.for: KKGRAPHIC, EXGRAPH, PLOT

graphmin.for: GRINIT, GRDRAW, GRAXIS, GRPRINT, GRCLOSE, WAIT, COORD, SMLPAUSE

kkadd.for: ADDFILES, VALUEAT

kkparmmenu.for: PARMMENU, TAKEALL

kkwave.for: WAVEMENU, WAVETAKE, WAVESCAL, ZSCALE

kkextra.for: LEXTRA, REXTRA, LFUNCT, RFUNCT

kkcalc.for: KKCALC, LINT, CUBINT, RINT, SPLINE

8 References

Maron N. 1990, ApSS 172, 21

Purcell E.M. 1969, ApJ 158, 433

Ossenkopf V., Henning Th., Mathis J.S. 1992, A&A 261, 567

Rouleau F., Martin P.G. 1991, ApJ377, 526