Документ взят из кэша поисковой машины. Адрес оригинального документа : http://hea-www.harvard.edu/PINTofALE/pro/nenh.pro Дата изменения: Mon Apr 14 22:07:56 2003 Дата индексирования: Mon Oct 1 23:42:40 2012 Кодировка: Поисковые слова: mdi

function nenh,logT,abund=abund,eqfile=eqfile,chidir=chidir,$
elem=elem,xelem=xelem,nproton=nproton, _extra=e
;+
;function nenH
; calculate the ratio N(e)/N(H) (or N(e)/N(p) -- see keyword
; NPROTON below) in a plasma of specified temperature
;
;syntax
; r=nenh(logT,abund=abund,eqfile=eqfile,chidir=chidir,elem=elem,$
; xelem=xelem,/nproton)
;
;parameters
; logT [INPUT; required] temperature(s) (log_10(T [K])) at which
; to compute ratio.
;
;keywords
; abund [I/O] abundances of elements. default is Allen abundances.
; eqfile [INPUT] pathname, relative to CHIDIR, of file containing
; ionization equilibrium values, passed w/o comment to RD_IONEQ
; chidir [INPUT] path name to the CHIANTI installation, passed w/o
; comment to RD_IONEQ
; elem [INPUT] return data for only these elements (default: all)
; xelem [INPUT] exclude these elements (overrides ELEM)
; * default: none
; (note that older CHIANTI ion balance files do not have
; Cu, Zn, so originally XELEM was set by default to [29,30])
; nproton [INPUT] if set, returns the ratio N(e)/N(p)
; _extra [JUNK] ignore -- here only to prevent crashing program
;
;subroutines
; GETABUND [RDABUND]
; RD_IONEQ [READ_IONEQ (a CHIANTI routine)]
; SYMB2ZION [LAT2ARAB]
; KILROY
;
;history
; vinay kashyap (Sep97)
; modified to have just one call to RD_IONEQ instead of one for
; each Z, changed default of XELEM to none (VK; Jun02)
;-

; usage
if n_elements(logT) eq 0 then begin
print,'Usage: r=nenh(logT,abund=abund,eqfile=eqfile,chidir=chidir,$'
print,' elem=elem,xelem=xelem,/nproton)'
print,' compute ratio of e/H number density at specified temperatures'
return,0.
endif

; initialize
tlog=[logT(*)] & nt=n_elements(tlog) ;at which temperatures..
atom=[ 'H','He','Li','Be','B', 'C','N','O','F','Ne','Na','Mg','Al','Si',$
'P','S','Cl','Ar','K','Ca','Sc','Ti','V','Cr','Mn','Fe','Co','Ni',$
'Cu','Zn'] ;elements from 1-30
nz=n_elements(atom) ;{H..Zn}
if n_elements(abund) ne nz then abund=getabund('allen')
reH=fltarr(nt) ;OUTPUT

; if ELEM and/or XELEM are set...
if not keyword_set(xelem) then xelem=[-1]
sze=size(elem) & szx=size(xelem) & nsze=n_elements(sze) & nszx=n_elements(szx)
etyp=sze(nsze-2) & xtyp=szx(nszx-2) & elm=intarr(nz)
if etyp eq 0 then elm=indgen(nz)+1 ;default -- all elements
if etyp gt 0 and etyp le 7 then begin
if etyp eq 7 then begin qut=!quiet & !quiet=1 & endif
nelm=sze(nsze-1)
for i=0,nelm-1 do begin
if etyp eq 7 then symb2zion,elem(i),z,ion else z=fix(elem(i))
if z gt 0 then elm([z-1])=z
endfor
if etyp eq 7 then !quiet=qut
endif
if xtyp gt 0 and xtyp le 7 then begin
if xtyp eq 7 then begin qut=!quiet & !quiet=1 & endif
nxlm=szx(nszx-1)
for i=0,nxlm-1 do begin
if xtyp eq 7 then symb2zion,xelem(i),z,ion else z=fix(xelem(i))
if z gt 0 then elm([z-1])=0
endfor
if xtyp eq 7 then !quiet=qut
endif
;
oz=where(elm gt 0,mz) ;select elements

; compute electron number
ieq=rd_ioneq(elm(oz),tlog,eqfile=eqfile,chidir=chidir)
nprot=fltarr(nt)
for iz=0,mz-1 do begin ;{for each element in list
zz=elm(oz(iz)) & nion=zz+1L > 1
kilroy; was here
if mz eq 1 then ioneq=ieq else ioneq=reform(ieq(*,0L:nion-1L,iz),nt,nion)
numelec=findgen(zz+1) ;# of electrons for each ionization stage
for it=0,nt-1 do begin ;{for each temperature
reH(it)=reH(it)+abund(zz-1)*total(numelec*ioneq(it,*))
if zz eq 1 then nprot(it)=abund[oz[iz]]*ioneq[it,1]
endfor ;IT=0,NT-1}
endfor ;IZ=0,MZ-1}

; If the proton density were required

; this part because of CHIANTI routine PROTON_DENS()
; -- note that the default operation of NENH() is equivalent
; to PROTON_DENS(/hydrogen) and the default operation of
; PROTON_DENS() is equivalent to NENH(/nproton))
;
; example:
; r=nenh(!logt,/nproton) & rc=proton_dens(!logT)
; plot,!logT,rc,/ynoz & oplot,!logT,1./r,col=2
; example:
; r=nenh(!logt) & rc=proton_dens(!logT,/hydrogen)
; plot,!logT,1./rc & oplot,!logT,r,col=2

oH=where(nprot gt 0,moH)
if moH eq 0 then begin
; was H included in XELEM? tch tch.
ieq=rd_ioneq(1,tlog,eqfile=eqfile,chidir=chidir)
nprot(*)=abund(0)*ioneq(*,1)
endif
oH=where(nprot gt 0,moH)
if moH eq 0 then message,'BUG?!'
if keyword_set(nproton) then begin
tmp=0.*reH+1.
tmp[oH]=reH[oH]/nprot[oH] & reH=tmp
endif

return,reH
end