; Lineshapes
;
; INPUTS
;  ad = doppler width (cm^-1)
;  al = lorentz width (cm^-1)
;  dnu = distance from line center (cm^-1)
;  
; OUTPUTS
;  lineshape, in 1/(cm^-1)
;
; DIMENSIONS
;  If ad, al, and dnu are scalers, lineshape* returns a scalar
;  If ad and al are non scalar, they must have the same dimension and dnu must be scalar
;  If dnu is nonscalar, then ad and al must be scalars
;
; RESTRICTIONS
;  This could do a more careful job of checking against underflow
; MODIFICATION HISTORY:
;  2004 June LAY call voigt, not goodvoigt4; a is double


function lineshape, ad, al, dnu

   a= double(al/ad)
   u = double(dnu/ad)
   voig=(1/(ad*sqrt(!pi)))*voigt(a, u)
   return, voig

end

function lineshape_dop, ad, dnu

  x = double(dnu/ad)
  return, 1/(ad * sqrt(!pi)) * exp( -x^2)

end

function lineshape_lor, al, dnu

  return, (al/!pi) * (1/(dnu^2 + al^2))
end

pro test_lineshape

print, '% Testing lineshape'

print, '%  test #1: plot the doppler and lorentz line shapes'
  ad = 1.
  al = 1.
  dnu = findgen(601)/100.-3
  f_lor = lineshape_lor(al, dnu)
  f_dop = lineshape_dop(ad, dnu)
  yr = [0, max(f_dop)]
  plot, dnu, f_dop, xr=[-3,3], /xs, yr=yr, $
     xtit='(nu-nu0)/a', ytit='f(nu-nu0)', tit='Following Goody & Yung Fig 2.31'
  oplot, dnu, f_lor
  xyouts, 0, yr[1], 'Doppler', ali=0.5
  xyouts, 0, max(f_lor), 'Lorentz', ali=0.5
  read, dummy, prompt='(hit return to continue)', format='(A)'
  
print, '% test #2: Plotting Voigt, Lorentz, and Doppler for ad << al'
  ad = 1.e-1
  al = 1.
  dnu = findgen(41)/10.
  f_lor = lineshape_lor(al, dnu)
  f_dop = lineshape_dop(ad, dnu/10.)
  f = lineshape(ad, al, dnu)
  yr = [0,max(f)*1.1]
  plot, dnu, f, xr=[0,4], /xs, yr=yr, /ys, $
     xtit='(nu-nu0)/a', ytit='f(nu-nu0)', tit='ad = 0.1, al = 1'
  oplot, dnu/10, f_dop*(max(f)/max(f_dop)), li=1
  oplot, dnu, f_lor, ps=4
  xyouts, 0.3, max(f), $
          'Following Goody & Yung discussion, p. 112 for ad<<al!c'+$
          'Doppler (dotted, scaled) is a delta function,!c' + $
          'so the Voigt profile (solid) resembles Lorentz (diamonds) everywhere' 
  read, dummy, prompt='(hit return to continue)', format='(A)'

print, '% test #3: Plotting Voigt, Lorentz, and Doppler for ad >> al'
  ad = 1.
  al = 1.e-2
  dnu = findgen(601)/100.
  f_lor = lineshape_lor(al, dnu)
  f_dop = lineshape_dop(ad, dnu)
  f = lineshape(ad, al, dnu)
  yr = max(f)*[1e-4,1]
  plot, dnu, f, xr=[0,max(dnu)], /xs, yr=yr, /ylo, $
     xtit='(nu-nu0)/a', ytit='f(nu-nu0)', tit='ad = 1, al = 0.01'
  oplot, dnu, f_dop, li=1
  oplot, dnu, f_lor, li=2
  xyouts, 1.3,0.5, $
          'Following Goody & Yung discussion, p. 113 for al<<ad!c'+$
          'Voigt profile (solid) resembles Doppler (dotted) near line center, !c'+$
          'and Lorentz (dashed) in the wings' 
  read, dummy, prompt='(hit return to continue)', format='(A)'

print, '% test #4: Checking normalization for al >> ad'
  ad = 1.e-2
  al = 1.
  nuwgt_overkill, 0.001, ad, al, dnu, nuwgt
  f_dop = lineshape_dop(ad, dnu)
  f_lor = lineshape_lor(al, dnu)
  f = lineshape(ad, al, dnu)
  print, ' totals for Doppler, Lorentz, and Voigt = ', $
    [ total( f_dop * nuwgt), total( f_lor * nuwgt), total( f * nuwgt) ]
  
print, '% test #5: Checking normalization for ad >> al'
  ad = 1.
  al = 1.e-2
  nuwgt_overkill, 0.001, ad, al, dnu, nuwgt
  f_dop = lineshape_dop(ad, dnu)
  f_lor = lineshape_lor(al, dnu)
  f = lineshape(ad, al, dnu)
  print, ' totals for Doppler, Lorentz, and Voigt = ', $
    [ total( f_dop * nuwgt), total( f_lor * nuwgt), total( f * nuwgt) ]
  
end