; $Id: cerrorf.pro,v 1.2 2002/01/14 20:59:26 amunter Exp $
FUNCTION cerrorf, z, xlimit=xlimit, ylimit=ylimit
;+
; NAME:		cerrorf
; PURPOSE: This is the complex error function:
;
;	w(z) = exp(-z^2)*erfc(-iz)
;
;	See: Nuclear Instruments and Methods Vol. 137 (1976) 179-181
;	     for its definition and use.
;	See: Communications of the ACM Vol 12, No. 11, (1969) 635 for the
;	     algorithm used.
;
; CATEGORY: "simple" function?
; CALLING SEQUENCE:	result=cerrorf(z)
; INPUTS: z = any real or complex number or array
; OPTIONAL INPUT PARAMETERS:
;	xlimit, ylimit define a box (with one corner at the origin)
;	in the first complex quadrant, outside which
;	a simple Gauss-Hermite quadrature algorithm is used, rather than the full 
;	algorithm of continued fractions. See the above references for more info.
;	If not set, these default to xlimit=5.33 and ylimit=4.29
; COMMON BLOCKS: none
; SIDE EFFECTS: no idea
; RESTRICTIONS: haven't got a clue yet.
; PROCEDURE:
;	If you are trying to fit data, and your fitting function is the
;	convolution of a lorentzian (your intrinsic line shape) and a 
;	gaussian (your instrumental resolution) then this is what you need.
;
;	Your fitting function is then:
;
;		D(x) = Imax*(GammaL/GammaG)*Sqrt(Pi*ln2)Re[w(v+ia)]
;	where
;		D(x) is the intensity
;	   	  Imax is obvious
;		  GammaL = Full Width Half Max. (FWHM) of the Lorentzian
;		  GammaG = FWHM of the Gaussian (=2*Sqrt(2*ln2)*sigma)
;	   	  a  = (GammaL/GammaG)Sqrt(ln2)
;	   	  v  = (x-x0)2Sqrt(ln2)/GammaG
;	   	  x0 = center of the lorentzian peak.
;
;	The derivatives are also extremely easy to get, when you are trying to
;	minimize Chi-squared:
;
;		Re[dw(z)/da] = -2/Sqrt(Pi)+2vIm[w(z)]+2aRe[w(z)]
;		Re[dw(z)/dv] = -2vRe[w(z)]+2aIm[w(z)]
; MODIFICATION HISTORY:
;	Paul Schleger, Riso National Lab., September 1993 	
;
;-
IF keyword_set(xlimit) THEN xlimit = xlimit ELSE xlimit= 5.33
IF keyword_set(ylimit) THEN ylimit = ylimit ELSE ylimit = 4.29
;****** Get Real and Imaginary parts ******
x=FLOAT(z)
y=IMAGINARY(z)
;****** Find out where we are not in the first quadrant ******
xandyneg = where((x LT 0) AND (y LT 0),xycount)
xneg = where((x LT 0) AND (y GE 0),xcount)
yneg = where((x GE 0) AND (y LT 0),ycount)
;****** Just do the calculation in the 1st Quadrant, then map back at end *****
x=ABS(x)
y=ABS(y)
;****** Initialize some local variables ******
s=x*0. & capn=LONG(s) & nu = LONG(s)+8 & np1 = LONG(s) & lambda = s
r1 = s & r2 = s & s1 = s & s2 = s & np1 = s
t1 = s & t2 = s & c = s
;**** Find out where the "box" region
index = WHERE((y LT ylimit) AND (x LT xlimit),count)
IF (count NE 0) THEN s(index) = (1.-y(index)/ylimit)*SQRT(1.-(x(index)/xlimit)^2)
h = 1.6*s
h2 = 2.*h
IF (count NE 0) THEN BEGIN
	capn(index) = LONG(6 + 23.*s(index))
	nu(index) = LONG(9. + 21.*s(index))
END
lambda = (h NE 0)*h2^capn
b = ((h EQ 0) OR (lambda EQ 0))

FOR n=MAX(nu),0,-1 DO BEGIN
	i = WHERE(n LE nu)
	np1(i) = DOUBLE(n + 1)
	t1(i) = y(i) + h(i) + np1(i)*r1(i)
	t2(i) = x(i) - np1(i)*r2(i)
	c(i) = 0.5/(t1(i)*t1(i)+t2(i)*t2(i))
	r1(i) = c(i)*t1(i)
	r2(i) = c(i)*t2(i)
	i= WHERE((h GT 0) AND (n LE capn),count)
	IF (count NE 0) THEN BEGIN
		t1(i) = lambda(i)+s1(i) 
		s1(i) = r1(i)*t1(i) - r2(i)*s2(i)
		s2(i) = r2(i)*t1(i) + r1(i)*s2(i)
		lambda(i) = lambda(i)/h2(i)
	END
END
echkx = WHERE(ABS(x) LE 9.00,xchkcount)
xandyneg = WHERE(ABS(z) LE 9.00 AND xandyneg,xycount)
yneg = WHERE(ABS(z) LE 9.00 AND yneg,ycount)
xneg = WHERE(ABS(z) LE 9.00 AND xneg,xcount)
rew=fltarr(n_elements(x))
IF (xchkcount NE 0) THEN rew(echkx) = (y(echkx) EQ 0.)*exp(-x(echkx)^2) 
rew = rew+ (y NE 0.)*1.12837916709551*((b NE 0)*r1+(b EQ 0)*s1)
imw = 1.12837916709551*((b NE 0)*r2+(b EQ 0)*s2)
w = COMPLEX(rew,imw)
;****** Have the folowing rules to figure out other quadrant values:
;	w(-z) = 2exp(-z^2)-w(z)
;	w(CONJ(z)) = CONJ(w(-z))
;****** So do it! ******
IF (xycount NE 0) THEN $
	w(xandyneg) = 2.*exp(-(COMPLEX(x,y))(xandyneg)^2)-w(xandyneg)
IF ( ycount NE 0) THEN $
	w(yneg) = CONJ(2.*exp(-(COMPLEX(x,y))(yneg)^2)-w(yneg))
IF ( xcount NE 0) THEN $
	w(xneg) = CONJ(w(xneg))
RETURN, w
END
