Users guide, Part I: Optics programs

2.a The Kramers-Kronig analysis of transmission spectra

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

Abstract

The Kramers-Kronig analysis of transmission spectra (KKTRANS) is an application of the more general Kramers-Kronig analysis routines of KKTOOL for the derivation of optical constants from transmission spectra of samples with inclusions of separated small spherical particles. The sample must be optically thin and the embedding matrix material is characterized by a constant real refractive index. A convenient extrapolation of the transmission spectrum outside of the measured range has to be done by the user. Several standard functions for the extrapolations are available and the continuation can be visually controlled. The obtained optical constants may be stored in various forms.

1 Basic algorithms

Kramers-Kronig relations do not only exist for optical constants of dielectric functions but also for various other complex functions, e.g. the composed expression f = (ϵ − 1)/(ϵ + 2). Its imaginary part is obtained from transmission spectra of diluted powder samples by

Q_{abs} = \frac{8πan_{env}}{λ} Im\frac{ϵ − 1}{ϵ + 2}  (1)

where ϵ is taken relative to the embedding medium. Then, one can obtain its real part from the Kramers-Kronig integral

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

The following technical solutions are used to carry out the Kramers-Kronig integration: The user has to give analytic functions for an extrapolation of the data outside of the measured range. 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 measured data values the imaginary part of f 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 details of the integration are identical to the according routines in KKTOOL.

Knowing real and imaginary part of f, the dielectric function of the diluted particles ϵ_{p} follows from:

ϵ_{p} = (\frac{1 + 2f}{1 − f} ) × ϵ_{env} .  (3)

The whole formalism is only valid as long as the diluted particles are considered to be spheres.

A correction for non-constant ϵ_{env} might be done later by hand. It was omitted in this program in favour of a Kramers-Kronig consistent set of dielectric functions which might be destroyed by an arbitrary dielectric function of the embedding medium.

In addition to the Kramers-Kronig integration for the selected wavelength points, a test integration for the behaviour in the static case is carried out. Here, the function f must be real and below 1 for physically reasonable materials with positive static dielectric constants (Purcell 1969). Hence, the integral of f′′/λ over the whole measured wavelength range has to be below 1.

f(0) > \frac{2}{π} ∫_{λ_{1}}^{λ_{2}} \frac{f′′(λ)}{λ}dλ .  (4)

This integral is determined parallel to the Kramers-Kronig integration for the wavelength range of measured data.

2 User interface

According to the linear structure of the program, the user has exactly to follow the dialog sequence of the program. There is only one loop where the extrapolations may be adjusted and corrected until the best visual result is obtained.

2.1 Terminal input

At first the user is asked for the name of the input file containing the extinction spectrum. Files containing Q_{abs}, Q_{abs}/a, C_{abs}/V and mass absorption coefficients are recognized. After the input of the real dielectric constant of the embedding medium the extrapolation menus are called:

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)

Here, one has to define functions for the extrapolation of the absorption spectrum outside of the range covered by the given data before the Kramers-Kronig integral can be computed. The Lorentz behaviour is a special case of the power low extrapolation in the menu points (2). In the short wavelength range it behaves as aλ^{3}. Furthermore, the contributions from wavelengths much shorter than the considered interval to the real part of the function f have to be determined. This value is computed from the real dielectric constant ϵ_{p} of the particle material at wavelengths small compared to the lowest measured wavelength and the constant of the embedding medium ϵ_{env}. In the long wavelength extrapolation the Lorentz behaviour is according to aλ^{−1}. Chosing a power law behaviour for the extrapolation the exponent must be above 0 for the short wavelength range and below 2 in the long wavelength range. The highest exponents lead, however, to an unphysical static behaviour of such materials when the particles are really spherical (see Eq. (4)). 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.

After the choice of the extrapolation functions a graphics plot showing the original and the extrapolated Im(f) data. 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. From this visual control the quality of the extrapolation may be well estimated. In case of an insufficient bahaviour of the extrapolation, both extrapolation functions may be chosen again.

When a good extrapolation is defined, the user is asked whether to use a special wavelength scale. Instead of the computation of the dielectric functions at the same wavelength points as in the given spectrum, one can chose another wavelength scaling in the wavelength modus menu:

Wavelengths modus:

     Use a linear scaling :              (1) in wavelengths (2) in wavenumbers
     Use a logarithmic scaling :         (3) in wavelengths (4) in wavenumbers

Here, one can select the best wavelength scaling for the new material. One has to 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.

After the selection of the extrapolating functions and the wavelength modus, the Kramers-Kronig integration is started. At first the test integral for the static behaviour is computed. In case of an unphysically high value, a warning is displayed that the absorption values are too high. The Kramers-Kronig analysis is, nevertheless, carried out. Finally the real and imaginary part of the computed dielectric constants are plotted graphically and a file name to store the resulting optical constants is aked for. The whole procedure may be repeated for the same or another spectrum.

2.2 Data files

All data files treated as input or output are ASCII files containing one line per wavelength point. The files with the absorption spectra consist of lines with two, the files with the optical constants of three numbers separated by blanks. The files have to consist of at least 2 lines. The contents of files are interpreted according to the file extension. A number of different file types is recognized automatically according to their 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 following columns. For absorption spectra “qa” stands for Q_{abs} (the radius a has to be given by hand), “qr” means Q_{abs}/a (in μm^{−1}), “cv” stands for C_{abs}/V (in cm^{−1}), and “kp” means mass absorption coefficients C_{abs}/(V ρ) (in cm^{2}/g). For files of optical constants the extension letters “nk” state that the second and third column 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 160 and 180 kB for the PC versions and 650 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 kktrans to enable a correct handling of the graphics window and to start the Kramers-Kronig analysis of transmission spectra 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. The files containing the measured spectra 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 (KKTRANS)

