High rate line-by-line calculation method and program by forked line type resolution

Abstract

The invention offers a calculation method and program capable of performing line-by-line calculations using a Voigt function at speeds of 50-100 times what is conventional. The Voigt function is divided into a first range around the peak and a skirt portion not contained in the first range. The first range is replaced by a cubic function, and the skirt portion is taken as the Voigt function to perform calculations in predetermined ranges of equal intervals. Furthermore, the peak area of the first range is replaced by a cubic function, and the skirt portion is taken as a function representing the difference between the Voigt function and the cubic function to perform calculations in second predetermined intervals smaller than the aforementioned first predetermined intervals. This is repeated until the desired level of precision is reached. Additionally, interpolation is performed by dividing these predetermined intervals into four or five parts.

Description

TECHNICAL FIELD

The present invention relates to high-speed line-by-line calculation of the absorption coefficient spectra of gas molecules necessary for atmospheric observations.

BACKGROUND ART

While the Earth's atmosphere receives electromagnetic radiation from the sun, it is itself releasing electromagnetic radiation mainly in the infrared region into space. This electromagnetic radiation interacts such as by absorption, radiation and scattering with gas molecules, clouds, aerosols and the like present in the atmosphere. By measuring changes in the spectral distribution of electromagnetic radiation caused by such interactions with the atmosphere (theoretical calculations of such distributions are referred to as radiative transfer calculations) using sensors on land or in airplanes or satellites, it is possible to study the air temperature, air pressure and concentrations of atmospheric trace gas constituents.

While a technique known as line-by-line calculation is used for radiative transfer calculations of the atmosphere, line-by-line calculations involve a enormous number of calculations. For this reason, the calculation time becomes extremely long, making it practically difficult to perform calculations over a wide spectral range at high speeds using personal computers or the like. Additionally, there are methods that make use of data obtained in line-by-line calculations beforehand involving table lookup routines as methods for avoiding the problem discussed above, but such methods have the drawback of requiring large data files. Additionally, other algorithms not requiring prior calculation have fared no better than about 20 times the speed of line-by-line calculations.

Non-Patent Publication 1: Akihiro Uchiyama, “Line-by-Line Computation of the Atmospheric Absorption Spectrum Using the Decomposed VOIGT Line Shape”, J. Quant. Spectrosc. Radiat. Transfer, Vol. 47, No. 6, pp. 521-532, 1992.

DISCLOSURE OF THE INVENTION

Problems to be Solved by the Invention

Due to the problems mentioned above, calculation methods capable of line-by-line calculations at high speed have been sought, and the present invention offers such a calculation method and program.

Means for Solving the Problems

The present invention relations to a Voigt function approximation calculating program used for line-by-line calculations, said program performing:

• • (1) a step of dividing a domain of a Voigt function into a first range around the peak of the Voigt function and a skirt portion not contained in the first range, replacing the first range with a cubic function, calculating the values and derivatives of said cubic function and the Voigt function in the skirt portion for each of first predetermined intervals, and connecting said cubic function and said Voigt function at points of connection thereof using the values and derivatives of both functions;
  • (2) a step of adding together the results of step (1) for a plurality of absorption lines;
  • (3) a step of calculating the function values and derivatives for the results of step (2) by interpolation over intervals smaller than said first predetermined intervals;
  • (4) a step of dividing said first range into a second range near the peak and a skirt portion not contained in the second range, replacing said second range of a “function representing the difference between the Voigt function and said cubic function” with a cubic function, and calculating values and derivatives of said cubic function and said “function representing the difference between the Voigt function and said cubic function” in the skirt portion for each of second predetermined intervals;
  • (5) a step of connecting said cubic function and said “function representing the difference between the Voigt function and said cubic function” at points of connection thereof using the values and derivatives of both functions;
  • (6) a step of adding the results of steps (4) and (5) to the results of step (3) for a plurality of absorption lines;
  • (7) a step of calculating the function values and derivatives for the results of step (6) by interpolation over intervals smaller than said second predetermined intervals; and
  • (8) a step of adding the values of the “function representing the difference between the Voigt function and said cubic function” to the results of step (7) for a plurality of absorption lines in said second range.

A method using the above steps allows a Voigt function to be calculated at high speed using a cubic function.

In addition thereto, the steps (5) through (8) may be repeated one or more times until a third predetermined interval is reached. As a result, it is possible to obtain detailed results with small calculation intervals. The first through third predetermined intervals can be determined as described below.

The first predetermined interval for the widest sub-function is j^kmadv. Here, j is a single-digit natural number, dv is the increment in wave number, and kmax is the largest natural number satisfying the relationship j^kma+2pdv≦Vmax. However, Vmax represents the maximum calculation range from the center of the absorption line, and p is a natural number (p=1, 2, 3) for controlling the calculation precision

The most specific third predetermined interval is j^kmimdv. Here, j is a single-digit natural number of about 2-6, dv is the increment of the wave number, and kmin is the maximum non-negative decimal fraction (0 when non-existent) satisfying the relationship j^kminpdv≦α (a being approximately γ/4). However, γ is an approximate value of the full-width at half-max of the absorption line, and p is a natural number (p=1, 2, 3) for controlling the calculation precision.

