110 lines
2.8 KiB
Fortran
110 lines
2.8 KiB
Fortran
SUBROUTINE TODENS(ID,T,AN,ANE)
|
|
C ==============================
|
|
C
|
|
C determines AN (and ANP, AHTOT, and AHMOL) from T and ANE
|
|
C
|
|
C Input parameters:
|
|
C T - temperature
|
|
C ANE - electron number density
|
|
C
|
|
C Output:
|
|
C AN - total particle density
|
|
C ANP - proton number density
|
|
C AHTOT - total hydrogen number density
|
|
C AHMOL - relative number of hydrogen molecules with respect to the
|
|
C total number of hydrogens
|
|
C
|
|
INCLUDE 'PARAMS.FOR'
|
|
INCLUDE 'MODELP.FOR'
|
|
common/hydmol/anhmi,ahmol
|
|
parameter (un=1.d0,two=2.d0,half=0.5d0)
|
|
C
|
|
QM=0.
|
|
Q2=0.
|
|
QP=0.
|
|
Q=0.
|
|
DQN=0.
|
|
TK=BOLK*T
|
|
THET=5.0404D3/T
|
|
C
|
|
C Coefficients entering ionization (dissociation) balance of:
|
|
C atomic hydrogen - QH;
|
|
C negative hydrogen ion - QM
|
|
C hydrogen molecule - QP
|
|
C ion of hydrogen molecule - Q2
|
|
C
|
|
QM=1.0353D-16/T/SQRT(T)*EXP(8762.9/T)
|
|
QH=EXP((15.38287+1.5*LOG10(T)-13.595*THET)*2.30258509299405)
|
|
c
|
|
if(t.gt.16000.) then
|
|
ih2=0
|
|
ih2p=0
|
|
else
|
|
QP=TK*EXP((-11.206998+THET*(2.7942767+THET*
|
|
* (0.079196803-0.024790744*THET)))*2.30258509299405)
|
|
Q2=TK*EXP((-12.533505+THET*(4.9251644+THET*
|
|
* (-0.056191273+0.0032687661*THET)))*2.30258509299405)
|
|
ih2=1
|
|
end if
|
|
|
|
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(ID,T,ANE,Q)
|
|
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
|
|
G2=QH/ANE
|
|
G3=0.
|
|
G4=0.
|
|
G5=0.
|
|
D=0.
|
|
E=0.
|
|
G3=QM*ANE
|
|
A=UN+G2+G3
|
|
D=G2-G3
|
|
IF(IT.LE.1) THEN
|
|
IF(IH2.EQ.0) THEN
|
|
F1=UN/A
|
|
FE=D/A+Q
|
|
ELSE
|
|
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
|
|
END IF
|
|
AH=ANE/FE
|
|
ANH=AH*F1
|
|
END IF
|
|
AE=ANH/ANE
|
|
GG=AE*QP
|
|
E=ANH*Q2
|
|
B=ANH*QM
|
|
C
|
|
c S(1)=AN-ANE-YTOT(ID)*AH
|
|
c S(2)=ANH*(D+GG)+Q*AH-ANE
|
|
c S(3)=AH-ANH*(A+TWO*(E+GG))
|
|
c
|
|
hhn=A+TWO*(E+GG)
|
|
anh=ane/(d+gg+q*hhn)
|
|
ah=anh*hhn
|
|
an=ane+ytot(id)*ah
|
|
C
|
|
AHTOT=AH
|
|
AHMOL=TWO*ANH*(ANH*Q2+ANH/ANE*QP)/AH
|
|
ANP=ANH/ANE*QH
|
|
RETURN
|
|
END
|