This routine is the most outer surface of the program. It calls the file input routines to read the absorption spectrum, the routines for the definition of the extrapolation and the wavelength points. Then it calls the computational routines, the graphic routine to visualize the resulting optical constants and the file output routines to store the results.

6.2 File operations (READCNST, READTABLE, NAMEXT, SAVECNST, SAVEFILE, NAMECHEC)

The file input reads ASCII files with two numbers per line. separated by at least one blank. The contents of a file are interpreted according to the file extension. A number of different file types is recognized automatically. When the first letter of the file extension is an “l”, the first column of the file contains the wavelengths (in cm). A letter “z” means wavenumbers (in cm^{−1}) in the first column. The next two letters of the file extension stand for the content of the second column: qa - Q_{abs} qr - Q_{abs}/a (in μm^{−1}) cv - C_{abs}/V (in cm^{−1}) , kp - C_{abs}/(V ρ) (in cm^{2}/g).

The output produces files with lines containing three ASCII numbers separated by blanks. The extension determines which kind of data shall be written to the file. The first letter again characterizes the content of the first column. The following two letters detemine what is written to the next two columns: nk - n k ee - ϵ′ ϵ′′. Both for file input and output all files with other extensions are rejected.

In principle the routines for file input and output work in the same way. An interface routine (READCNST and SAVECNST) adapts the array of data read or to be written to the data structure used by the other parts of the program. This includes the selection of the appropriate data arrays and the organization of the data in an ascending order with wavelength. The real file operations are done by READTABLE and SAVEFILE. They ask for a filename and transform the data depending on the result of NAMECHEC or NAMEEXT in the required format. The data exchange between READTABLE and READCNST only contains a field of wavelengths and of Q_{abs}/a (in cm^{−1}), the exchange between SAVEFILE and SAVECNST only wavelength n and k. The translation into these quantities is done within READTABLE and SAVEFILE. They are the procedures reading and writing the data on the disks. NAMECHEC and NAMEEXT analyse the filename given by the user, remove the path from it and extracts its extension. The extension is translated into an integer characterizing the content of the file. NAMECHEC recognizes the extensions for files of optical constants, NAMEEXT those for files of absorption spectra.

6.3 Graphics (GRAPH, EXGRAPH , PLOT, GRAPHMIN)

The graphics routines allow to show the original absorption spectrum together with the functions chosen for extrapolation and the resulting optical constants of the particle material n and k. 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.

EXGRAPH is the subroutine for the demonstration of the extrapolation. It is called by PARMMENU after the definition of the extrapolation. It uses the measured transmission data 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 which shows the behaviour of th extrapolation graphically. Only the quantity f′′ is plotted versus wavelength. The routine GRAPH is called when the Kramers-Kronig integration is finished. It selects the field of refractive indices n and absorption coefficients k and calls the subroutine PLOT twice to show these data. 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 Extrapolation routines (PARMMENU, LEXTRA, REXTRA, LFUNCT, RFUNCT)

Before the Kramers-Kronig integration is carried out, continuations of the function outside of the given range have to be defined. The extrapolation routines provide the menus for the selection of an appropriate extrapolation function at the ends of the wavelength range where the absorption spectrum was measure. 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 f′′ is 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 2 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, one has to define additionally the short wavelength limit of the real dielectric function of the particle material. The according f′ value is stored in the fifth parameter.

PARMMENU is the shell which calls the routines for the definition of the extrapolations. It then calls EXGRAPH to enable a visual control of the continuation. When the user is content with the extrapolation, the wavelength scale menu WAVEMENU and the Kramers-Kronig integration CALC are called.

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.5 Wavelength scale (WAVEMENU, WAVESCAL, ZSCAL, VALUEAT)

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, linear or logarithmic scales. If the wavelength range covered by the new scale exceeds the range given by the source material the extrapolated functions ser used outside.

WAVEMENU is the selection menu to chose the four different treatments. 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. The data values at the new wavelength points are computed from a linear interpolation between the given data points using VALUEAT or from the analytical extrapolation formulae. Finally, 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. The procedure VALUEAT interpolates the optical constants of one materials at the selected wavelength points by a linear function.

6.6 The integration routines (KKCALC, INFINT, 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, it calls INFINT which computes the Kramers-Kronig integral relative to 0 wavelenths using the cubic spline parameters. The integration includes only the wavelength range of measured data. When the integral value exceeds π/2 a warning message is displayed that the measured absorption values are unphysically high. Now, the Kramers-Kronig integral is computed for each wavelength point. KKCALC determines the small interval around the wavelength point to be taken out of the integration. 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.7 Used data structures and the data exchange between different parts

According to the data structure of KKTOOL, there is on main COMMON block OPTCST containing the optical data. In the first field the wavelength points are stored, the next two fields contain the resulting optical constants and the last field contains the absorption spectrum. All routines which use these optical data in their standard form include this COMMON block. This does not hold for routines as READTABLE, SAVEFILE or PLOT which use only single fields of optical data. The COMMON block ANZ contains the number of points defined in each field of OPTCST. The COMMON block KBR contains the refractive index of the embedding medium. These quantities remain static during the run of the program unless the user changes them explicitiely. 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:

kktrans.for: KKTRANS, READCNST, READTABLE, NAMEEXT, SAVECNST, GRAPH, EX-GRAPH, KKCALC, INFINT

kkdefine.for: PARMMENU, LEXTRA, REXTRA, WAVEMENU, WAVESCAL, ZSCAL kkcloned.for: SAVEFILE, NAMECHEC, PLOT, LFUNCT, RFUNCT, VALUEAT, CUBINT, LINT,

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

All routines within kkcloned.for are completely identical to the routines with the same names in KKTOOL.

8 References

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