A method in accordance with any one of claims 1-4, wherein the second predetermined intervals are determined using the following equation.

For the sub-function with the (k−kmin+1)-th smallest width, the predetermined interval is j^kdv. Here, j is a single-digit natural number, dv is the increase in wave number and k is such that kmin≦k<kmax.

For a quartered interpolation, the interpolation is calculated with j set to 4, using the function values y₀, y₁ and function derivative values y₀′, y₁′ at x₀, X₁ in the interpolation interval (x₀, x₁), with the below-given Equation (1) as a function value interpolation equation, Equation (2) as a function derivative value interpolation equation and ε as a non-negative decimal fraction. $\begin{array}{cc}\left[\mathrm{Eq}.1\right]& \\ \left(\begin{array}{c}{y}_a\\ y_b\\ y_c\end{array}\right)=\frac{1}{64}\left(\begin{array}{cccc}54-6\varepsilon & 10+6\varepsilon & 9\left(1-\varepsilon \right)& -3\left(1-\varepsilon \right)\\ 32& 32& 8\left(1-\varepsilon \right)& -8\left(1-\varepsilon \right)\\ 10+6\varepsilon & 54-6\varepsilon & 3\left(1-\varepsilon \right)& -9\left(1-\varepsilon \right)\end{array}\right)\left(\begin{array}{c}{y}_0\\ y_1\\ \left(x_1-x_0\right)y_0^{\prime }\\ \left(x_1-x_0\right)y_1^{\prime }\end{array}\right)& \left(1\right)\\ \left[\mathrm{Eq}.2\right]& \\ \left(\begin{array}{c}{y}_a^{\prime }\\ y_b^{\prime }\\ y_c^{\prime }\end{array}\right)=\frac{1}{16\left(x_1-x_0\right)}\left(\begin{array}{cccc}-18+2\varepsilon & 18-2\varepsilon & 3\left(1-\varepsilon \right)& -5\left(1-\varepsilon \right)\\ -24+8\varepsilon & 24-8\varepsilon & -4\left(1-\varepsilon \right)& -4\left(1-\varepsilon \right)\\ -18+2\varepsilon & 18-2\varepsilon & -5\left(1-\varepsilon \right)& 3\left(1-\varepsilon \right)\end{array}\right)\left(\begin{array}{c}{y}_0\\ y_1\\ \left(x_1-x_0\right)y_0^{\prime }\\ \left(x_1-x_0\right)y_1^{\prime }\end{array}\right)& \left(2\right)\end{array}$

For an interpolation divided into five parts, the interpolation is calculated with j set to 5, using the function values y₀, y₁ and function derivative values y₀′, y₁′ at x₀, x₁ in the interpolation interval (x₀, x₁), with the below-given Equation (3) as a function value interpolation equation, Equation (4) as a function derivative value interpolation equation and E as a non-negative decimal fraction. $\begin{array}{cc}\left[\mathrm{Eq}.3\right]& \\ \left(\begin{array}{c}{y}_a\\ y_b\\ y_c\\ y_d\end{array}\right)=\frac{1}{125}\left(\begin{array}{cccc}112-12\varepsilon & 13+12\varepsilon & 16\left(1-\varepsilon \right)& -4\left(1-\varepsilon \right)\\ 81-6\varepsilon & 44+6\varepsilon & 18\left(1-\varepsilon \right)& -12\left(1-\varepsilon \right)\\ 44+6\varepsilon & 81-6\varepsilon & 12\left(1-\varepsilon \right)& -18\left(1-\varepsilon \right)\\ 13+12\varepsilon & 112-12\varepsilon & 4\left(1-\varepsilon \right)& -16\left(1-\varepsilon \right)\end{array}\right)\left(\begin{array}{c}{y}_0\\ y_1\\ \left(x_1-x_0\right)y_0^{\prime }\\ \left(x_1-x_0\right)y_1^{\prime }\end{array}\right)& \left(3\right)\\ \left[\mathrm{Eq}.4\right]& \\ \left(\begin{array}{c}{y}_a^{\prime }\\ y_b^{\prime }\\ y_c^{\prime }\\ y_d^{\prime }\end{array}\right)=\frac{1}{25\left(x_1-x_0\right)}\left(\begin{array}{cccc}-24-\varepsilon & 24+\varepsilon & 8\left(1-\varepsilon \right)& -7\left(1-\varepsilon \right)\\ -36+11\varepsilon & 36-11\varepsilon & -3\left(1-\varepsilon \right)& -8\left(1-\varepsilon \right)\\ -36+11\varepsilon & 36-11\varepsilon & -8\left(1-\varepsilon \right)& -3\left(1-\varepsilon \right)\\ -24-\varepsilon & 24+\varepsilon & -7\left(1-\varepsilon \right)& 8\left(1-\varepsilon \right)\end{array}\right)\left(\begin{array}{c}{y}_0\\ y_1\\ \left(x_1-x_0\right)y_0^{\prime }\\ \left(x_1-x_0\right)y_1^{\prime }\end{array}\right)& \left(4\right)\end{array}$

