; file: gasiter_nhtot.pro = gas mix iteration to constant NHtot ; init: Sep 16 - Oct 16 2013 Rob Rutten ; last: Mar 26 2015 Rob Rutten Deil pro gasiter_nhtot,nhtot,temp, $ nhtotnew,elpress,pgas,phneutral,phminus,phplus,ph2,ph2plus,rho, $ verbose=verbose,niter=niter,precision=precision ;+ ; NAME: ; gasiter_nhtot ; PURPOSE: ; iterate to find the LTE partitioning (Saha, molecular) equilibria ; for all hydrogen species and the major electron-donor elements ; for solar element-mix gas at a given total hydrogen density amd ; temperature. Outputs partial pressures and gas density. ; CALL: ; gasiter_nhtot,nhtot,temp, $ ; nhtotnew,elpress,pgas,phneutral,phminus,phplus,ph2,ph2plus,rho, $ ; verbose=verbose,niter=niter,precision=precision ; INPUTS: ; nhtot = total hydrogen density, in cgs. May be constant or array. ; temp = temperature, in K. May be constant or array but equal dim if both. ; KEYWORDS: ; verbose = print out iteration steps, final ratios, error messages ; niter = maximum number of iterations (default = 1000) ; precision = desired precision for nhtot, default = 1E-3 ; OUTPUTS: ; nhtotnew = new value after the iteration ; elpress = partial electron pressure ; pgas = total gas pressure ; phneutral = partial neutral hydrogen pressure ; phminus = partial H-minus pressure ; phplus = partial free proton pressure ; ph2 = partial molecular hydrogen pressure ; ph2plus = partial ionized molecular hydrogen pressure ; rho = gas density ; RESTRICTION: ; calling program needs to contain: ; common common_atomdata_wsa,at_weight,at_abu,at_ionerg1,at_ionerg2,at_ionerg3,at_sym ; atomdata_wsa,1 ; initialize this common ; HISTORY: ; Sep 16 2013 RR: written following p 145-146 RTSA 2003. ; Using LTE equilibrium routines labeled _wsa, found at ; https://github.com/aasensio/lte/tree/master/idl = Andres ; Asensio Ramos. Appear to be old Fortran routines of Wittmann ; described in 1974SoPh...35...11W and translated or rewritten in IDL by ; Jorge Sanchez Almeida in the 1990s. ; Sep 27 2013: RJR changed convergence tests to permit array input to ; make the operation much faster than value-by-value unrolled. ; Oct 10 2013: added rho to output ;- ; no-parameter reply if (n_params() lt 10) then begin ; nr required parameters print,' gasiter_nhtot,nhtot,temp, [input]' print,' nhtotnew,elpress,pgas,phneutral,phminus,phplus,ph2,ph2plus, [output]' print,' verbose=verbose,niter=niter,precision=precision' return endif ; check keywords if (not keyword_set(verbose)) then verbose=0 if (not keyword_set(niter)) then niter=1000 if (not keyword_set(precision)) then precision=1E-3 ; set input parameters to equal array size if scalar plus array nnhtot=size(nhtot,/n_elements) ntemp=size(temp,/n_elements) narr=max([nnhtot,ntemp]) if (ntemp ne nnhtot) then begin if (ntemp eq 1) then temp=fltarr(narr)+temp if (nnhtot eq 1) then nhtot=fltarr(narr)+nhtot if (nnhtot gt 1 and ntemp gt 1) then begin print,' #### dim nhtot = ',ntostr(nnhtot),$ ' but dim temp = ',ntostr(ntemp),' gasiter_nhtot.pro ABORTED' return endif endif ; check temp if (min(temp) lt 2000) then $ print,' ### gasiter_nhtot: at min(temp) = ',ntostr(min(temp)),$ ' result likely unreliable!' ; physics constants kcgs=1.380658D-16 ; Boltzmann constant (erg/deg) ; H and He ionization energies for electron density estimation hi_ion=13.598 ; HI hydrogen ionization energy in eV hei_ion=24.590 ; He I heii_ion=54.403 ; He II ; He abundance heabund=8.51E-02 ; initial estimate of hydrogen ionization and electron density ; using schematic lower-limit metal-donor electron contribution abudonor=5E-3 ;RR abundance of fake collective metal donor nel=abudonor*nhtot ; Saha for fake metal saha_mc,temp,nel,6.0,2,1,n1_metal,n0_metal ;RR 6.0 = ionization energy in eV for fake collective metal donor nelectron=abudonor*n1_metal*nhtot ; iterate electron density only solving HI > HII ionization hiparfun=partfunc_nv(temp,1,0) nel=nelectron for isaha=0,1000 do begin saha_mc,temp,nel,hi_ion,hiparfun,1.0,n1_ntot,n0_ntot nelnew=nelectron+n1_ntot*nhtot if (max(abs(nelnew-nel)/nel) lt abudonor) then break nel=nelnew endfor nelectron=nelnew nhneutral=n0_ntot*nhtot elpress=nelectron*kcgs*temp ; refine by iteration using Wittmann-Sanchez-Asensio codes which ; include all hydrogen species plus electrons from HeI>HeII and HeII>HeIII ; plus electrons from ionizing the major neutral and once-ionized metals for iiter=0,niter do begin ; get all hydrogen species gaspress_wsa,elpress,temp,phneutral,phminus,phplus,ph2,ph2plus,pgas,rho nhtotnew=(phneutral+phminus+phplus+2*ph2+2*ph2plus)/(kcgs*temp) if (max(abs(nhtot-nhtotnew)/nhtot) lt precision) then break elpress=nhtot/nhtotnew*elpress endfor ; print results if verbose and only single variable input if (verbose and narr eq 1) then begin print,' input: N_Htot = '+ntostr(nhtot)+' temp = '+ntostr(temp) if (iiter GE niter) then print,' #### gasiter_nhtot: not converged' print,' ==== nr steps initial iteration = '+ntostr(isaha) print,' ==== initial result; nhneutral = '+ntostr(nhneutral)+$ ' nelectron = '+ntostr(nelectron) print,' ==== nr steps full iteration = '+ntostr(iiter) print,' ==== result iteration: N_Htot = '+ntostr(nhtotnew)+$ ' N_Hneutral = '+ntostr(phneutral/(kcgs*temp))+$ ' N_e = '+ntostr(elpress/(kcgs*temp)) print,' ==== H fractions of NH_tot: '+$ ' H-neutral='+ntostr(phneutral/(kcgs*temp)/nhtotnew)+$ ' H-minus='+ntostr(phminus/(kcgs*temp)/nhtotnew)+$ ' H_p='+ntostr(phplus/(kcgs*temp)/nhtotnew)+$ ' H2='+ntostr(2*ph2/(kcgs*temp)/nhtotnew)+$ ' H2+='+ntostr(2*ph2plus/(kcgs*temp)/nhtotnew) print,' ==== N_HI/N_Hmin = '+ntostr(phneutral/phminus)+$ ' N_e/N_p='+ntostr(elpress/phplus) endif ; array input and verbose mode: less print but precision plot if (verbose and narr gt 1) then begin sv,alog10(abs(nhtotnew-nhtot)/nhtot) print, ' gasiter_nhtot: see plot; isaha='+ntostr(isaha)+$ ' iiter='+ntostr(iiter) endif end ; ================= test per IDLwave Hyper-C =========================== ; common with atomic data in the Wittmann-Sanchez-Asensio programs common common_atomdata_wsa,at_weight,at_abu,at_ionerg1,at_ionerg2,at_ionerg3,at_sym atomdata_wsa,1 ; initialize common with these numbers ; choose params nhtot=1.E17 temp=10000. temp=6000. temp=50000. ; get values (NB: answer 1E17, 6000 = 4E7 HI/H^min ratio RTSA p178 below (8.3) gasiter_nhtot,nhtot,temp,$ nhtotnew,elpress,pgas,phneutral,phminus,phplus,ph2,ph2plus,/verbose end