; file: exthneubf_gray.pro = HI bound-free following Gray 1992, 2005 ; init: Jul 8 2010 ; last: Oct 2 2013 Rob Rutten function exthneubf_gray,wav,temp,precision=precision ;+ ; ext=exthneubf_gray(wav,temp,precision=precision) ; in: wavelength [Angstrom] number or array but not both arrays ; temperature [K] number or array but not both arrays ; precision (default 1E-3) ; out: sigma H_bf cm^2 per neutral H atom assuming Saha-Boltzmann ; NB: includes stimulated-emission correction already ;- ; answer no-parameter query if(n_params(0) lt 2) then begin print,'ext=exthneubf_gray(wav,temp,precision=precision)' print,' wav = number or array [AA]' print,' temp = number or array [K]' print,' but not both arrays, that would be asking too much' print,' precision=precision (default 1E-3)' return,-1 endif ; keyword defaults if (n_elements(precision) eq 0) then precision=1.E-3 ; read physics constants (cgs) @constants_cgs.idl ; bound-free: Gray but with explicit summation loop over levels ;RR essential to make level and wav double precision floats lambda=1.D*float(wav) nlevels=200 edgewav=fltarr(nlevels+1) boltzpop=fltarr(nlevels+1) parfunh=partfunc_rr(temp,1,0) sigmabf=0 ; initiate summation ; loop over hydrogen edges beyond this wavelength for ilevel=1,nlevels do begin ;RR real n, not IDL index with zero level=ilevel*1.D0 edgewav(ilevel)=level^2/rydbergaa ; Gray 2005 (8.5) gauntbf=1.0-0.3456/ $ ((lambda*rydbergaa)^0.3333)*(lambda*rydbergaa/level^2-0.5) ; set gauntbf to zero for wavelengths beyond this edge (with Krijger trick) if (size(lambda,/n_elements) gt 1) then begin if (max(lambda) ge edgewav(ilevel)) then $ gauntbf[where(lambda gt edgewav(ilevel))]=0. endif else if (lambda ge edgewav(ilevel)) then gauntbf=0. ; summation with Boltzmann populations as Gray 2005 8.7 2nd version glevel=2.*level^2 ; statistical weight chiexc=rydbergev*(1.-1./(level^2)) ; in eV add=1.0449E-26*lambda^3/level^3*gauntbf*exp(-chiexc/(kev*temp)) sigmabf=sigmabf+add if (max(add/sigmabf) lt precision/10.) then break ;; print,' --- ',ntostr([ilevel,edgewav(ilevel),gauntbf,glevel,sigmabf]) endfor ; correction stimulated emission although this is still per neutral H atom alphabf=sigmabf*(1.-exp(-hcgs*ccgs/(lambda*1E-8*kcgs*temp))) ;; print,' ----- ', alphabf,alphabf/sigmabf ; STOP return,alphabf end