Additionally, the invention offers a method for achieving high speeds, when assuming the Voigt function to be K(x, y) and the difference from the Voigt profile of the absorption line to be K(x, y)+f(x), by replacing $\begin{array}{cc}\tilde{K}\left(x,y\right)=AK\left(x,y\right)+\mathrm{Bf}\left(x\right)\mathrm{and}& \left[\mathrm{Eq}.5\right]\\ \frac{\partial K\left(x,y\right)}{\partial x}& \left[\mathrm{Eq}.6\right]\end{array}$

Furthermore, the invention offers a method for achieving high speeds, when assuming the difference from the Voigt profile to be K(x, y)f(x), by replacing $\begin{array}{cc}\tilde{K}\left(x,y\right)=K\left(x,y\right)f\left(x\right)\mathrm{and}& \left[\mathrm{Eq}.7\right]\\ \frac{\partial \tilde{K}\left(x,y\right)}{\partial x}=\frac{\partial K\left(x,y\right)}{\partial x}f\left(x\right)+K\left(x,y\right)\frac{\partial f\left(x\right)}{\partial x}& \left[\mathrm{Eq}.8\right]\end{array}$

Additionally, in line-mixing correction: $\begin{array}{cc}\tilde{K}\left(x,y\right)=AK\left(x,y\right)+\mathrm{BL}\left(x,y\right)\mathrm{and}& \left[\mathrm{Eq}.9\right]\\ \frac{\partial \tilde{K}\left(x,y\right)}{\partial x}=-2\left[\left(\mathrm{Ax}+\mathrm{By}\right)K\left(x,y\right)-\left(\mathrm{Ay}-\mathrm{Bx}\right)L\left(x,y\right)-\frac{B}{\sqrt{\pi }}\right]& \left[\mathrm{Eq}.10\right]\end{array}$
Here, L(x, y) is the imaginary component of the function w(z) (the real part is the Voigt function), where complex number z=x+iy, defined by the following equation. $\begin{array}{cc}w\left(z\right)=\frac{i}{\pi }{\int }_{-\infty }^{\infty }\frac{\exp \left(-t^2\right)}{z-t}dt=\exp \left(-z^2\right)\mathrm{erfc}\left(-iz\right)=K\left(x,y\right)+iL\left(x,y\right)& \left[\mathrm{Eq}.11\right]\end{array}$
(erfc(z) is a complex complementary error function).

Furthermore, the invention offers a device for performing the method of the present invention by loading the method described above as a program.

BRIEF DESCRIPTION OF THE DRAWINGS

[FIG. 1] FIG. 1 is a schematic view of a Voigt function.

[FIG. 2] FIG. 2 is a diagram showing how a Voigt function is divided into sub-functions according to the method of the present invention.

[FIG. 3] FIG. 3 is a diagram showing the error between a Voigt function and a cubic function.

EXPLANATION OF REFERENCE NUMBERS

• • 1 Voigt function
  • 2 first range
  • 3 skirt portion
  • 4 function contained in first range
  • 5 portion not contained in first range

BEST MODE FOR CARRYING OUT THE INVENTION

An embodiment of the present invention shall be described below. FIG. 1 is a diagram showing the shape of a Voigt function 1. According to the present invention, this Voigt function 1 is first divided into sub-functions. These sub-functions are in a first range 2 around the peak and a skirt portion 3 outside the first range 2. This first range 1 is replaced by a cubic function, and the portion not contained in the first range, in other words, the skirt portion, is calculated with the Voigt function.

The details of the calculation method are shown in FIG. 2. As shown in FIG. 2(a), the function 4, which is the portion of the Voigt function 1 contained in the first range 2, is replaced by a cubic function. At this time, the cubic function is determined so as to smooth the connections by matching the derivatives of the functions at the points of connection between the cubic function and the Voigt function. Additionally, the portion 5 not contained in the first range is calculated by the Voigt function. The calculation is performed over a first predetermined interval, indicated by j^kmaxdv where dv is the change in wave number. Here, j is a single-digit natural number, and kmax is the largest natural number satisfying the relationship j^kmax+2pdv≦Vmax. Vmax may, for example, be 25 cm⁻¹. Additionally, j is normally a number of about 2-6, and is preferably 4. Also, Vmax represents the calculation range from the center of the absorption line, and p is a natural number for controlling the calculation precision, p preferably having a value of 1, 2 or 3.

Results calculated in this way are added together for a plurality of absorption lines. At this stage, the result is an approximate value calculated in the first predetermined interval.

Next, unevaluated points are calculated over an interval smaller than the first predetermined interval. That is, the unevaluated points are interpolated over intervals obtained by quartering the first predetermined interval. The calculation of the point of division is performed by the following interpolation formula.

