1 files changed, 71 insertions, 0 deletions
diff --git a/noao/onedspec/splot/voigt.x b/noao/onedspec/splot/voigt.x
new file mode 100644
index 00000000..08a44c78
--- /dev/null
+++ b/noao/onedspec/splot/voigt.x
@@ -0,0 +1,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