Viewing contents of file '../idllib/ssw/allpro/acgaunt.pro'
function acgaunt, wave, te_6, $
G1=gaunt_ff, G2=gaunt_fb, G3=gaunt_2p
;+
; NAME: ACGAUNT
; PURPOSE: Calculate continuum gaunt factor using approximations
; of R. Mewe (18-JUN-85) to full calculations of
; paper VI (Arnaut and Rothenflug for ion balances).
; CATEGORY:
; CALLING SEQUENCE: cgauntf = acgaunt( wave, te_6)
; cgauntf = acgaunt( wave, te_6, G1=gff, G2=gfb, G3=g2p)
;
; INPUTS: wave = Wavelength in Angstrom (1-d vector or scalar)
; te_6 = Temperature in 10^6 K (1-d vector or scalar)
; OPTIONAL INPUTS: none.
; OUTPUTS: Function result
; = cgauntf(n_elements(te_6), n_elements(wave))
; = array of approximate continuum gaunt factors.
; OPTIONAL OUTPUTS:
; G1 = Free-free Gaunt factor
; G2 = Free-bound Gaunt factor
; G3 = 2-photon Gaunt factor
;
; COMMON BLOCKS: none
; SIDE EFFECTS: none
; RESTRICTIONS: none
; PROCEDURE: see a paper of R. Mewe et al. (A & Ap)
; MODIFICATIONS: written by N.Nitta from a Fortran version, March 1991.
; 31-jul-93, JRL, Added a check on the exponent to prevent floating underflow message
; 23-Jun-94, DMZ, made Gaunt factors double precision
;
;-
;
; elements of te_6 and wave
;
nte_6=n_elements(te_6)
if(nte_6 le 1) then te_6=te_6*replicate(1.,1)
nwave=n_elements(wave)
if(nwave le 1) then wave=wave*replicate(1.,1)
;
; constituents of gaunt factor as arrays
;
gaunt_2p=dblarr(nte_6,nwave)
gaunt_ff=dblarr(nte_6,nwave)
gaunt_fb=dblarr(nte_6,nwave)
;
;
; Local variables for 2 photon calculation
;
n_edg = 6 ; Number of edges
lamda = fltarr(n_edg) ; Wavelength of edges
temp = fltarr(n_edg) ; Max temperature of edges
CC = fltarr(n_edg) ; Noramlization coefficient
DD = fltarr(n_edg) ; exponent coefficient
; OVIII, OVII, NVII, NVI , CVI, CV
lamda=[ 19, 22, 25, 29, 34, 41]
temp= [2.45, 0.9, 1.7, 0.6, 1.05, 0.37]
CC = [ 3.2, 11.3, 0.69, 2.1, 3.6, 10.3]
DD = [ 10.3, 2.5, 10.3, 2.5, 10.3, 2.5]
;
; Local variables for Free-bound calculation for low temperatures
n_t_l = 5 ; Number of temperature regions
n_w_l = 4 ; Number of wavelength regions
; integer*4 il,jl ! Column and row number of case table
al = fltarr(n_t_l,n_w_l) ; al coefficients
bl = fltarr(n_t_l,n_w_l) ; bl coefficients
cl = fltarr(n_t_l,n_w_l) ; cl coefficients
dl = fltarr(n_t_l,n_w_l) ; dl coefficients
tem_liml = fltarr(n_t_l-1) ; Limits of temperature cases
wav_liml = fltarr(n_w_l-1) ; Limits of wavelength cases
; Normalization Coefficients
al(*,0) = [ .248, 5.42e-12, 3.68e-4, 1.86e+3, 6.5e-3]
al(*,1) = [ .248, 5.42e-12, 3.68e-4, .176, .176]
al(*,2) = [ .248, 5.42e-12, .323, .323, .323]
al(*,3) = [.0535, .0535, .0535, .0535, .0535]
; Exponent Coefficients
bl(*,0) = [-1., -9.39, -4.9, -.686, -5.41]
bl(*,1) = [-1., -9.39, -4.9, -1., -1.]
bl(*,2) = [-1., -9.39, -1., -1., -1.]
bl(*,3) = [-1., -1., -1., -1., -1.]
; Exponent Coefficients
cl(*,0) = [.158, 0.0, 0.0, 0.0, 0.0]
cl(*,1) = [.158, 0.0, 0.0, .233, .233]
cl(*,2) = [.158, 0.0, .16, .16, .16]
cl(*,3) = [.046, .046, .046, .046, .046]
; Exponent Coefficients
dl(*,0) = [-1., 0.0, 0.0, 0.0, 0.0]
dl(*,1) = [-1., 0.0, 0.0, -1., -1.]
dl(*,2) = [-1., 0.0, -1., -1., -1.]
dl(*,3) = [-1., -1., -1., -1., -1.]
tem_liml=[.015, .018, .035, .07] ; Boundaries btwn temp. bins
wav_liml=[227.9, 504.3, 911.9] ; Boundaries btwn wave bins
;#####################################################################
;
; Local variables for Free-bound calculation for high temperatures
n_t = 10 ; Number of Temperature regions
n_w = 16 ; Number of wavelength regions
a = fltarr(n_t,n_w) ; a coefficients
b = fltarr(n_t,n_w) ; b coefficients
c = fltarr(n_t,n_w) ; c coefficients
d = fltarr(n_t,n_w) ; d coefficients
tem_lim = fltarr(n_t-1) ; Limits of temperature cases
wav_lim = fltarr(n_w-1) ; Limits of wavelength cases
; Normalization Coefficients
dumarr=fltarr(10)
a(*, 0) = [0.68, 3.73, 5.33, 14.0, 49.0, 49.0, 49.0, 4.2, 4.2, 4.2]
a(*, 1) = [0.68, 3.73, 5.33, 14.0, 49.0, 49.0, 49.0, 4.2, 4.2, 18.4]
a(*, 2) = [0.68, 3.73, 5.33, 14.0, 49.0, 49.0, 49.0,5.08,5.08, 18.4]
a(*, 3) = [0.68, 3.73, 5.33, 14.0, 49.0, 49.0, 49.0,3.75,3.75,3.75]
a(*, 4) = [0.68, 3.73, 5.33, 14.0, 49.0, 22.4, 22.4,2.12,2.12,2.12]
a(*, 5) = [0.68, 3.73, 5.33, 14.0, 49.0, 22.4, 46.3, 46.3,6.12,6.12]
a(*, 6) = [0.68, 3.73, 5.33, 14.0, 12.3, 12.3, 12.3, 12.3,6.12,6.12]
a(*, 7) = [0.68, 3.73, 5.33, 14.0, 12.3, 12.3, 10.2, 10.2,4.98,4.98]
a(*, 8) = [0.68, 3.73, 5.33, 10.2, 10.2, 10.2, 10.2, 10.2,4.98,4.98]
a(*, 9) = [0.68, 3.73, 5.33, 10.2, 3.9, 3.9, 2.04, 2.04, 1.1,1.1]
a(*,10) = [0.68, 3.73,.653, 1.04, 1.04, 1.04, 1.04, 1.04,1.04,1.04]
a(*,11) = [0.68, 3.73,.653, .653, .653, .653, 1.04, 1.04,1.04,1.04]
a(*,12) = [0.68,3.73, .653, .653, .653, .653, .653, .653,.653,.653]
a(*,13) = dumarr + 0.6
a(*,14) = dumarr + .37
a(*,15) = dumarr + .053
; Exponent Coefficients
b(*, 0) = [-1.,-1.,-1.595,-.543,-1.572,-1.572,-1.572,-.82,-.82,-.82]
b(*, 1) = [-1.,-1.,-1.595,-.543,-1.572,-1.572,-1.572,-.82,-.82,-1.33]
b(*, 2) = [-1.,-1.,-1.595,-.543,-1.572,-1.572,-1.572, -1., -1.,-1.33]
b(*, 3) = [-1.,-1.,-1.595,-.543,-1.572,-1.572,-1.572, -1., -1., -1.]
b(*, 4) = [-1.,-1.,-1.595,-.543,-1.572, -1.2, -1.2, -1., -1., -1.]
b(*, 5) = [-1.,-1.,-1.595,-.543,-1.572,-1.2,-3.06,-3.06,-1.556,-1.556]
b(*, 6) = [-1.,-1.,-1.595,-.543,-2.09,-2.09,-2.09,-2.09,-1.556,-1.556]
b(*, 7) = [-1.,-1.,-1.595,-.543,-2.09,-2.09,-2.19,-2.19,-1.556,-1.556]
b(*, 8) = [-1.,-1.,-1.595,-2.19,-2.19,-2.19,-2.19,-2.19,-1.556,-1.556]
b(*, 9) = [-1.,-1.,-1.595,-2.19,-2.763,-2.763,-1.31,-1.31,-1.,-1.]
b(*,10) = dumarr -1.
b(*,11) = dumarr -1.
b(*,12) = dumarr -1.
b(*,13) = dumarr -1.
b(*,14) = dumarr -1.
b(*,15) = dumarr -1.
; Exponent Coefficients
c(*, 0) = [0.55,0.21,0.0, 0.0,-.826,-.826,-.826, 4.,4.,4.]
c(*, 1) = [0.55,0.21,0.0, 0.0,-.826,-.826,-.826, 4.,4.,0.0]
c(*, 2) = [0.55,0.21,0.0, 0.0,-.826,-.826,-.826, 3.9,3.9,0.0]
c(*, 3) = [0.55,0.21,0.0, 0.0,-.826,-.826,-.826, 4.2,4.2,4.2]
c(*, 4) = [0.55,0.21,0.0, 0.0,-.826, 0.0, 0.0, 5.6,5.6,5.6]
c(*, 5) = [0.55,0.21,0.0, 0.0,-.826, 0.0, 0.0, 0.0,0.0,0.0]
c(*, 6) = [0.55,0.21,0.0, 0.0,-.208,-.208,-.208,-.208,0.0,0.0]
c(*, 7) = [0.55,0.21,0.0, 0.0,-.208,-.208,-.208,-.208,0.0,0.0]
c(*, 8) = [0.55,0.21,0.0,-.208,-.208,-.208,-.208,-.208,0.0,0.0]
c(*, 9) = [0.55,0.21,0.0,-.208, 0.0, 0.0, 0.0, 0.0,0.58,0.58]
c(*,10) = [0.55,0.21,0.72,0.58, 0.58, 0.58, 0.58, 0.58,0.58,0.58]
c(*,11) = [0.55,0.21,0.72,0.72, 0.72, 0.72, 0.58, 0.58,0.58,0.58]
c(*,12) = [0.55,0.21,0.72,0.72, 0.72, 0.72, 0.72, 0.72,0.72,0.72]
c(*,13) = dumarr + .55
c(*,14) = dumarr + .158
c(*,15) = dumarr + .05
; Exponent Coefficients
d(*, 0) = [-1.,-1.,0.0,0.0,-1.,-1.,-1.,-1.,-1.,-1.]
d(*, 1) = [-1.,-1.,0.0,0.0,-1.,-1.,-1.,-1.,-1.,0.0]
d(*, 2) = [-1.,-1.,0.0,0.0,-1.,-1.,-1.,-1.,-1.,0.0]
d(*, 3) = [-1.,-1.,0.0,0.0,-1.,-1.,-1.,-1.,-1.,-1.]
d(*, 4) = [-1.,-1.,0.0,0.0,-1.,0.0,0.0,-1.,-1.,-1.]
d(*, 5) = [-1.,-1.,0.0,0.0,-1.,0.0,0.0,0.0,0.0,0.0]
d(*, 6) = [-1.,-1.,0.0,0.0,-2.,-2.,-2.,-2.,0.0,0.0]
d(*, 7) = [-1.,-1.,0.0,0.0,-2.,-2.,-2.,-2.,0.0,0.0]
d(*, 8) = [-1.,-1.,0.0,-2.,-2.,-2.,-2.,-2.,0.0,0.0]
d(*, 9) = [-1.,-1.,0.0,-2.,0.0,0.0,0.0,0.0,-1.,-1.]
d(*,10) = dumarr - 1.
d(*,11) = dumarr - 1.
d(*,12) = dumarr - 1.
d(*,13) = dumarr - 1.
d(*,14) = dumarr - 1.
d(*,15) = dumarr - 1.
; Boundaries between temperature bins
tem_lim = [.2,.258,.4,.585,1.,1.5,3.,4.5,8.]
; Boundaries between wavelength bins
wav_lim = [1.4,4.6,6.1,9.1,14.2,16.8,18.6,22.5,25.3,31.6, $
51.9,57.0,89.8,227.9,911.9]
;
;
;
;#####################################################################
; Calculate Free-Free Gaunt factor : gaunt_ff
; (good for 1.e4 < te_6*1.e6 < 1.e9 K; 1 < wave < 1000 Ang)
;#####################################################################
;
low = where (te_6 le 1.,lcount) ; low temperature
high = where (te_6 gt 1.,hcount) ; high temperature
;
if(lcount gt 0) then begin
dummy=mkdarr(wave,te_6(low))
te_dummy=dummy(1,*)
wave_dummy=1. > dummy(0,*) < 1000.
gaunt_ff(low,0:nwave-1)=0.29*wave_dummy^(0.48*(wave_dummy^(-0.08)))$
* te_dummy^(0.133*alog10(wave_dummy)-0.2)
endif
;
if(hcount gt 0) then begin
dummy=mkdarr(wave,te_6(high))
te_dummy=dummy(1,*)
wave_dummy=dummy(0,*)
gaunt_ff(high,0:nwave-1)=1.01*wave_dummy^(0.355*(wave_dummy^(-0.06)))$
* (te_dummy/100.)^(0.3*(wave_dummy^(-0.066)))
endif
;
;#####################################################################
; Calculate 2-photon Gaunt factor
;#####################################################################
;
negl_2p = where((wave le 19.) or (wave gt 200.),ncount)
finit_2p = where((wave gt 19.) and (wave le 200.),fcount)
;
; gaunt_2p(0:nte_6-1,negl_2p) = 0.
; gaunt_2p(0:nte_6-1,finit_2p) = 0. ; 19 < wave <= 200
;
if(fcount gt 0) then begin
dum_2p=mkdarr(wave(finit_2p),te_6)
te_2p=dum_2p(1,*)
wv_2p=dum_2p(0,*)
for i=0,5 do begin
alpha = 106. / (lamda(i)) * (te_2p^(-0.94))
; print,'i=',i,' alpha=',alpha
gaunt_2p(0:nte_6-1,finit_2p)=gaunt_2p(0:nte_6-1,finit_2p) $
+ (wv_2p gt lamda(i)) * CC(i) $
* ((lamda(i)/wv_2p)^alpha) $
* sqrt(abs(cos(!pi*((lamda(i)/wv_2p)-0.5))))$
* ((temp(i) / te_2p)^0.45) $
* 10.^((-DD(i)*(alog10(te_2p/temp(i)))^2)>(-37)) ; (31-jul-93 JRL)
; print,'gaunt_2p=',gaunt_2p
endfor
endif
;
;#####################################################################
; Calculate Free-Bound Gaunt factor
;#####################################################################
;
low0=where(te_6 lt 0.01,lcount0)
low1=where(((te_6 ge 0.01) and (te_6 lt 0.1)),lcount1)
high=where(te_6 ge 0.1,hcount)
; if(lcount0 gt 0) then gaunt_fb(low0,0:nwave-1) = 0.
;
if(lcount1 gt 0) then begin
;-- Low temperature case --- : te_6 < 0.1 M K
il = indd( tem_liml,Te_6(low1)) ; return il
jl = indd( wav_liml,wave ) ; return jl
;
; print,'il,jl -->',il,jl
;
ijl=mkdarr(jl,il)
; print,'ijl -->',ijl
il=ijl(1,*)
jl=ijl(0,*)
;
tel=reform(rebin(te_6(low1),n_elements(te_6(low1)),$
n_elements(wave)),1, $
n_elements(te_6(low1))*n_elements(wave))
;
;
; print,'tel ---> ',tel
;
gaunt_fb(low1,0:nwave-1) = al(il,jl)*(tel^bl(il,jl)) $
* exp( cl(il,jl)*(tel^dl(il,jl)))
endif
;
;
if(hcount gt 0) then begin
;-- High temperature case --- : Te_6 >= 0.1 M K
i = indd( tem_lim,Te_6(high)) ; return i
j = indd( wav_lim,wave ) ; return j
;
; print,'i,j -->',i,j
;
ij=mkdarr(j,i)
; print,'ij -->',ij
i=ij(1,*)
j=ij(0,*)
;
teh=reform(rebin(te_6(high),n_elements(te_6(high)),$
n_elements(wave)),1, $
n_elements(te_6(high))*n_elements(wave))
;
; print,'teh ---> ',teh
;
gaunt_fb(high,0:nwave-1) = a(i,j)*(teh^b(i,j)) $
* exp( c(i,j)*(teh^d(i,j)))
endif
return, gaunt_ff + gaunt_2p + gaunt_fb
end