The function values y₀, y₁ and function derivative values y₀′, y₁′ at x₀, x₁ in the interpolation interval (x₀, x₁) can be determined by the following two equations: $\begin{array}{cc}\left[\mathrm{Eq}.12\right]& \\ \left(\begin{array}{c}{y}_a\\ y_b\\ y_c\end{array}\right)=\frac{1}{64}\left(\begin{array}{cccc}54-6\varepsilon & 10+6\varepsilon & 9\left(1-\varepsilon \right)& -3\left(1-\varepsilon \right)\\ 32& 32& 8\left(1-\varepsilon \right)& -8\left(1-\varepsilon \right)\\ 10+6\varepsilon & 54-6\varepsilon & 3\left(1-\varepsilon \right)& -9\left(1-\varepsilon \right)\end{array}\right)\left(\begin{array}{c}{y}_0\\ y_1\\ \left(x_1-x_0\right)y_0^{\prime }\\ \left(x_1-x_0\right)y_1^{\prime }\end{array}\right)& \left(1\right)\\ \left[\mathrm{Eq}.13\right]& \\ \left(\begin{array}{c}{y}_a^{\prime }\\ y_b^{\prime }\\ y_c^{\prime }\end{array}\right)=\frac{1}{16\left(x_1-x_0\right)}\left(\begin{array}{cccc}-18+2\varepsilon & 18-2\varepsilon & 3\left(1-\varepsilon \right)& -5\left(1-\varepsilon \right)\\ -24+8\varepsilon & 24-8\varepsilon & -4\left(1-\varepsilon \right)& -4\left(1-\varepsilon \right)\\ -18+2\varepsilon & 18-2\varepsilon & -5\left(1-\varepsilon \right)& 3\left(1-\varepsilon \right)\end{array}\right)\left(\begin{array}{c}{y}_0\\ y_1\\ \left(x_1-x_0\right)y_0^{\prime }\\ \left(x_1-x_0\right)y_1^{\prime }\end{array}\right)& \left(2\right)\end{array}$
Here, ε is a non-negative decimal fraction. With regard to ε, even if the Voigt function and cubic function match at the interpolation points x₀, x₁ as shown in FIG. 3, there will be an error in the function values due to differences in the curvature of the curves in the region between the interpolation points. While adequate precision can be obtained even if E=0, the error can be further reduced by about half by changing the value of ε to adjust the curvature of the cubic function. From these two equations it is possible to determine the function values y_a, y_b, y_c and the derivative values y_a′, y_b′, y_c′ at points where the interpolation interval (x₀, x₁) of the aforementioned first predetermined interval has been quartered. Here, x₀, x₁ are wave numbers in the present invention. The aforementioned adjustment of ε is performed by comparing the calculation results when changing ε with the calculation results for a normal line-by-line method for a small spectral range which is only a part of the spectral range to be calculated.

Additionally, the present invention proposes an interpolation in which the first predetermined interval is divided into five parts. The interpolation equations for dividing into five parts are as follows: $\begin{array}{cc}\left[\mathrm{Eq}.14\right]& \\ \left(\begin{array}{c}{y}_a\\ y_b\\ y_c\\ y_d\end{array}\right)=\frac{1}{125}\left(\begin{array}{cccc}112-12\varepsilon & 13+12\varepsilon & 16\left(1-\varepsilon \right)& -4\left(1-\varepsilon \right)\\ 81-6\varepsilon & 44+6\varepsilon & 18\left(1-\varepsilon \right)& -12\left(1-\varepsilon \right)\\ 44+6\varepsilon & 81-6\varepsilon & 12\left(1-\varepsilon \right)& -18\left(1-\varepsilon \right)\\ 13+12\varepsilon & 112-12\varepsilon & 4\left(1-\varepsilon \right)& -16\left(1-\varepsilon \right)\end{array}\right)\left(\begin{array}{c}{y}_0\\ y_1\\ \left(x_1-x_0\right)y_0^{\prime }\\ \left(x_1-x_0\right)y_1^{\prime }\end{array}\right)& \left(3\right)\\ \left[\mathrm{Eq}.15\right]& \\ \left(\begin{array}{c}{y}_a^{\prime }\\ y_b^{\prime }\\ y_c^{\prime }\\ y_d^{\prime }\end{array}\right)=\frac{1}{25\left(x_1-x_0\right)}\left(\begin{array}{cccc}-24-\varepsilon & 24+\varepsilon & 8\left(1-\varepsilon \right)& -7\left(1-\varepsilon \right)\\ -36+11\varepsilon & 36-11\varepsilon & -3\left(1-\varepsilon \right)& -8\left(1-\varepsilon \right)\\ -36+11\varepsilon & 36-11\varepsilon & -8\left(1-\varepsilon \right)& -3\left(1-\varepsilon \right)\\ -24-\varepsilon & 24+\varepsilon & -7\left(1-\varepsilon \right)& 8\left(1-\varepsilon \right)\end{array}\right)\left(\begin{array}{c}{y}_0\\ y_1\\ \left(x_1-x_0\right)y_0^{\prime }\\ \left(x_1-x_0\right)y_1^{\prime }\end{array}\right)& \left(4\right)\end{array}$
The symbols used in the equations are the same as in the case of quartering described above. Additionally, those skilled in the art would recognize that y_d, y_d′ clearly represent a function value and a derivative value. Here, the interpolation equation described above can be used not only for the first predetermined interval, but also for other predetermined intervals, and this should also be obvious to those skilled in the art.

