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