;+
; NAME:
;	approxwidth.pro
; PURPOSE: 
;	Calculate equivalent widths assuming appoximate forms of 
;	Doppler broadening, Lorentx broadening and Voigt profile. 
;	The equations are from Ch. 4 of Goody & Yung.
; DESCRIPTION:
;	a function to help determine the approx. equivalent width
;	of spectral lines. This was used to check widthv.pro and
;	widtha.pro. 
; CATEGORY: 
; CALLING SEQUENCE: approxwidth, S, N, adop, aL, amf, Wd, Wl, Wv
; INPUTS: 
;	S	-- strength of the line shape, given in cm^-1/(molecules*cm^-2)
;	N	-- Column abundance of molecules which cause the feature,
;		   given in molecules/cm^2.
;	adop	-- Doppler line shape, given in cm^-1
;	aL	-- Lorentzian line shape, given in cm^-1
;	amf	-- air mass factor
; OPTIONAL INPUT PARAMETERS: none
; KEYWORD PARAMETERS:
; OUTPUTS: 
;	Wd	-- the approx. equivalent width using a Doppler profile.
;	Wl	-- the approx. equivalent width using a Lorentz profile.
;	Wv	-- the approx. equivalent width using a Voigt profile.
; SIDE EFFECTS: none
; RESTRICTIONS: 
; PROCEDURE: 
;  Wd from Goody & Yung 4.20
; MODIFICATION HISTORY: Written by Jason Cook, May 4, 1999, BU
; 2009 Nov 25 LAY, SwRI. 
;   Remove print statement, 
;   multiply by amf for optically-thin case
;   replace !pi w/ sqrt(!pi) for w
; 2009 Dec 5 LAY
;   change low-w doppler case from S*m for w<1 to
;   Eq. 4.19 for w < 3.
;-

function approxwidth_dop, S, N, adop, amf

m = N * amf

w=S*m/(adop*sqrt(!pi)) 
;print, w

;if alog(w) lt 0 then begin 
if w lt .1 then begin 

   sum = 1.d
   fact = 1.d
   wi = 1.d
   i = 1.d
   repeat begin
      fact *= i
      wi *= (-w)
      addend = (wi)/(fact*sqrt(i))
      sum += addend
;      sum += ((-1)^i)*(w^i)/(factorial(i)*sqrt(i))
   endrep until i ge 30 or abs(addend/sum) lt 1e-6
   Wd = sqrt(!pi) * w * adop * sum
   
endif else begin

   if w lt 2.5 then begin
      ; 4.23 with al = 0
      nu0 = 1.25 * adop
      Wd = 2*nu0*(1 - exp(-S*m/(2*nu0) ) )            

   endif else begin
      lnw = alog(w)
      Wd=2*adop* (lnw^(1/2.)+0.2886*lnw^(-1/2.)-0.2473*lnw^(-3/2.)+0.3403*lnw^(-5/2.))
   endelse
endelse

return, Wd

end