; Following Drayson et al. 1976. "Rapid computation of the Voigt profile." JSQRT 16, 661-614.
;
; Added check for y = 0 (pure doppler)
;
; Calling sequence is the same as for the IDL fumction VOIGT.
;   a= al/ad
;   u = double(dnu/ad)
; To normalize, 
;   voig=(1/(ad*sqrt(!pi)))*goodvoigt4(a, u)
;

function goodvoigt4, y, x2

common share_goodvoigt4, b,h,xn,yn,xx,hh,nby2,c,ri,hn,d0,d1,d2,d3,d4

if n_elements(x2) ne 1 then begin
  n = n_elements(x2)
  VGT2list = dblarr(n)
  for i=0d,n-1 do begin
    VGT2list[i] = goodvoigt4(y, x2[i])
  end
  return, VGT2list
end

if n_elements(hn) ne 25 then begin

     b=dindgen(22)
     b(1-1)=0.0d
     b(2-1)=0.00000007093602d
     h=0.201
	xn = [10.0, 9.0, 8.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0d]
	yn = [ 0.6, 0.6, 0.6, 0.5, 0.4, 0.4, 0.3, 0.3, 0.3, 0.3, 1.0, 0.9, 0.8, 0.7, 0.7d]
	nby2 =[9.5, 9.0, 8.5, 8.0, 7.5, 7.0, 6.5, 6.0, 5.5, 5.0, 4.5, 4.0, 3.5, 3.0, 2.5, $
	 2.0, 1.5, 1.0, 0.5d]
	xx = [0.5246476, 1.65068, 0.7071068d]
	hh = [0.2562121, 0.02588268, 0.2820948d]
	c = [0.00000007093602, -0.0000002518434, 0.0000008566874, -0.000002787638, $
	0.000008660740, -0.00002565551, 0.00007228775, -0.0001933631, $
	0.0004899520, -0.001173267, 0.002648762, -0.005623190, 0.01119601, $
	-0.02084976, 0.03621573, -0.05851412, 0.08770816, -0.121664, 0.15584, $
	-0.184, 0.2d]
	
	ri=dindgen(15)
	for i=1,15 do begin
		ri[i-1] = -i/2.
	endfor
	
	hn=dindgen(25)
	d0=dindgen(25)
	d1=dindgen(25)
	d2=dindgen(25)
	d3=dindgen(25)
	d4=dindgen(25)
	for i=1,25 do begin
		hn[i-1]=h*(i-0.5d)
		co=4.0*hn[i-1]*hn[i-1]/25.0d -2.0d
	
		for j=2,21 do begin
			b[j+1-1]=co*b[j-1]-b[j-2]+c[j-1]
		endfor
		
		d0[i-1]=hn[i-1]*(b[22-1]-b[21-1])/5.d
		d1[i-1]=1.0-2.0*hn[i-1]*d0[i-1]
		d2[i-1]=(hn[i-1]*d1[i-1]+d0[i-1]-1.d)/ri[2-1]
		d3[i-1]=(hn[i-1]*d2[i-1]+d1[i-1]-1.d)/ri[3-1]
		d4[i-1]=(hn[i-1]*d3[i-1]+d2[i-1]-1.d)/ri[4-1]
	endfor
              
end

x = abs(x2)

y2=y*y

; Test for y = 0
if (y eq 0) then begin
  VGT2 = exp(- ((x^2) < 700.d) )
  return, VGT2
end

; Test for region IIIB
if (y ge 11.0d -0.6875d * x) then begin
	u=x-xx[3-1]
	v=x+xx[3-1]
	VGT2=y*(hh[3-1]/(y2+u*u)+hh[3-1]/(y2+v*v))
	return, VGT2
endif

; Test for region IIIA
if( (y gt 1.d and x gt 1.85d * (3.6d - y)) or (x gt 5.0d - 0.8d * y) ) then begin
	u=x-xx[1-1]
	v=x+xx[1-1]
	uu=x-xx[2-1]
	vv=x+xx[2-1]
	VGT2=y*(hh[1-1]/(y2+u*u)+hh[1-1]/(y2+v*v)+hh[2-1]/(y2+uu*uu)+hh[2-1]/(y2+vv*vv))
	return, VGT2
endif

; Test for region II
if (y ge 1.d) then begin
; Continued fraction.  Compute number of terms needed
	if(y ge 1.45d) then begin
		i=y+y
	endif else begin
		i=11.0d * y
	endelse
	
	j=floor(x+x+1.85d)

	imax=floor(xn[j-1]*yn[i-1]+0.46d)
	imin= 16 < (21-2*imax)

; Evaluate continued fraction
	uu=y
	vv=x
	for j=imin, 19	do begin
		u=nby2[j-1]/(uu*uu+vv*vv)
		uu=y+u*uu
		vv=x-u*vv
	endfor

	VGT2=uu/(uu*uu+vv*vv)/1.772454d
	return, VGT2
endif

; All that remains must be in region I
; Region I:  Dawson's function at x from Taylor series

n=floor(x/h)
dx=x-hn[n+1-1]
u=(((d4[n+1-1]*dx+d3[n+1-1])*dx+d2[n+1-1])*dx+d1[n+1-1])*dx+d0[n+1-1]
v=1.0d - 2.0d*x*u

; Taylor series expansion about y=0.0

vv=exp(y2-x*x)*cos(2.0d*x*y)/1.128379d - y*v
uu=-y
imax=5.0+(12.5d - x)*0.8d*y

for i=2,imax,2 do begin
	u=(x*v+u)/ri[i-1]
	v=(x*u+v)/ri[i+1-1]
	uu = -uu*y2
	vv=vv+v*uu
endfor

VGT2=1.128379d*vv
return, VGT2
end