The interpolation points are calculated as described above. Next, as shown in FIG. 2(b), further calculations are performed by interpolation of the region around the peak in the first range. As described above, the area around the peak in the first range is called the second range 6, and divided from a skirt portion not contained in the second range 6. Then the second range is represented by a second cubic function, and the skirt portion not contained in the second range is represented by a “function representing the difference between the Voigt function and the first cubic function” 7. This second cubic function represents the difference between the new cubic function and the first cubic function. Here, this “function representing the difference between the Voigt function and the first cubic function” shall henceforth be referred to simply as the “difference function”. As described above, the second cubic function and the “difference function” 7 are again joined under the condition that at the points of connection, the function values and derivative values are matched.

After connecting the functions in this way, the function value and derivative value are calculated in the second predetermined interval for the second range. This second predetermined interval is smaller than the first predetermined interval and is calculated by the equation given below.

For the sub-function with the (k−kmin+1)-th smallest width, the predetermined interval is j^kdv. Here, j is a single-digit natural number, dv is the increase in wave number and k is such that kmin≦k<kmax.

Additionally, these calculations are performed for all absorption lines. Furthermore, interpolation is performed between the calculated points. The equations (1) and (2) or (3) and (4) are used for the interpolation.

Here, the above method for calculating a “difference function” for a skirt portion not contained in the second range and a cubic function for the second range that connects smoothly therewith for a plurality of absorption lines has been explained, but it is also possible to use a method of calculating the Voigt function of the skirt portion not contained in the second range and a cubic function connecting smoothly therewith for a plurality of absorption lines, then subtracting the first cubic function from the entire first range for a plurality of absorption lines.

With regard to this second range, the above calculations may be repeated until a predetermined interval is reached. That is, the second range may be further divided into a peak area and a skirt portion which are respectively calculated by a cubic function and difference function, then connected, and the calculation performed with k in the equation for calculating the second predetermined interval incremented by 1. This calculation is performed for all absorption lines, and these are added to the results calculated in FIG. 2(b). The steps described above are repeated until the predetermined interval is reached.

Further, a calculation is performed after replacing the peak area with a cubic function and making the skirt portion a “difference function”, over a third predetermined interval smaller than the second interval. This third predetermined interval is represented by the equation given below. The third predetermined interval is j^kmindv. Here, j is a single-digit natural number, dv is the increment of the wave number, and kmin is the maximum non-negative decimal fraction satisfying the relationship j^kminpdv≦α(if non-existent, it is set to 0). However, a has a value of about γ/4 (γ being an estimate of the full width at half max of the absorption line), and p is a natural number (p=1, 2, 3) controlling the calculation precision.

By following the above steps, it is possible to calculate a Voigt function precisely and at high speed. Here, it is known that there is a slight difference from the Voigt function of carbon dioxide and the like, and a method for correcting for this difference is known. The method of the present invention can also be applied to such corrections. These corrections shall be explained next.

General profile correction shall be explained first.

The Voigt function is represented by K(x, y). When the difference from the Voigt profile is represented by the shape of K(x, y)+f(x), the high-speed technique can be applied by setting:

[Eq. 16]
{tilde over (K)}(x,y)
as the function to be substituted, and replacing it with: $\begin{array}{cc}\tilde{K}\left(x,y\right)=AK\left(x,y\right)+\mathrm{Bf}\left(x\right)& \left[\mathrm{Eq}.17\right]\\ \frac{\partial \tilde{K}\left(x,y\right)}{\partial x}=A\frac{\partial K\left(x,y\right)}{\partial x}+B\frac{\partial f\left(x\right)}{\partial x}& \left[\mathrm{Eq}.18\right]\end{array}$
Similarly, if the difference from the Voigt profile is represented by the shape of K(x, y)f(x), then the high-speed technique can be performed by replacing it with $\begin{array}{cc}\tilde{K}\left(x,y\right)=K\left(x,y\right)f\left(x\right)& \left[\mathrm{Eq}.19\right]\\ \frac{\partial \tilde{K}\left(x,y\right)}{\partial x}=\frac{\partial K\left(x,y\right)}{\partial x}f\left(x\right)+K\left(x,y\right)\frac{\partial f\left(x\right)}{\partial x}& \left[\mathrm{Eq}.20\right]\end{array}$

Explaining the above correction method with a specific example, in a sub-Lorentzian correction: $\begin{array}{cc}K\left(x,y\right)& \left[\mathrm{Eq}.21\right]\\ \frac{\partial K\left(x,y\right)}{\partial x}& \left[\mathrm{Eq}.22\right]\end{array}$
will respectively become: $\begin{array}{cc}\tilde{K}\left(x,y\right)=K\left(x,y\right)A\exp \left(-Bx\right)& \left[\mathrm{Eq}.23\right]\\ \frac{\partial \tilde{K}\left(x,y\right)}{\partial x}=\frac{\partial K\left(x,y\right)}{\partial x}A\exp \left(-Bx\right)+K\left(x,y\right)\left[-\mathrm{sgn}\left(x\right)\mathrm{AB}\exp \left(-Bx\right)\right]& \left[\mathrm{Eq}.24\right]\end{array}$
Here, A and B are correction coefficients. Additionally, those skilled in the art will recognize that sgn(x) is the sign function.

