;+
; Project     : YOHKOH-BCS
;    
; Name        : PVOIGT
;
; Purpose     : Compute voigt function with partial derivatives
;
; Explanation : Calculates the Voigt function after  Humlicek (1982) 
;               JQSRT 27,437. Derivatives are from Humlicek (1979) 
;               JQSRT 21, 309. This method is fast and VERY accurate.
;
; Category    : fitting
;
; Syntax      : pvoigt,a,v,h,f
;
; Inputs      : a = damping parameter
;               v = frequency vector
;
; Outputs     : h=H(a,v)
;               f=2F(a,v)
;
;               dH(a,v)/dv = 2.0(af-vh) 
;               dH(a,v)/da = 2.0(ah+vf-.56418958)
;
;               d2 H(a,v)/d2v = 4[(v^2-a^2)h - 2avf - h/2 + a*0.56418958]
;
;
;               d(2F)/dv = -2.0(ah+vf-.5641896) (Cauchy-Riemann equation)
;
;               d2 (2F(a,v))/d2v = -4[(a^2-v^2)f - 2avh + f/2 +v*0.56418958]
;
;               where .5641896 is 1/sqrt(pi)
;
; Opt Outputs : None
;
; Keywords    : None
;
; Restrictions: None
;
; Side effects: None
;
; History     : Version 1,  17-July-1993,  D M Zarro.  Written
;               Modified from T.R Metcalf (IAH)
;
; Contact     : DZARRO@SOLAR.STANFORD.EDU
;-
 
  pro pvoigt, a, v, h, f

  check = check_math(1,1)   ; turn off math errors to avoid underflow reports

  sizev = size(v)

  if sizev(0) GT 1 then begin
     print,'ERROR: voigt: v has too many dimensions: scalar or vector only'
     return
  endif

  t = complex(float(a),float(-v))
  w1 = complexarr(sizev(n_elements(sizev)-1))
  u = t*t
  absv = abs(v)
  s = absv+abs(a)
  s2 = 0.195*absv-0.176

  region1 = where(s GE 15)
  region2 = where(s GE 5.5 and s LT 15)
  region3 = where(s LT 5.5 and abs(a) GE s2) 
  region4 = where(s LT 5.5 and abs(a) LT s2)
 
  ;if region1(0) GE 0 then print,'region1: ',n_elements(region1)
  ;if region2(0) GE 0 then print,'region2: ',n_elements(region2)
  ;if region3(0) GE 0 then print,'region3: ',n_elements(region3)
  ;if region4(0) GE 0 then print,'region4: ',n_elements(region4)

; Region 1

  if region1(0) GE 0 then $
     w1(region1) = t(region1)*(1.4104739589+u(region1)* $
                   0.56418958355)/(0.74999999999+u(region1)*(3.0+u(region1)))

; Region 2

  if region2(0) GE 0 then $
     w1(region2) = t(region2)*(4.6545640642+u(region2)*(3.9493270848+u(region2)*0.56418958355)) / $
                   (1.8750000000+u(region2)*(11.250000000+u(region2)*(7.5000000000+u(region2))))

; Region 3

  if region3(0) GE 0 then $
     w1(region3) = (179.0714766+t(region3)*(289.4827444+t(region3)*(231.1440618+t(region3)*(111.5306766+t(region3)*(34.03818681+t(region3)*(6.269921295+t(region3)*0.5641900381)))))) / $
     (179.0714766+t(region3)*(491.5434301+t(region3)*(606.7189137+t(region3)*(439.3066119+t(region3)*(203.2455046+t(region3)*(60.83060715+t(region3)*(11.11317476+t(region3))))))))

; Region 4

  if region4(0) GE 0 then $
     w1(region4) =  (exp(u(region4))-t(region4)*(36183.30536-u(region4)*(3321.990492-u(region4)*(1540.786893-u(region4)*(219.0312964-u(region4)*(35.76682780-u(region4)*(1.320521697-u(region4)*0.5641900381)))))) / $
       (32066.59372-u(region4)*(24322.84021-u(region4)*(9022.227659-u(region4)*(2186.181081-u(region4)*(364.2190727-u(region4)*(61.57036588-u(region4)*(1.841438936-u(region4)))))))))

  check = check_math(0,0)   ; turn on math errors and check if any occurred
  if (check NE 0) and (check NE 32) then begin
    if (check and 1) GT 0 then $
       print,'WARNING: math error in voigt.pro: Integer divide by zero'
    if (check and 2) GT 0 then $
       print,'WARNING: math error in voigt.pro: Integer overflow'
    if (check and 16) GT 0 then $
       print,'WARNING: math error in voigt.pro: Floating point divide by zero'
    ; skip 32 since that is the underflow we are trying to avoid
    if (check and 64) GT 0 then $
       print,'WARNING: math error in voigt.pro: Floating point overflow'
    if (check and 128) GT 0 then $
       print,'WARNING: math error in voigt.pro: Floating point operand error'
  endif

  h = float(w1)      ; The real part
  f = imaginary(w1)

  if (sizev(0) EQ 0) then begin    ; v is a scalar in this case - not an array
    h=h(0)  ; Return a scalar not an array in this case
    f=f(0)
  endif

end