1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
|
# VOIGT -- Compute the real (Voigt function) and imaginary parts of the
# complex function w(z)=exp(-z**2)*erfc(-i*z) in the upper half-plane
# z=x+iy. The maximum relative error of the real part is 2E-6 and the
# imaginary part is 5E-6.
#
# From: Humlicek, J. Quant. Spectrosc. Radiat. Transfer, V21, p309, 1979.
procedure voigt (xarg, yarg, wr, wi)
real xarg #I Real part of argument
real yarg #I Imaginary part of argument
real wr #O Real part of function
real wi #O Imaginary part of function
int i
real x, y, y1, y2, y3, d, d1, d2, d3, d4, r, r2
real t[6], c[6], s[6]
data t/.314240376,.947788391,1.59768264,2.27950708,3.02063703,3.8897249/
data c/1.01172805,-.75197147,1.2557727e-2,1.00220082e-2,-2.42068135e-4,
5.00848061e-7/
data s/1.393237,.231152406,-.155351466,6.21836624e-3,9.19082986e-5,
-6.27525958e-7/
begin
x = xarg
y = abs (yarg)
wr = 0.
wi = 0.
y1 = y + 1.5
y2 = y1 * y1
# Region II
if (y < 0.85 && abs(x) > 18.1*y+1.65) {
if (abs(x) < 12)
wr = exp (-x * x)
y3 = y + 3
do i = 1, 6 {
r = x - t[i]
r2 = r * r
d = 1 / (r2 + y2)
d1 = y1 * d
d2 = r * d
wr = wr + y * (c[i] * (r * d2 - 1.5 * d1) + s[i] * y3 * d2) /
(r2 + 2.25)
r = x + t[i]
r2 = r * r
d = 1 / (r2 + y2)
d3 = y1 * d
d4 = r * d
wr = wr + y * (c[i] * (r * d4 - 1.5 * d3) - s[i] * y3 * d4) /
(r2 + 2.25)
wi = wi + c[i] * (d2 + d4) + s[i] * (d1 - d3)
}
# Region I
} else {
do i = 1, 6 {
r = x - t[i]
d = 1 / (r * r + y2)
d1 = y1 * d
d2 = r * d
r = x + t[i]
d = 1 / (r * r + y2)
d3 = y1 * d
d4 = r * d
wr = wr + c[i] * (d1 + d3) - s[i] * (d2 - d4)
wi = wi + c[i] * (d2 + d4) + s[i] * (d1 - d3)
}
}
end
|