For line-mixing correction, the following relationships arise: $\begin{array}{cc}\frac{\partial K\left(x,y\right)}{\partial x}=2\left[yL\left(x,y\right)-xK\left(x,y\right)\right]& \left[\mathrm{Eq}.25\right]\\ \frac{\partial L\left(x,y\right)}{\partial x}=2\left[xL\left(x,y\right)+yK\left(x,y\right)-\frac{1}{\sqrt{\pi }}\right]\mathrm{Therefore},& \left[\mathrm{Eq}.26\right]\\ K\left(x,y\right)& \left[\mathrm{Eq}.27\right]\\ \frac{\partial K\left(x,y\right)}{\partial x}& \left[\mathrm{Eq}.28\right]\end{array}$
will respectively be $\begin{array}{cc}\tilde{K}\left(x,y\right)=\mathrm{AK}\left(x,y\right)+\mathrm{BL}\left(x,y\right)& \left[\mathrm{Eq}.29\right]\\ \frac{\partial \tilde{K}\left(x,y\right)}{\partial x}=-2\left[\left(\mathrm{Ax}+\mathrm{By}\right)K\left(x,y\right)-\left(\mathrm{Ay}-\mathrm{Bx}\right)L\left(x,y\right)-\frac{B}{\sqrt{\pi }}\right]& \left[\mathrm{Eq}.30\right]\end{array}$
Here, L(x, y) is the imaginary component of the function w(z), where complex number z=x+iy, defined by the following equation. The real part is the Voigt function. $\begin{array}{cc}w\left(z\right)=\frac{i}{\pi }{\int }_{-\infty }^{\infty }\frac{\exp \left(-t^2\right)}{z-t}dt=\exp \left(-z^2\right)\mathrm{erfc}\left(-iz\right)=K\left(x,y\right)+iL\left(x,y\right)& \left[\mathrm{Eq}.31\right]\end{array}$
erfc(z) is a complex complementary error function.

According to the calculation methods of the present invention described above, it is possible to perform line-by-line calculations using Voigt functions with good precision and at high speed. By further loading these methods as computer programs, it is possible to result in specific devices. Such a program may, for example, be a stand-alone type program operating on a personal computer. It may also be provided in the form of a plug-in for other software.

Claims

1. A method for increasing the speed of line-by-line calculations of multiple overlapping absorption lines, said method comprising:

(1) a step of dividing a domain of a Voigt function representing the shape of an absorption line into a first range around the peak of the Voigt function and a skirt portion not contained in the first range, replacing the first range with a cubic function whose function values and derivative values match with those of said Voigt function at the connection points with said Voigt function, and calculating the values and derivative values of said cubic function and the Voigt function in the skirt portion for each of first predetermined intervals;

(2) a step of adding together the results of step (1) for a plurality of absorption lines;

(3) a step of interpolating the results of step (2) to calculate the function values and derivative values in second predetermined intervals smaller than said first predetermined intervals;

(4) a step of dividing said first range into a second range near the peak and a skirt portion not contained in the second range, replacing said second range of a “function representing the difference between the Voigt function and said cubic function” with a second cubic function whose function values and derivative values match with those of said “function representing the difference between the Voigt function and said cubic function” at the connection points with said “function representing the difference between the Voigt function and said cubic function”, and calculating the values and derivative values of said second cubic function replacing said second range and the “function representing the difference between the Voigt function and said cubic function” in the skirt portion for said second predetermined intervals;

(5) a step of adding the results of step (4) to the results of step (3) for a plurality of absorption lines;

(6) a step of interpolating the results of step (5) to calculate the function values and derivative values for intervals smaller than said second predetermined intervals; and

(7) a step of adding the values and derivative values of the “function representing the difference between the ““function representing the difference between the Voigt function and said cubic function”” and said second cubic function” to the results of step (6) for a plurality of absorption lines in said second range.

2. A method in accordance with claim 1, wherein (8) the step of calculating function values and derivative values by interpolation is repeated while narrowing the interval until it becomes an interval of a minimum unit.

3. A method in accordance with claim 1, wherein said steps (4) through (6) are repeated one or more times until the “intervals smaller than said second predetermined intervals” in step (6) reach third predetermined intervals.

4. A method in accordance with claim 1, wherein said predetermined interval is determined by using the following equation.

Said first predetermined intervals are jkmaxdv. Here, j is a single-digit natural number, dv is the increment in wave number, and kmax is the largest natural number satisfying the relationship jkmax+2pdv≦Vmax, where Vmax represents the maximum calculation range from the center of the absorption line, and p is a natural number for controlling the calculation precision.

5. A method in accordance with claim 3, wherein said predetermined interval is determined using the following equation.

The third predetermined intervals are jkmindv. Here, j is a single-digit natural number, dv is the increment of the wave number, and kmin is the maximum non-negative decimal number satisfying the relationship jkminpdv≦α (α being approximately γ/4) where γ is an approximate value of the full-width at half-maximum of the absorption line, and p is a natural number for controlling the calculation precision.

6. A method in accordance with claim 1, wherein said second predetermined interval is determined by using the following equation.

