SUBROUTINE EOS (Pin,t, X, Y, Z, tcut,q,cp,delad, 1 RHO, chirho, chit, u, cv, eta,fmu) C C FINDS RHO GIVEN PRESSURE AND TEMPERATURE. C HYDROGEN IS ASSUMED TO HAVE C ONE STATE WHEREAS HELIUM HAS TWO STATES (TWO ELECTRONS C IN GROUND.) THE ROUTINE PASSES THE ELCTRON PRESSURE C (TO BE USED IN THE OPACITY) AND DELAD. THE EQUATION C OF STATE IS THAT OF AN IDEAL MONATOMIC GAS PLUS C RADIATION. IONIZATION EFFECTS ARE INCLUDED IN C THE CALCULATION OF THE THERMODYNAMIC DERIVATIVES. C *********************************************** C THE PRIMARY INARDS OF THIS ROUTINE ARE DUE TO C W. D. PESNELL ... WHOM WE THANK. c *********************************************** c IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION YGO(4) PARAMETER (ACCUR=1.0D-6) pwant=pin VGUESS=1.d0/(Pin/T/8.428d7) DO 1 ITRY=1,50 CALL GETEOS (T,VGUESS,tcut,X,Y,Z,P,PE,PV,PT,CV,U,OOMU) PV=VGUESS*PV/P FP=DLOG(P)-dlog(PWANT) DELV=-FP/PV IF (DABS(DELV) .GT. 0.30) DELV=0.30*DELV/DABS(DELV) VGUESS=DLOG(VGUESS)+DELV VGUESS=DEXP(VGUESS) c IF (DABS(DELV) .LT. ACCUR) THEN CALL GETEOS (T,VGUESS,tcut,X,Y,Z,P,PE,PV,PT,CV,U,OOMU) PWANT=pin FP=P-PWANT DELV=-FP/PV VGUESS=VGUESS+DELV RHO=1.0D0/VGUESS c These are chit, chirho, and Gamma3-1. CHIT=PT*T/PWANT CHIRHO=-PV/(PWANT*RHO) G3=PWANT*CHIT/(RHO*T*CV) DELAD=CHIT+CHIRHO/G3 DELAD=1.0D0/DELAD c other quantities used in evolution c Q=-chirho/chit Q=chirho/chit g1=chirho+chit*g3 cp=cv*g1/chirho fmu=1.d0/oomu c guess at degeneracy parameter (assuming only mild degeneracy fne=pe/(1.38d-16*t) tl=dlog10(t) eta=-36.113d0+dlog(fNe)-3.4539d0*TL c RETURN c ENDIF 1 CONTINUE WRITE (6,1000) T, DEXP(PWANT),X,Y 1000 FORMAT (' EOS DID NOT CONVERGE. UGH!. T=', * 1PD11.3,' P=', D11.3/' X= ',0pf7.3,' Y=',f7.3) WRITE (6,1001) 1001 FORMAT (' CALCULATION STOPPED') STOP 'SORRY' END c *********************************************** SUBROUTINE GETEOS (TIN,VIN,tcut,X,Y,Z,P,PE,PV,PT,ET,E,TOTLN) C *********************************************** c This routine and the other EOS routines are due c to W. Dean Pesnell. C C EQUATION OF STATE C ARGUMENTS...TIN (DEGREES K), VIN=1/RHO (CM**3/GM) C METALS... C NA,AL ALWAYS IONIZED C MG,SI,FE INCLUDED AS SINGLE ELEMENT C ALL OTHERS IGNORED c c set tcut to be log of temperature above which full c ionization is enforced C C TABLE OF RETURNED QUANTITIES. C C FUNCTION NAME DERIVATIVE WITH RESPECT TO C TEMP. SP. VOL. C------------------------------------------------------- C PRESSURE I P I PT I PV I C INT. ENERGY I E I ET I EV I C ELEC. PRS. I PE I PET I PEV I C MOLAR DENSITYI BP I BPT I BPV I C------------------------------------------------------- C C C X = HYDROGEN MASS FRACTION C Y = HELIUM MASS FRACTION C Z = METALLIC MASS FRACTION C PARGRM = MEAN MOLECULAR MOLAR DENSITY WITHOUT C ELECTRONS C C R = GAS CONSTANT 8.31434E7 C A = STEFAN-BOLTZMAN CONSTANT 7.56471E-15 C BK = BOLTZMAN S CONSTANT 8.6170837E-5 C AVAGD = AVAGADRO S NUMBER 6.02217E23 C AD3 = A/3 C IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (ZERO = 0.D0,ONE = 1.D0, TWO=2.D0, * THRE=3.D0, FOR=4.D0, TEN=10.D0, AHF=0.5D0, * QRT=0.25D0 ) PARAMETER ( R = 8.31434D7, A = 7.56471D-15, * BK = 8.6170837D-5, AVAGD = 6.02217D23, * AD3 = A/3.D0 ) DATA T3OUT,T4OUT/1.665795163D-25,3.802592017D-28/ DATA T2OUT,T5OUT/ 5.347896D-35,6.614536D-34/ C C IONIZATION POTENTIALS FOR HYDROGEN AND HELIUM C DATA XH,XHE,XHE2/13.595D0,24.581D0,54.403D0/ DATA C1,C2,C3/4.0092926D-9,1.00797D0,4.0026D0/ DATA XM,CM,ZPZP/7.9D0,0.7D0,0.12014D0/ DATA PREC / 1.D-10 / DATA ONHLF/1.5D0/ tl=dlog10(tin) V = VIN T = TIN IF( V .LE. ZERO ) GO TO 11 IF( T .LE. ZERO ) GO TO 10 FRE = ZERO ENT = ZERO PARGRM = X/C2 + Y/C3 RMUC = ONE/PARGRM RT = R*T TT4 = T**4 TK = ONE/(T*BK) SQT = DSQRT(T) C C1=ORIGINAL(C1(0.33334622))/R T1 = V*SQT**3*C1 T2 = T2OUT IF( T .GT. 2.D3 ) T2 = DEXP(-XH*TK) T3 = T3OUT IF( T .GT. 5.D3 ) T3 = DEXP(-XHE*TK) T4 = T4OUT IF( T .GT. 1.D4 ) T4 = DEXP(-XHE2*TK) T5 = T5OUT IF( T .GT. 1.2D3 ) T5 = DEXP(-XM*TK) c c check for tcut and therefore enforce complete ionization c smoothly interpolate between full ionization at tcut+0.2 c and unionized at tcut and below. Force full ionization c by setting t1 equal to about 1000. c t1max=10000d0 dlt1max=4.d0 c if (tl.lt.tcut+0.2d0.and.tl.gt.tcut) then factor=(tl-tcut)/0.2d0 t1=10.d0**(dlog10(t1)+factor*dlt1max) elseif (tl.ge.tcut+0.2d0) then t1=t1max endif c D = T1*T2 B = FOR*T1*T3 C = B*T1*T4 DD = TWO*CM*T1*T5 ZNA = Z*2.48D-3/24.969D0 ZMG = Z*ZPZP/45.807D0 C C CONVERGE ON ELECTRON DENSITY USING THE SAHA C EQUATION. C C GES IS THE MOLAR DENSITY OF ELECTRONS. C GES = (X+Y*AHF)/(ONE+Y/(FOR*C)) IF( GES .LT. X ) GES = AHF*(DSQRT(D*(D+FOR*X))-D) IF( GES .LT. 1.D-6*Z ) GES = 1.D-6*Z XC2 = X/C2 YC3 = Y/C3 C C NEWTON METHOD FOR ELECTRON DENSITY. C DO 1 I=1,25 T2 = C/GES+GES+B GEP = XC2*D/(GES+D)+YC3*(B+TWO*C/GES)/T2 @ + ZMG*DD/(GES+DD) + ZNA T1 = ONE+XC2*D/(D+GES)**2+YC3/T2* @ (TWO*C/GES**2+(B+TWO*C/GES)*(ONE-C/GES**2)/T2) @ + ZMG*DD/(GES+DD)**2 DGES = (GEP-GES)/T1 GES = GES+DGES IF( DABS(DGES)/GES .LT. PREC ) GOTO 3 1 CONTINUE GOTO 12 3 CONTINUE C C ELECTRON PRESSURE C PE = RT*GES/V C C TOTLN = 1/MU = X/C2+Y/C4+Z/C3+GES C TOTLN = PARGRM+GES XX = D/(GES+D) T2 = GES+B+C/GES YY = B/T2 ZZ = C/(GES*T2) WW = DD/(GES+DD) C C DERIVATIVES OF THE SAHA EQUATION FOR THE C PRESSURE AND INTERNAL ENERGY TEMPERATURE AND C DENSITY DERIVATIVES. C T1 = YC3*(B+TWO*C/GES) QC0 = ONE+XC2*XX/(GES+D)+ZMG*WW/(GES+DD)+YC3/T2* @ (TWO*C/GES**2+(B+TWO*C/GES)*(ONE-C/GES**2)/T2) QC1 = XC2*(ONE-XX)/(GES+D) QC4 = ZMG*(ONE-WW)/(GES+DD) QC2 = (YC3-T1/T2)/T2 QC3 = (YC3*TWO-T1/T2)/(GES*T2) QGV = (QC1*D+QC2*B+QC3*TWO*C+QC4*DD)/(QC0*V) QP1 = D*(ONHLF+XH*TK)/T QP2 = B*(ONHLF+XHE*TK)/T QP3 = C*(THRE+(XHE+XHE2)*TK)/T QP4 = DD*(ONHLF+XM*TK)/T QGT = (QC1*QP1+QC2*QP2+QC3*QP3+QC4* QP4)/QC0 C C ELECTRON PRESSURE DERIVATIVES. C PET = ONE + QGT/GES PEV =-ONE + QGV/GES C C PRESSURE DUE TO THE IDEAL GAS C P = RT*TOTLN/V PT = P/T+RT*QGT/V PV = RT*QGV/V-P/V C C BP IS R/MU C BP = R*TOTLN BPV = R*QGV BPT = R*QGT C C ADD THE RADIATION PRESSURE C P = P+AD3*TT4 PT = PT+FOR*AD3*TT4/T C C IONIZATION ENERGY C EI = (R/BK)*(XH*XX*XC2+YC3*(XHE*YY+(XHE+XHE2)*ZZ) @ +ZMG*XM*WW + ZNA*5.524D0 ) C TOTAL INTERNAL ENERGY E = ONHLF*RT*TOTLN + A*V*TT4 + EI EV = T*PT-P QXT = ((ONE-XX)*QP1-XX*QGT)/(GES+D) DT2 = QGT*(ONE-C/GES**2)+QP2+QP3/GES QYT = (QP2-B*DT2/T2)/T2 QZT = (QP3-C*QGT/GES-C*DT2/T2)/(T2*GES) QWT = ((ONE-WW)*QP4-WW*QGT)/(GES+DD) EIT = (R/BK)*(XH*QXT*XC2+YC3*(XHE*QYT+(XHE+XHE2)*QZT)+ @ ZMG*XM*QWT) ET = ONHLF*R*(TOTLN+T*QGT)+FOR*A*V*TT4/T+EIT RETURN 10 WRITE(6,100) T,V 100 FORMAT(1H0,27HNEGATIVE TEMP IN EOS T=,1PE10.3,2X, * 2HV=,E10.3) STOP 11 WRITE(6,101) T,V 101 FORMAT(1H0,27HNEGATIVE DENSITY IN EOS T=,1PE10.3,2X, * 2HV=,E10.3) STOP 12 WRITE(6,102) T,V 102 FORMAT(1H0,27HNO CONVERGENCE IN EOS T=,1PE10.3,2X, * 2HV=,E10.3) STOP END