SUBROUTINE ELDENS(ID,T,AN,ANE,ENRG,ENTT,WM,IPRI) C =============================================== C C Evaluation of the electron density and the total hydrogen C number density for a given total particle number density C and temperature; C by solving the set of Saha equations, charge conservation and C particle conservation equations (by a Newton-Raphson method) C C Input parameters: C T - temperature C AN - total particle number density C C Output: C ANE - electron density C ANP - proton number density C AHTOT - total hydrogen number density C AHMOL - relativer number of hydrogen molecules with respect to the C total number of hydrogens C ENERG - part of the internal energy: excitation and ionization C INCLUDE 'IMPLIC.FOR' INCLUDE 'BASICS.FOR' INCLUDE 'MODELQ.FOR' INCLUDE 'ATOMIC.FOR' DIMENSION R(3,3),S(3),P(3) common/terden/rhoter,anta,entrp common/eospar/anmol(600,mdepth), * anato(100,mdepth), * anion(100,mdepth) c if(ioptab.lt.-1) return C if(anerel.le.0.) then anerel=0.5 if(t.lt.9000.) anerel=0.4 if(t.lt.8000.) anerel=0.1 if(t.lt.7000.) anerel=0.01 if(t.lt.6000.) anerel=0.001 if(t.lt.5500.) anerel=0.0001 if(t.lt.5000.) anerel=1.e-5 if(t.lt.4000.) anerel=1.e-6 end if c if(ifmol.gt.0.and.t.lt.tmolim) then aein=an*anerel call moleq(id,t,an,aein,ane,enrg,entt,wm,ipri) anerel=ane/an return end if c QMI=0. Q2=0. QP=0. Q=0. DQN=0. TK=BOLK*T THET=5.0404D3/T anta=an C C Coefficients entering ionization (dissociation) balance of: C atomic hydrogen - QH; C negative hydrogen ion - QM (considered only if IHM>0); C hydrogen molecule - QP (considered only if IH2>0); C ion of hydrogen molecule - Q2 (considered only if IH2P>0). C IF(T.LE.9000.) THEN QMI=1.0353D-16/T/SQRT(T)*EXP(8762.9/T) QP=TK*EXP((-11.206998+THET*(2.7942767+THET* * (0.079196803-0.024790744*THET)))*2.30258509299405) call mpartf(1,1,0,t,uh,duh) uh=max(uh,two) call mpartf(0,0,2,t,uh2,duh2) q2=1.47e-20/(t*sqrt(t))*uh2/uh/uh*exp(51951.8/t) END IF tkln15=log(bolk*t)*1.5 QH0=EXP((15.38287+1.5*LOG10(T)-13.595*THET)*2.30258509299405)*two C ANE=AN*ANEREL IT=0 C C Basic Newton-Raphson loop - solution of the non-linear set C for the unknown vector P, consistiong of AH, ANH (neutral C hydrogen number density) and ANE. C 10 IT=IT+1 C C procedure STATE determines Q (and DQN) - the total charge (and its C derivative wrt temperature) due to ionization of all atoms which C are considered (both explicit and non-explicit), by solving the set C of Saha equations for the current values of T and ANE C CALL STATE(1,ID,T,ANE) C C Auxiliary parameters for evaluating the elements of matrix of C linearized equations. C Note that complexity of the matrix depends on whether the hydrogen C molecule is taken into account C Treatment of hydrogen ionization-dissociation is based on C Mihalas, in Methods in Comput. Phys. 7, p.10 (1967) C IF(IATREF.EQ.IATH.or.ioptab.ge.-1) THEN qh=qh0/pfhyd G2=QH/ANE G3=0. G4=0. G5=0. D=0. E=0. G3=QMI*ANE A=UN+G2+G3 D=G2-G3 IF(IT.GT.1) GO TO 60 IF(T.GT.9000.) THEN F1=UN/A FE=D/A+Q AH=ANE/FE ANH=AH*F1 else if(t.gt.4000.) then E=G2*QP/Q2 B=TWO*(UN+E) GG=ANE*Q2 C1=B*(GG*B+A*D)-E*A*A C2=A*(TWO*E+B*Q)-D*B C3=-E-B*Q F1=(SQRT(C2*C2-4.*C1*C3)-C2)*HALF/C1 FE=F1*D+E*(UN-A*F1)/B+Q AH=ANE/FE ANH=AH*F1 else c1=q2*(two*ytot(id)-un) c2=ytot(id) c3=-an anh=(sqrt(c2*c2-4.*c1*c3)-c2)*half/c1 ah=anh*(un+two*anh*q2) c1=un+qmi*anh c2=-q*ah c3=-qh*anh ane=(sqrt(c2*c2-4.*c1*c3)-c2)*half/c1 end if 60 AE=ANH/ANE GG=AE*QP E=ANH*Q2 B=ANH*QMI C C Matrix of the linearized system R, and the rhs vector S C if(ifmol.eq.0.or.t.ge.tmolim) then R(1,1)=YTOT(ID) R(1,2)=0. R(1,3)=UN S(1)=AN-ANE-YTOT(ID)*AH else R(1,1)=YTOT(ID)-UN R(1,2)=A+E+GG R(1,3)=UN S(1)=AN-ANE-ANH*(A+E+GG)-(YTOT(ID)-UN)*AH end if c R(2,1)=-Q R(2,2)=-D-TWO*GG R(2,3)=UN+B+AE*(G2+GG)-DQN*AH R(3,1)=-UN R(3,2)=A+4.*(E+GG) R(3,3)=B-AE*(G2+TWO*GG) S(2)=ANH*(D+GG)+Q*AH-ANE S(3)=AH-ANH*(A+TWO*(E+GG)) C C Solution of the linearized equations for the correction vector P C CALL LINEQS(R,S,P,3,3) C C New values of AH, ANH, and ANE C AH=AH+P(1) ANH=ANH+P(2) DELNE=P(3) ANE=ANE+DELNE C C hydrogen is not the reference atom C ELSE C C Matrix of the linearized system R, and the rhs vector S C IF(IT.EQ.1) THEN ANE=AN*HALF AH=ANE/YTOT(ID) END IF R(1,1)=YTOT(ID) R(1,2)=UN R(2,1)=-Q-QREF R(2,2)=UN-(DQN+DQNR)*AH S(1)=AN-ANE-YTOT(ID)*AH S(2)=(Q+QREF)*AH-ANE C C Solution of the linearized equations for the correction vector P C CALL LINEQS(R,S,P,2,3) AH=AH+P(1) DELNE=P(2) ANE=ANE+DELNE END IF C C Convergence criterion C IF(ANE.LE.0.) ANE=1.D-6*AN IF(ABS(DELNE/ANE).GT.1.D-3.AND.IT.LE.10) GO TO 10 C C ANEREL is the exact ratio betwen electron density and total C particle density, which is going to be used in the subseguent C call of ELDENS C ANEREL=ANE/AN AHTOT=AH AHMOL=ANH*ANH*Q2 ANP=ANH/ANE*QH ANHM=ANH*ANE*QMI RHOTER=WMY(ID)*AH*HMASS if(ipri.gt.0) then dens(id)=rhoter elec(id)=ane wmm(id)=dens(id)/(an-ane) end if C c internal energy and entropy c call entene(t,ah,anh,anp,ane,energ,entrop) ener=energ entr=entrop c if(id.eq.1) write(6,602) id,t,an,ener,entr c c energy and entropy of H_2 c if(t.lt.9000..and.ahmol.gt.0..and.uh2.gt.0.) then ener=ener+(duh2-51951.8/t)*tk*ahmol entr=entr+ahmol*(tkln15-log(ahmol)+log(uh2)+1.0487+ * 103.973)*bolk end if C enrg=ener entt=entr wm=rhoter/an/hmass c if(ifmol.le.0.or.t.ge.tmolim) then if(n0hn.gt.0) then anato(1,id)=popul(n0hn,id) else anato(1,id)=dens(id)/wmm(id)/ytot(id) end if if(iathe.gt.0) then anato(2,id)=popul(n0a(iathe),id) else anato(2,id)=dens(id)/wmm(id)/ytot(id)*abndd(2,id) end if end if c RETURN END