PROGRAM NEOT c this program computes the equilibirum ionization balance for a c gas with Ptot and T implicit none include 'neoT.h' include 'const.dek' integer i,j,k,NTiter external func double precision a,b,tol,ne,zbrent double precision sumj,sumk,sumkj,m_amu,donate double precision Tmin,Tmax,dT double precision rhototN,rhototmu,rhoN,rhoe,beta double precision Ptot,mutot,muN,mue,XkoAk(Nelem) double precision nN,nk(Nelem),njk(Jmax,Nelem),fnenjk(Jmax,Nelem) include 'neoT.com' c X,Y,Z = mass fractions c Jk(k) = number of ionization stages for species k c Ak(k) = atomic weight for species k c alpha(k) = abundance fraction for species k c elem(k) = name of element of species k c species(j,k) = name of ionization species j of species k c set the rootsolving parameters for zbrent a = 1.0d5 ! lower limit to ne for rootsolving zbrent call b = 1.0d25 ! upper limit to ne for rootsolving zbrent call tol = 1.0d-5 ! tolerance for ne value from zbrent c set constants, compute combined constants m_amu = 1.66053878d-24 Cphi = 2.0d0*(2.0d0*pi*me*kerg/(h*h))**(1.50d0) c open the data files and populate arrays... c file='XYZ.tab' contains X,Y, and Z c file='Ak.tab' contains elem(k), Jk(k), and Ak(k) c file='UofT.tab' contains Ttab(i) and Utab(j,k,i) CALL readdata c compute the abundance fractions alpha(1) = (1.0d0+(Ak(1)/Ak(2))*(Y/X)+(Ak(1)/Ak(3))*(Z/X))**(-1) alpha(2) = alpha(1)*(Ak(1)/Ak(2))*(Y/X) alpha(3) = alpha(1)*(Ak(1)/Ak(3))*(Z/X) C WRITE(6,*) alpha(1),alpha(2),alpha(3) c compute the nuclear mean molecular weight XkoAk(1) = X/Ak(1) XkoAk(2) = Y/Ak(2) XkoAk(3) = Z/Ak(3) sumk = 0.0d0 DO 09 k=1,Nelem sumk = sumk + XkoAk(k) 09 CONTINUE muN = 1.0d0/sumk c set the boundary conditions Ptot = 100.0d0 Tmin = 1500.0d0 Tmax = 25000.0d0 dT = 250.0d0 NTiter = nint((Tmax-Tmin)/dT) + 1 c open the output data files containing the results c file='neoT.out' contains total, nuclear, e- densities, mus, rhos, beta c file='neoT.hydrogen' contains hydrogen results c file='neoT.helium' contains hydrogen results c file='neoT.calcium' contains hydrogen results OPEN(unit=4,file='neoT.out',status='unknown') OPEN(unit=1,file='neoT.hydrogen',status='unknown') OPEN(unit=2,file='neoT.helium',status='unknown') OPEN(unit=3,file='neoT.calcium',status='unknown') c stamp the pressure, write the column headers for all output files WRITE(4,600) 'Ptot',Ptot CALL headers c *** MAIN TEMPERATURE LOOP *** , index i c step through temperature starting at Tmin in steps of dT to Tmax DO 01 i=1,NTiter T = Tmin + float(i-1)*dT c compute ntot, the boundary condition on the gas for this T c populate the partition functions U(j,k)(T) (requires interpolation) c populate the Saha Phi(j,k)(T) ntot = Ptot/(kerg*T) CALL pop_UofT CALL pop_PhiofT c conserve particle and charge density; root solve for ne at this T ne = zbrent(func,a,b,tol) c apply particle conservation, and obtain the respective densities c use appropriate sums to compute the electron mean molecular weight c compute the total mean molecular weight nN = ntot - ne rhoe = me*ne beta = ne/ntot rhoN = 0.0d0 ! initialize for summation of nk*Ak*m_amu sumkj = 0.0d0 ! initialize for summation of (Xk/Ak)*sumj DO 11 k=1,Nelem nk(k) = alpha(k)*nN sumj = 0.0d0 ! initialize for summation of j*fjk DO 13 j=1,Jk(k) njk(j,k) = f(j,k)*nk(k) IF ((k.eq.1).AND.(j.eq.1)) donate = -1.0d0 !H- removes an electron IF ((k.eq.1).AND.(j.eq.2)) donate = 0.0d0 !H0 does not donate IF ((k.eq.1).AND.(j.eq.Jk(k))) donate = 1.0d0 !H+ donates 1 IF (k.gt.1) donate = float(j-1) !general case fnenjk(j,k) = donate*njk(j,k)/ne sumj = sumj + donate*f(j,k) 13 CONTINUE rhoN = rhoN + nk(k)*Ak(k)*m_amu sumkj = sumkj + XkoAk(k)*sumj 11 CONTINUE mue = 1.0d0/sumkj mutot = 1.0d0/((1.0d0/muN)+(1.0d0/mue)) c compute the mass densities: c (1) from sum of rhoN and rhoe c (2) from mutot, the total mean molecular weight rhototN = rhoN + rhoe rhototmu = ptot*mutot*m_amu/(kerg*T) c communicate results c report to the 'neoT.out' file c report element results c k=1 'neoT.hydrogen' c k=2 'neoT.helium' c k=3 'neoT.calcium' WRITE(4,602) T,ntot,nN,ne,muN,mue,mutot,rhoN,rhoe,rhototN, & rhototmu,beta DO 25 k=1,Nelem WRITE(k,602) T,nk(k),(njk(j,k),j=1,Jk(k)), & (f(j,k),j=1,Jk(k)),(fnenjk(j,k),j=1,Jk(k)), & (Phi(j,k),j=1,Jk(k)-1),(log10(U(j,k)),j=1,Jk(k)) 25 CONTINUE c increment to next temperature 01 CONTINUE c close up shop CLOSE(unit=1) CLOSE(unit=2) CLOSE(unit=3) CLOSE(unit=4) STOP c formats 600 FORMAT(1x,a9,1x,1pe12.5) 602 FORMAT(2x,26(1x,1pe11.4)) END c include numerical recipe codes include 'splinecodes.f' include 'zbrent.f' c begin subroutines c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: DOUBLE PRECISION FUNCTION func(ne) c func(ne) is the application of particle and charge density c conservation solved at the current T c this function is evaluated iteratively via routine zbrent until c the value of func(ne)=0 within the specified tolerace, tol. zbrent c sends down the guess ne from which we evaluate func(ne). when c ne solves func(ne)=0, the routine provides the equilibrium conditions c of the gas, in particular the ionization fractions f(j,k) c in this routine, we compute the ionization fractions, f(j,k) c current (guess) ne and the current T, the computations requires c the appropriate products and summations for application of the c recursive Saha equation c it is required that the Phi(j,k) are supplied as global variables c the temperature T is passed globally to all routines c all global (non-local) variables are included in "neoT.com" c.............................................................................. implicit none include 'neoT.h' integer j,k double precision sumj,sumkj,donate double precision ne,P(Jmax,Nelem),S(Nelem) include 'neoT.com' c to conserve charge density for the current ne we must compute c P(j,k), S(k), and f(j,k) for each ionization stage of each element DO 61 k=1,Nelem ! loop over elements k c compute the P(j,k); since j=1 is neutral using the array indices, c note the run of the index j is reduced by j-1 so that j-1=0 for c the neutral species; the P(j,k) are the products of Phi(j-1,k)/ne, c where P(1,k)=1 by definition P(1,k) = 1.0d0 DO 15 j=2,Jk(k) P(j,k) = P(j-1,k)*Phi(j-1,k)/ne 15 CONTINUE c compute the S(k), which is the sum of the P(j,k) for species k S(k) = 0.0d0 DO 13 j=1,Jk(k) S(k) = S(k) + P(j,k) 13 CONTINUE c compute the f(j,k), which is the ratio P(j,k)/S(k) DO 11 j=1,Jk(k) f(j,k) = P(j,k)/S(k) 11 CONTINUE 61 CONTINUE ! next element k c armed with the ionization fractions, we now compute the sum of c electron density donaed by ionized nuclear donors; again, since j=1 c index is for neutrals, the number of donated electrons is j-1 sumj = 0.0d0 sumkj = 0.0d0 DO 41 k=1,Nelem DO 43 j=1,Jk(k) IF ((k.eq.1).AND.(j.eq.1)) donate = -1.0d0 !H- removes an electron IF ((k.eq.1).AND.(j.eq.2)) donate = 0.0d0 !H0 does not donate IF ((k.eq.1).AND.(j.eq.Jk(k))) donate = 1.0d0 !H+ donates 1 IF (k.gt.1) donate = float(j-1) !general case sumj = sumj + donate*f(j,k) 43 CONTINUE sumkj = sumkj + alpha(k)*sumj 41 CONTINUE c conserve charge density; zbrent root solves this function func = (ntot-ne)*sumkj - ne RETURN END c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: SUBROUTINE pop_UofT c populate the partition function arrays U(j,k)(T) c the temperature T is passed globally to all routines c all global (non-local) variables are included in "neoT.com" c.............................................................................. implicit none include 'neoT.h' integer i,j,k,n double precision lg10U,theta double precision xT(NT),yU(NT),dy0,dy2(NT) include 'neoT.com' c the data tables are a function of Theta, convert T theta = 5040.0d0/T c if we are at low T, above the maximum theta=2.0 in the Ttab(i) table, c then use the statistical weight for the partition function. the c statistical weight is stored as the last entry i=NT c if theta is bounded within the Ttab(i) table, then interpolate to c obtain the partition functions for this T c interpolation subroutines spline and splint require that we pass c the Ttab and Utab data via the argument list, since arguments passed c to subroutines cannot be global variables, we must store Ttab and Utab c as local variables prior to the call to the subroutines' we must also c pass the number of entries in the tables as a local variable n = NT-1 ! local variable of number of entires in the tables dy0 = 1.0d30 ! use natural spline derivitives at the egdes c for each element k, loop over the ionization species j DO 11 k=1,Nelem DO 09 j=1,Jk(k) c are we off the Ttab table? if so, then use the statistical weight IF (theta.gt.Ttab(n)) then U(j,k) = 10.0d0**Utab(j,k,n+1) GOTO 09 ! skip to next ionization stage for this k END IF c if we are here, then we are on the Ttab(i) table, we must c store global tables as local arrays for spline/splint calls DO 07 i=1,n xT(i) = Ttab(i) yU(i) = Utab(j,k,i) 07 CONTINUE c compute the cubic spline second derivitives for this j,k; they c are returned as the array dy2(i) CALL spline(xT,yU,n,dy0,dy0,dy2) c interpolate on the tables at the current theta, c returns lg10U=log10[U(j,k)(T)] CALL splint(xT,yU,dy2,n,theta,lg10U) c store the result in base 10 U(j,k) = 10.0d0**lg10U 09 CONTINUE ! next j 11 CONTINUE ! next k RETURN END c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: SUBROUTINE pop_PhiofT c compute the quantity Phi(j,k)(T) for all ionization stages of c all elements for the current T c the temperature T is passed globally to all routines c prior to calling this routine, we have populated the partition c functions U(j,k)(T) c all global (non-local) variables are included in "neoT.com" c.............................................................................. implicit none include 'neoT.h' include 'const.dek' integer j,k double precision echikT,Urat,T32 include 'neoT.com' c the loop over ionization stages only goes to Jk(k)-1 because the c Phi(j,k) sorresponds to the lower ionization stage j DO 11 k=1,Nelem DO 09 j=1,Jk(k)-1 echikT = exp(-chi(j,k)/(kerg*T)) ! Boltzmann factor T32 = T**(1.50d0) ! T to the 3/2 power Urat = U(j+1,k)/U(j,k) ! ratio of partition functions Phi(j,k) = Cphi*T32*Urat*echikT ! Saha Equation 09 CONTINUE 11 CONTINUE RETURN END c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: SUBROUTINE readdata c read in the data files for this gas c all global (non-local) variables are included in "neoT.com" c.............................................................................. implicit none include 'neoT.h' include 'const.dek' integer i,j,k character*80 dumstr include 'neoT.com' c *** X,Y,Z c X = hydrogen mass fraction c Y = helium mass fraction c Z = metals mass fraction OPEN(unit=1,file='XYZ.tab',status='old') READ(1,*) dumstr,X READ(1,*) dumstr,Y READ(1,*) dumstr,Z CLOSE(unit=1) c *** Jk(k) and Ak(k) data c Jk(k) = number of ionization stages for species k c Ak(k) = automic weight for speiceis k OPEN(unit=1,file='Ak.tab',status='old') DO 21 k=1,Nelem READ(1,*,end=25) elem(k),Jk(k),Ak(k) 21 CONTINUE 25 CLOSE(unit=1) c *** chi(j,k) data c chi(j,k) = ionization potential of ionization stage j species k c convert from [eV] to [ergs] on the fly OPEN(unit=1,file='ChiI.tab',status='old') DO 11 k=1,Nelem DO 13 j=1,Jk(k)-1 READ(1,*,end=15) dumstr,chi(j,k) chi(j,k) = ev2erg*chi(j,k) 13 CONTINUE 11 CONTINUE 15 CLOSE(unit=1) c Ttab(j,k) and Utab(j,k,i) tabular data c read in the temperature partition function table, Ttab(i) and c the partition function data table grid Utab(j,k,i) c for interpolation, column 1 contains the ionization species name c there are 10 columns of data plus an 11th column that is the c statistical weight for low T (theta > 2.0) OPEN(unit=1,file='UofT.tab',status='old') READ(1,*) dumstr,(Ttab(i),i=1,NT) ! actually Theta DO 03 k=1,Nelem DO 01 j=1,Jk(k) READ(1,*,end=05) species(j,k),(Utab(j,k,i),i=1,NT) 01 CONTINUE 03 CONTINUE 05 CLOSE(unit=1) RETURN END c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: SUBROUTINE headers c write the headers to the output files c.............................................................................. implicit none include 'neoT.h' include 'const.dek' integer j,k include 'neoT.com' c the headers for unit=4, file='neoT.out' WRITE(4,602) 'T','ntot','nN','ne','muN','mue','mutot','rhoN', & 'rhoe','rhototN','rhototmu','beta' c the headers for unit=1, file='neoT.hydrogen' c the headers for unit=2, file='neoT.helium' c the headers for unit=3, file='neoT.calcium' DO 11 k=1,Nelem WRITE(k,604) 'T','n',k,('n',j-1,k,j=1,Jk(k)), & ('f',j-1,k,j=1,Jk(k)),('fne',j-1,k,j=1,Jk(k)), & ('Phi',j-1,k,j=1,Jk(k)-1),('U',j-1,k,j=1,Jk(k)) 11 CONTINUE RETURN c the formats for the headers, 602 assumes Jk(j)=2 c whereas 603 and 604 assume Jk(j)=3 602 FORMAT(1x,t5,a1,t17,a4,t29,a2,t41,a2,t53,a3,t65,a3,t77,a5, & t89,a4,t101,a4,t113,a7,t125,a8,t137,a4) 604 FORMAT(1x,4x,a1,10x,a1,i1,1x,3(9x,a1,i1,i1), & 3(9x,a1,i1,i1),2x,3(7x,a3,i1,i1), & 3(7x,a3,i1,i1),3(9x,a1,i1,i1)) END