For the sub-function with the (k−kmin+1)-th smallest width, the predetermined interval is jkdv. Here, j is a single-digit natural number, dv is the increase in wave number and k is such that kmin≦k<kmax. kmax is a natural number satisfying the relationship jkmax+2pdv=Vmax, where Vmax represents the maximum calculation range from the center of the absorption line, and p is a natural number for controlling the calculation precision. kmin is the maximum non-negative decimal number satisfying the relationship jkminpdv=α (α being approximately γ/4), where γ is an approximate value of the full-width at half-maximum of the absorption line, and p is a natural number for controlling the calculation precision.

7. A method in accordance with any one of claims 1-3, wherein said interpolation is calculated by quartering the interpolation interval, with the function values represented by ya, yb, yc and the derivative values represented by ya′, yb′, yc′ at the points where the interpolation interval (x0, x1) is divided, using the function values y0, y1 and function derivative values y0′, y1′ at x0, x1, with the below-given Equation (1) as a function value interpolation equation, Equation (2) as a function derivative value interpolation equation and ε as a non-negative decimal fraction. ( y a y b y c ) = 1 64 ⁢ ( 54 - 6 ⁢ ɛ 10 + 6 ⁢ ɛ 9 ⁢ ( 1 - ɛ ) - 3 ⁢ ( 1 - ɛ ) 32 32 8 ⁢ ( 1 - ɛ ) - 8 ⁢ ( 1 - ɛ ) 10 + 6 ⁢ ɛ 54 - 6 ⁢ ɛ 3 ⁢ ( 1 - ɛ ) - 9 ⁢ ( 1 - ɛ ) ) ⁢ ( y 0 y 1 ( x 1 - x 0 ) ⁢ y 0 ′ ( x 1 - x 0 ) ⁢ y 1 ′ ) ( 1 ) ( y a ′ y b ′ y c ′ ) = 1 16 ⁢ ( x 1 - x 0 ) ⁢ ( - 18 + 2 ⁢ ɛ 18 - 2 ⁢ ɛ 3 ⁢ ( 1 - ɛ ) - 5 ⁢ ( 1 - ɛ ) - 24 + 8 ⁢ ɛ 24 - 8 ⁢ ɛ - 4 ⁢ ( 1 - ɛ ) - 4 ⁢ ( 1 - ɛ ) - 18 + 2 ⁢ ɛ 18 - 2 ⁢ ɛ - 5 ⁢ ( 1 - ɛ ) 3 ⁢ ( 1 - ɛ ) ) ⁢ ( y 0 y 1 ( x 1 - x 0 ) ⁢ y 0 ′ ( x 1 - x 0 ) ⁢ y 1 ′ ) ( 2 )

8. A method in accordance with any one of claims 1-3, wherein said interpolation is calculated by dividing the interpolation interval into five parts, with the function values represented by ya, yb, yc, Yd and the derivative values represented by ya′, yb′, yc′, yd′ at the points where the interpolation interval (x0, x1) is divided, using the function values y0, y1 and function derivative values y0′, y1′,x0, x1, with the below-given Equation (3) as a function value interpolation equation, Equation (4) as a function derivative value interpolation equation and ε as a non-negative decimal fraction. ( y a y b y c y d ) = 1 125 ⁢ ( 112 - 12 ⁢ ɛ 13 + 12 ⁢ ɛ 16 ⁢ ( 1 - ɛ ) - 4 ⁢ ( 1 - ɛ ) 81 - 6 ⁢ ɛ 44 + 6 ⁢ ɛ 18 ⁢ ( 1 - ɛ ) - 12 ⁢ ( 1 - ɛ ) 44 + 6 ⁢ ɛ 81 - 6 ⁢ ɛ 12 ⁢ ( 1 - ɛ ) - 18 ⁢ ( 1 - ɛ ) 13 + 12 ⁢ ɛ 112 - 12 ⁢ ɛ 4 ⁢ ( 1 - ɛ ) - 16 ⁢ ( 1 - ɛ ) ) ⁢ ( y 0 y 1 ( x 1 - x 0 ) ⁢ y 0 ′ ( x 1 - x 0 ) ⁢ y 1 ′ ) ( 3 ) ( y a ′ y b ′ y c ′ y d ′ ) = 1 25 ⁢ ( x 1 - x 0 ) ⁢ ( - 24 - ɛ 24 + ɛ 8 ⁢ ( 1 - ɛ ) - 7 ⁢ ( 1 - ɛ ) - 36 + 11 ⁢ ɛ 36 - 11 ⁢ ɛ - 3 ⁢ ( 1 - ɛ ) - 8 ⁢ ( 1 - ɛ ) - 36 + 11 ⁢ ɛ 36 - 11 ⁢ ɛ - 8 ⁢ ( 1 - ɛ ) - 3 ⁢ ( 1 - ɛ ) - 24 - ɛ 24 + ɛ - 7 ⁢ ( 1 - ɛ ) 8 ⁢ ( 1 - ɛ ) ) ⁢ ( y 0 y 1 ( x 1 - x 0 ) ⁢ y 0 ′ ( x 1 - x 0 ) ⁢ y 1 ′ ) ( 4 )

9. A method in accordance with any one of claims 1-8, for increasing the speed of line-by-line calculations, when assuming the Voigt function to be K(x, y) and the shape of the absorption line displaced from the Voigt profile of the absorption line to be K(x, y)+f(x), by replacing K ⁡ ( x, y ) with ⁢: K ~ ⁡ ( x, y ) = A ⁢   ⁢ K ⁡ ( x, y ) + Bf ⁡ ( x ) and ∂ K ⁡ ( x, y ) ∂ x with ∂ K ~ ⁡ ( x, y ) ∂ x = A ⁢ ∂ K ⁡ ( x, y ) ∂ x + B ⁢ ∂ f ⁡ ( x ) ∂ x

10. A method in accordance with any one of claims 1-8, for increasing the speed of line-by-line calculations, when assuming the Voigt function to be K(x, y) and the shape of the absorption line displaced from the Voigt profile to be K(x, y)f(x), by replacing K ⁡ ( x, y ) with ⁢: K ~ ⁡ ( x, y ) = K ⁡ ( x, y ) ⁢ f ⁡ ( x ) and ∂ K ⁡ ( x, y ) ∂ x with ∂ K ~ ⁡ ( x, y ) ∂ x = ∂ K ⁡ ( x, y ) ∂ x ⁢ f ⁡ ( x ) + K ⁡ ( x, y ) ⁢ ∂ f ⁡ ( x ) ∂ x

11. A method in claim 10, for increasing the speed of line-by-line calculations, by using, in a sub-Lorentzian correction, K ~ ⁡ ( x, y ) = K ⁡ ( x, y ) ⁢ A ⁢   ⁢ exp ⁡ ( - B ⁢  x  ) and ∂ K ~ ⁡ ( x, y ) ∂ x = ∂ K ⁡ ( x, y ) ∂ x ⁢ A ⁢   ⁢ exp ⁡ ( - B ⁢  x  ) + K ⁡ ( x, y ) ⁡ [ - sgn ⁡ ( x ) ⁢ AB ⁢   ⁢ exp ⁡ ( - B ⁢  x  ) ]

12. A method in accordance with claim 9, for increasing the speed of line-by-line calculations, by replacing, in a line-mixing correction, K ⁡ ( x, y ) with ⁢: K ~ ⁡ ( x, y ) = A ⁢   ⁢ K ⁡ ( x, y ) + B ⁢   ⁢ L ⁡ ( x, y ) and ∂ K ⁡ ( x, y ) ∂ x with ∂ K ~ ⁡ ( x, y ) ∂ x = - 2 ⁡ [ ( A ⁢   ⁢ x + B ⁢   ⁢ y ) ⁢ K ⁡ ( x, y ) - ( A ⁢   ⁢ y - B ⁢   ⁢ x ) ⁢ L ⁡ ( x, y ) - B π ] Here, L(x, y) is the imaginary component of the function w(z) (the real part is the Voigt function), where complex number z=x+iy, defined by the following equation. w ⁡ ( z ) = i π ⁢ ∫ - ∞ ∞ ⁢ exp ⁡ ( - t 2 ) z - t ⁢   ⁢ ⅆ t ⁢   = exp ⁡ ( - z 2 ) ⁢ erfc ⁡ ( - i ⁢   ⁢ z ) = K ⁡ ( x, y ) + i ⁢   ⁢ L ⁡ ( x, y ) (erfc(z) is a complex complementary error function).

13. A program for increasing the speed of line-by-line calculations of multiple overlapping absorption lines, said program comprising:

(1) a step of dividing a domain of a Voigt function representing the shape of an absorption line into a first range around the peak of the Voigt function and a skirt portion not contained in the first range, replacing the first range with a cubic function whose function values and derivative values match with those of said Voigt function at the connection points with said Voigt function, and calculating the values and derivative values of said cubic function and the Voigt function in the skirt portion for each of first predetermined intervals;

(2) a step of adding together the results of step (1) for a plurality of absorption lines;

(3) a step of interpolating the results of step (2) to calculate the function values and derivative values in second predetermined intervals smaller than said first predetermined intervals;

(4) a step of dividing said first range into a second range near the peak and a skirt portion not contained in the second range, replacing said second range of a “function representing the difference between the Voigt function and said cubic function” with a second cubic function whose function values and derivative values match with those of said “function representing the difference between the Voigt function and said cubic function” at the connection points with said “function representing the difference between the Voigt function and said cubic function”, and calculating the values and derivative values of said second cubic function replacing said second range and the “function representing the difference between the Voigt function and said cubic function” in the skirt portion for said second predetermined intervals;

(5) a step of adding the results of step (4) to the results of step (3) for a plurality of absorption lines;

(6) a step of interpolating the results of step (5) to calculate the function values and derivative values for intervals smaller than said second predetermined intervals; and

(7) a step of adding the values and derivative values of the “function representing the difference between the ““function representing the difference between the Voigt function and said cubic function ”” and said second cubic function” to the results of step (6) for a plurality of absorption lines in said second range.

14. A program in accordance with claim 13, wherein the step of calculating function values and derivative values by interpolation is repeated while narrowing the interval until it becomes an interval of a minimum unit.

15. A program in accordance with claim 13, wherein said steps (4) through (6) are repeated one or more times until the “intervals smaller than said second predetermined intervals” in step (6) reach third predetermined intervals.