subroutine EOS(P,t,xin,yin,zin,tcut,q,cp,dela, 1 d,xd,xt,u,cv,eta,fmu) implicit double precision (a-h,o-z) c c C EQUATION OF STATE FOR A MIXTURE OF X (HYDROGEN) Y (HELIUM) C and metals. Pressure and Temperature c are independent variables. Saha equation solver assumes c non-degenerate material, Eggleton Faulkner and Flannery c technique for Fermi-Dirac statistics... c c INPUT: c P = pressure (cgs) c T = temperature (K) c XIN = hydrogen abundance by mass c YIN = helium abundance by mass c Zin = "metal" abundance by mass c Tcut = assumed fully ionized above tcut+0.2; between Tcut c and Tcut+0.2, interpolates between full ionization c and Saha c c OUTPUT: c q = - dln(rho)/dlnT at constant pressure c Cp = dU/dT at constant pressure c dela = adiabatic temperature gradient (=dlnT/dlnP at constant S) c Ro = density (g/cm3) c Xd = dlnP/dlnrho at constant T c Xt = dlnP/dlnT at constant rho c U = internal energy per gram c Cv = dU/dT at constant volume c eta = degeneracy parameter c fmu = mean atomic weight per particle C DOUBLE PRECISION NE,dne1 fln10=2.30258509299404d0 C x=xin y=yin c For starters, set metals as the rest - i.e. ignore input z. z=1.d0-xin-yin c ptl=dlog(p)/fln10 tl=dlog(t)/fln10 iful=0 ne=0.d0 c c First compute Saha-style E.O.S. if TL tcut+0.2, or if Saha says so c or if pressure greater than 1e20 c if (tl.gt.tcut+0.2d0.or.iful.ge.1.or.ptl.gt.20.d0) then q=qi cp=cpi dudtp=dudtpi dela=delai d=di dl=dli xd=xdi xt=xti u=ui ul=uli cv=cvi eta=etai fmu=fmui c c If log(T) between Tcut and Tcut+0.2d0, interpolate between full c ionization and partial ionization for the quantities. c elseif (Tl.gt.tcut.and.tl.lt.Tcut+0.2d0) then fac=(tl-tcut)/0.2d0 q=q+fac*(qi-q) Cp=Cp+fac*(Cpi-Cp) Dela=Dela+fac*(Delai-Dela) dl=dl+fac*(dlog10(di)-dl) c d=10.d0**dl d=dexp(2.302585093d0*dl) Xd=Xd+fac*(Xdi-Xd) Xt=Xt+fac*(Xti-Xt) Ul=Ul+fac*(dlog10(Ui)-Ul) c u=10.d0**Ul u=dexp(2.302585093d0*Ul) Cv=Cv+fac*(Cvi-Cv) eta=eta+fac*(etai-eta) fmu=fmu+fac*(fmui-fmu) endif c c Trouble? Report! c if (xt.lt.0.d0.or.DELA.lt.0.d0.or.cv.lt.0.d0.or.xd.lt.0.d0) then print 131, ptl,tl,x,y,c,o,dlog10(d),xd,q,xt,DELA,cv,cp 1 ,d1,d2,d2-d1,u1,u2,u2-u1,eta1,eta write (31,131) ptl,tl,x,y,c,o,dlog10(d),xd,q,xt,DELA,cv,cp 1 ,d1,d2,d2-d1,u1,u2,u2-u1,eta1,eta 131 format(' Trouble in EOS land!'/' log P=',f9.5,'log T=',f9.5, 1 ' X,Y,C,0 :',4f7.4,/,' log ro =',f9.5/ 1 ' Xd=',f9.5,' Q=',f9.5,' Xt=',f9.5/ 2 ' DelAd=',f9.5,' Cv=',1pe10.3,' Cp=',1pe10.3/ 3 ' d1 =',e10.3,' d2=',e10.3,' d2-d1=',e10.3/ 4 ' U1 =',e10.3,' U2=',e10.3,' U2-U1=',e10.3/ 5 ' ETA1=',e10.3,' ETA=',e10.3) endif 100 format(3f8.4,1p3e12.4,0p2f9.5) c c Okay, all done! c RETURN END C SUBROUTINE STATE(PTL,TL,AX1,AX2,AX3,IFUL,D,UT,ETA,MU,NE,ITNE) IMPLICIT DOUBLE PRECISION (A-H,O-Z) c c Calculates equation of state using Saha equation a la Mihalas for c partial ionization. c C Does not include pressure ionization (except via temperature cutoffs). c C INPUTS: PTL = LOG TOTAL PRESSURE C TL = LOG TEMPERATURE C AX1 = MASS FRACTION OF HYDROGEN C AX2 = MASS FRACTION OF HELIUM C AX3 = MASS FRACTION OF 'METAL' c NE = FIRST GUESS AT ELECRON DENSITY c C OUTPUTS: C D = DENSITY C UT = INTERNAL ENERGY PER GRAM C ETA = DEGENERACY PARAMETER c Mu = MEAN ATOMIC WEIGHT PER PARTICLE C NE = NUMBER DENSITY OF ELECTRONS C ITNE = NUMBER OF ITERATIONS ON NE C IFUL = SET TO 1 IF FINDS COMPLETE IONIZATION C DOUBLE PRECISION Y(0:10,5),NE,MUI,MUE,X(3),MU DOUBLE PRECISION Netot COMMON/ATDAT/G(0:2,3),XI(0:2,3),A(3),IZ(3) c Common block CONST contains fundamental physical constants COMMON/CONST/CMS,CRS,CLS,CG,CSB,CPI,CP4,CCL,CME,CH,CNA,CEV,CK,CSYR DATA Y/55*0.D0/ c CSAHA=(2.D0*CPI*CK*CME)**(3.D0/2.D0) c CSAHA=CSAHA/CH/CH/CH csaha=2.414703187d15 C NSPEC=3 X(1)=AX1 X(2)=AX2 X(3)=AX3 c c T=10.D0**(TL) T=dexp(2.302585093d0*TL) c Ptot=10.D0**(PTL) Ptot=dexp(2.302585093d0*PTL) Prad=2.5219D-15*(T**4) Pgas=PTOT-PRAD if (pgas.le.0.d0) then print *,' Total pressure less than gas pressure?!' print '(A,f10.5,a,f10.5)','log P=',ptl,'log T=',tl print *,' Setting it equal to 0.01*P, but BEWARE' pgas=0.01d0*ptot endif PgasL=dlog10(Pgas) THRES=1.d-5 c If X+Y+Z doesn't add up to one, fudge with a heavy, atomic weight c of 20 and charge of 1. sumx=x(1)+x(2)+x(3) if (sumx.lt.1.d0) then X(3)=x(3)+1.d0-sumx endif C C COMPUTE . . . C AVZOA, the mass-weighted c MUI, the mean atomic weight per nucleus C MUE, the mean atomic weight per electron, free or bound OOMUI=0.D0 AVZOA=0.D0 DO 5 I=1,NSPEC OOMUI=OOMUI+X(I)/A(I) AVZOA=AVZOA+X(I)*IZ(I)/A(I) 5 CONTINUE c MUI=1.D0/OOMUI MUE=1.d0/avzoa C C COMPUTE NE C Use the Saha equation to find the new NE without pressure ionization c CALL SAHA(PGASL,DL,TL,NSPEC,X,NE,Y,ITNE) D=MUI*(Pgas/cna/ck/t-Ne/cna) dl=dlog10(d) c Now is a good time to compute Mu, the mean atomic weight per free particle fmuefr=D*CNA/NE mu=1.d0/(1.d0/mui+1.d0/fmuefr) c c Check for "complete" ionization c NETOT=D/MUe*CNA c print '(10hNe/Netot ,1pe12.5)',ne/netot IF (NE/NETOT.gt.0.999D0) IFUL=1 C C CACULATE INTERNAL ENERGY PER GRAM C C Calculate properties of e- gas. PPGE=NE*CK*T UPGE=1.5d0*ppge/d C RADIATION URAD=3.D0*PRAD/D C TRANSLATIONAL PPGI=D*CNA*CK*T/MUI UGT=UPGE+1.5D0*PPGI/D C IONIZATION ENERGY FOR EACH SPECIES UGI=0.D0 DO 50 K=1,NSPEC IF (X(K).GE.1.D-10) THEN UGII=0.D0 DO 60 IJ=1,IZ(K) DO 70 JJ=0,IJ-1 UGII=UGII+X(K)*Y(IJ,K)*XI(JJ,K)*CEV 70 CONTINUE 60 CONTINUE UGII=UGII*CNA/A(K) UGI=UGI+UGII ENDIF 50 CONTINUE c c Add in dissosiation energy of H2, assume always dissociated XIH2=4.477d0 UH2=0.5d0*X(1)*CEV*XIH2*CNA UGI=UGI+UH2 C TOTAL INTERNAL ENERGY PER GRAM c print '(1pe12.5/e12.5/e12.5/e12.5)',upge,1.5d0*ppgi/d,ugt,ugi UT=URAD+UGT+UGI c print '(1p5e14.6)',upge,1.5d0*ppgi/d,ugt,ugi,urad c Degeneracy parameter in nondegenerate limit ETA=-36.113d0+dlog(Ne)-3.4539d0*TL c RETURN END C subroutine SAHA(PGL,DL,TL,NSPEC,X,NE,Y,ITNE) implicit double precision (a-h,o-z) c c Subroutine to compute Ne for an ensemble of NSPEC atomic species c NSTAGE in the parameter statement is the maximum number of ionization c stages expected. Can be called with either (Pgas,T) or (Rho,T) c as the independent variables. c c If the sum of the abundances is less than 1, it adds a final species c that is a one-electron metal c c Follows algorithm in Mihalas (around p. 118) but with (Pg,T) as the c known quantities. Thus the linearization is slightly different, c but all notation is the same as Mihalas. c c Input: PGL = log gas pressure c note: if PGL is input as greater than 90 c then this routine c uses DL as the first independent variable c DL = log density c TL = log temperature c NSPEC = number of species to consider c X = array of mass fractions c c Output: Ne = number density of free electrons c Y(Nstage,NSPEC)= fraction of given stage that is ionized c ITNE = number of iterations performed c parameter(NSTAGE=10) double precision Ne,Ne0,Ntot,Netot,NNuc dimension Y(0:10,5),X(5) dimension phi(0:10,5),prophi(5) dimension P(0:nstage,5),S(5),alfa(5) COMMON/ATDAT/G(0:2,3),XI(0:2,3),A(3),IZ(3) c Common block CONST contains fundamental physical constants COMMON/CONST/CMS,CRS,CLS,CG,CSB,CPI,CP4,CCL,CME,CH,CNA,CEV,CK,CSYR c c first index in, i.e. XI(J,K) is ionization stage, second is species c c read in atomic and other parameters if this is the first call c CSAHA=CH**2/2.d0/CPI/CME/CK CSAHA=0.5d0*dexp(1.5d0*dlog(CSAHA)) CKEV=CK/CEV THRESH=1.d-5 C c now to work! c c T=10.d0**TL T=dexp(2.302585093d0*TL) FKT=CKEV*T c TM32=10.d0**(-3.d0*TL/2.d0) TM32=dexp(-3.d0*(TL*2.302585093d0)/2.d0) CSTM32=CSAHA*TM32 c compute total number fraction of each species UAm1=0.d0 nK=nspec do 5 k=1,nspec UAm1=UAm1+x(k)/A(k) 5 continue do 7 k=1,nK alfa(k)=X(k)/A(k)/UAm1 7 continue c print '(1p5e12.5)',(alfa(k),k=1,nK) c c if Pgas is an input variable, compute Ntot, the total c number density of all paritcles, from the Pressure c if (PGL.lt.90.d0) then c Pgas=10.d0**PGL Pgas=dexp(2.302585093d0*PGL) Ntot=Pgas/CK/T if (ntot.le.0.d0) ntot=1.d5 c and use this to make a first guess at Ne0 Ne0=Ntot/100.0 else c c otherwise, use the input density to determine c the number density of nuclei c rho=10.d0**DL rho=dexp(2.302585093d0*DL) Nnuc=rho*uam1*cna c and the total number of electrons, free or bound fnetot=0.d0 do 9 k=1,nk fnetot=fnetot+alfa(k)*iz(k) c print '(i4,2f10.4,1pe12.4,0pi4,1pe12.4)', c 1 k,pgl,tl,alfa(k),iz(k),fnetot 9 continue netot=nnuc*fnetot Ne0=Netot/100.d0 endif c c More preliminaries: phi's and stuff do 115 k=1,nk if (x(k).le.1.d-10) go to 115 prophi(k)=1.d0 do 125 j=iz(k)-1,0,-1 flpjk=dlog(g(j,k)/g(j+1,k)*csaha*tm32)+(xi(j,k)/fkt) c don't let the exponential factor get too big... causes trouble in unionized c material... if (flpjk.gt.28.d0) flpjk=28 phi(j,k)=dexp(flpjk) prophi(k)=prophi(k)*phi(j,k) 125 continue 115 continue C C ITERATE ON NE: c compute the P's and S's. Note Mihalas's convention that the c product terms involving the last ionization stage ie. P(iz(k),k) c are unity. c itne=0 verg=1.d0 do while (itne.le.40.and.verg.ge.thresh) itne=itne+1 sumbar=0.d0 c dsbdNe is the derivative of sumbar w.r.t. Ne0 dsbdNe=0.d0 do 15 k=1,nK S(k)=1.d0 if (x(k).le.1.d-10) go to 15 P(iz(k),k)=1.d0 sumjP=iz(k)*1.d0 sumjdP=0.d0 dSkdNe=0.d0 do 25 j=iz(k)-1,0,-1 P(j,k)=P(j+1,k)*ne0*phi(j,k) sumjP=sumjP+j*P(j,k) c dPjkdNe is the derivative of P(J,K) w.r.t. Ne0 dPjkdNe=(iz(k)-j)*ne0**(iz(k)-j-1)*prophi(k) sumjdP=sumjdP+j*dPjkdNe S(k)=S(k)+P(j,k) dSkdNe=dSkdNe+dPjkdNe 25 continue sumbar=sumbar+alfa(k)*sumjP/S(k) dsbdNe=dsbdNe+alfa(k)*(sumjdP/S(k)-sumjP/S(k)*dSkdNe/S(k)) c 15 continue c c the new computed value of Ne c if (PGL.lt.90.d0) then ne=(Ntot-ne0)*sumbar else ne=nnuc*sumbar endif c c print *,iter,' HII', P(1,1)/S(1),' HI',P(0,1)/S(1) c print *,iter,' HeII', P(1,2)/S(2),' HeI',P(0,2)/S(2) c print *,iter,' HeIII', P(2,2)/S(2) c print *,iter,' XMII', P(1,3)/S(3),'XMI',P(0,3)/S(3) c print * c c Compute a new (hopefully better) value for Ne: c Linear correction to ne0 c if (PGL.lt.90.d0) then dne=(ne-ne0)/(1.d0+sumbar-(Ntot-ne0)*dsbdNe) else dne=(ne-ne0)/(1.d0-nnuc*dsbdNe) endif c print '(I3,1p3e13.5,3x,2e12.5)',iter,ne0,ne,dne/ne0,netot c c CONVERGENCE CRITERIA: c difference between computed and guessed value of ne verg=dabs((ne-ne0)/ne0) c electron density/total particle density if (ntot.le.0.d0) then frace=0.5d0 else frace=dabs(ne0/ntot) endif verg=dmin1(frace*100.d0,verg) c Add the correction. c if the correction is too large, reduce it a tad. This is essential c when ne/netot is tiny and the correction tries to subtract the c value of ne0 from itself (to within machine accuracy) fac=dabs(dne/ne0) if (fac.ge.0.1d0) dne=(1.d0-1.d-6)*dne c and make the correction ne0=ne0+dne if (ne0.le.0.d0) ne0=1.d0 c If iteration counter is getting too large, then reduce convergence c threshold. if (itne.eq.20) then write (31,131) PGL,TL,X(1),X(2),X(3) write (31,'(I3,1p3e13.5,3x,e12.5)') iter,ne0,ne,dne/ne0 131 format(' Saha routine in trouble : PGL=',f9.4,' TL=',f9.4, 1 ' X=',f7.4,' Y=',f7.4,' C=',f7.4/ 'Reducing thresh') thresh=thresh*10.d0 endif if (itne.eq.30) then write (31,132) PGL,TL write (31,'(I3,1p3e13.5,3x,e12.5)') iter,ne0,ne,dne/ne0 132 format(' Saha routine in trouble again : PGL=',f9.4, 1 ' TL=',f9.4) thresh=thresh*10.d0 endif c ******** End of iteration loop ***************** enddo c c If electron density small, return with a guess at it based on slight c ionization if (frace.le.thresh) ne0=dsqrt(Pgas/ck/t) c store the new value of ne in the correct variable ne=ne0 c c compute fractional ionization do 28 k=1,nspec do 27 j=0,iz(k) y(j,k)=p(j,k)/s(k) 27 continue 28 continue c return end C block data c SUBROUTINE RATDAT C C ATOMIC DATA FOR THE CHOSEN SPECIES C IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/ATDAT/G(0:2,3),XI(0:2,3),A(3),IZ(3) DATA A/1.00782505d0,4.00260308d0,14.d0/ DATA IZ/1,2,1/ DATA G/2.d0,1.d0,0.d0, 1 1.d0,2.d0,1.d0, 2 2.d0,1.d0,0.d0/ DATA XI/13.598d0,0.d0,0.d0, 1 24.587d0,54.416d0,0.d0, 2 5.000d0,0.0d0,0.0d0/ c a stub c temp=1234.56d0 c return c END c subroutine IONIZED(PL,TL,X,Y,z,d,U,eta,fmu,Pe,Pi,ue, 1 Xd,q,dudtP) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C FULLY IONIZED EQUATION OF STATE: ITERATE ON DENSITY TO GIVE PROPER c TOTAL PRESSURE GIVEN THAT ALL ELECTRONS ARE FREE c C INPUT: PL = log(total pressure) c TL = log(temperature) c X = hydrogen mass fraction c Y = helium mass fraction c C = carbon mass fraction c c OUTPUT: d = density c U = total energy density c eta = degeneracy parameter c fmu = mean atomic weight per free particle c Pe = electron pressure c Pi = ion pressure c Ue = energy density in electron gas c Xd = dlnP/dlnrho at constant T c q = -dlnrho/dlnT at constant P c dudtP = dU/dT at constant P C c Common block CONST contains fundamental physical constants COMMON/CONST/CMS,CRS,CLS,CG,CSB,CPI,CP4,CCL,CME,CH,CNA,CEV,CK,CSYR DOUBLE PRECISION mue,mui c c Convergence threshold for fully ionized iterations c THRESH=1.d-5 c c Some preliminaries: subtract off radiation pressure to obtain c gas pressure alone c Ptot=10.d0**pl Ptot=dexp(2.302585093d0*pl) c T=10.d0**tl T=dexp(2.302585093d0*tl) Prad=2.5219d-15*t**4 Pgas=Ptot-Prad c Assume that anything left after x,y,c is o Z=1.d0-X-Y c Mean atomic weight per ion and per electron (free or bound) oomui=X/1.007825d0 + Y/4.002603d0 + Z/14.0d0 oomue=X/1.007825d0 + 2.d0*Y/4.002603d0 + 0.5d0*Z mue=1.d0/oomue mui=1.d0/oomui fmu=1.d0/(oomui+oomue) c Sum of ionization potentials in eV, for energy computation sumxix=13.598d0 sumxiy=79.003d0 sumxiz=1500.d0 C First guess at D0 and NE: assume nondegenerate gas with no Coulomb term D0=Pgas/CNA/CK/T/(1.D0/mui+1.D0/mue) c Don't allow the first guess of D0 to exceed 1e6, otherwise relativistic c corrections go haywire if (d0.gt.1.d6) d0=1.d6 c ************************************************************************ c C Iterations on D to get Pe+Pi=Pgas; continue iterating until density c changes by less than THRESH and derived Pgas matches true value c drho=d0 itfi=0 DO while ((dabs(drho/d0).gt.THRESH 1 .OR.dabs(dpgas/pgas).gt.THRESH) 1 .AND.itfi.lt.30) itfi=itfi+1 c Pressure of ions assumed to be a perfect gas Pi=d0/mui*CNA*CK*T C Compute pressure of electron gas given guess at domue via EFF DOmue0=D0/mue call EFF(domue0,T,eta,Pe,Ued,dpedom,dpedt,duedtr) c Compute new value of total gas pressure with this guess c at the density. Pgas1=Pe+Pi dPgas=Pgas-Pgas1 c Compute the corresponding linearized correction to the density c dPe/drho: dpedro=dpedom/mue c dPi/drho: dpidro=CNA*CK*T/mui c Use these derivatives to compute the linearized density correction fac=dpidro+dpedro drho=dpgas/fac c reduce the density correction if it is very large if (dabs(drho/d0).ge.1.d0) then c write (31,*)' Reducing drho by 50%, T=',t,' Rho=',d0,' Iter=' c 1 ,itfi drho=dsign(d0/2.d0,drho) elseif (dabs(drho/d0).gt.0.1d0) then c write (31,*)' Reducing drho by 10%, T=',t,' Rho=',d0,' Iter=' c 1 ,itfi drho=drho*0.9d0 endif c Apply the density correction d1=d0+drho c print *,itfi,d1,drho/d1 d0=d1 c enddo c C Done iterating on D for full ionization c ******************************************************************** c All done; set up for return d=d0 IF (ITFI.GT.29) then PRINT *, ' Fully Ionized Iterations failed to converge' print *, ' T=',dlog10(t),' P=',dlog10(PL) endif c c Internal energy c c Radiation pressure energy density Urad=3.d0*Prad/d c Energy DENSITY of electron gas Ue=Ued/d c Translational energy density of ions Ui=1.5d0*Pi/d c Excitation energy density Uex=cna*(X*sumxix+Y*sumxiy/4.d0+z*sumxiz/14.d0) Uex=Uex*cev c Total energy U=Urad+Ue+Ui+Uex c c Compute derivatives analytically c pressure derivatives dpidt=d*cna*ck/mui dpidro=cna*ck*t/mui dpedro=dpedom/mue dprdro=0.d0 dprdt=4.d0*prad/t c dpdro=dpidro+dpedro+dprdro dpdt=dpidt+dpedt+dprdt drodt=-dpdt/dpdro c c Energy derivative w.r.t. T at constant rho duidt=1.5d0*cna*ck/mui c c print *,"Gam=",gam c duedt=duedtr/d durdt=4.d0*3.d0*prad/t/d c dudtr=duidt+duedt++durdt dudtp=dudtr-t/d**2*dpdt*drodt+ptot/d**2*drodt c ********************************************************************** c Derivatives for output xd=dpdro/ptot*d xt=dpdt/ptot*T q=xt/xd return end c subroutine EFF(domue,T,psi,Pe,Ue,dPeDoM,dPedT,dUedTr) implicit double precision (a-h,o-z) c dimension rhat(4,4),phat(4,4),uhat(4,4) data rhat/2.315472d0,7.128660d0,7.504998d0,2.665350d0, 1 7.837752d0,23.507934d0,23.311317d0,7.987465d0, 2 9.215560d0,26.834068d0,25.082745d0,8.020509d0, 3 3.693280d0,10.333176d0,9.168960d0,2.668248d0/ data phat/2.315472d0,6.748104d0,6.564912d0,2.132280d0, 1 7.837752d0,21.439740d0,19.080088d0,5.478100d0, 2 9.215560d0,23.551504d0,19.015888d0,4.679944d0, 3 3.693280d0,8.859868d0,6.500712d0,1.334124d0/ data uhat/3.473208d0,10.122156d0,9.847368d0,3.198420d0, 1 16.121172d0,43.477194d0,37.852852d0,10.496830d0, 2 23.971040d0,60.392810d0,47.782844d0,11.361074d0, 3 11.079840d0,26.579604d0,19.502136d0,4.002372d0/ M=4 N=4 c c given Domue and T, compute rho/mue, Pe, Ue using the thermodynamically c consistent interpolation formulae of Eggleton Faulkner and Flannery (1973) c also returns analytic values of d(Pe)/d(rho/mue) at constant T, c d(Pe)/dT at constant rho, and d(Ue)/dT at constant rho c c note that Ue is internal energy per unit VOLUME c flcomp is the compton wavelength of the electron c flcomp=2.42580d-10 flcomp=2.42631d-10 fac=flcomp**3/8.d0/3.1415927d0/9.1094d-28/2.998d10**2 c beta=t/9.1094d-28/2.998d10**2*1.3807d-16 rhoh0=domue*flcomp**3/8.d0/3.141592654d0/1.66054d-24 c c Obtain f from domue via Newton's method using EFF eq. 15 c c f0=domue/beta**1.5*1.693d-7 f0=domue/dexp(1.5d0*dlog(beta))*1.693d-7 do 10 it=1,20 f=f0 g=beta*dsqrt(1.d0+f) dbdfg=-beta/2.d0/(1.d0+f) dbdgf=beta/g dgdfb=beta/2.d0/dsqrt(1.d0+f) c Now evaluate the polynomial for the rhohat corresponding to c the guessed f (and g) sumr=0.d0 sumdrf=0.d0 sumdrg=0.d0 do 20 im=1,M do 30 in=1,N fmgn=f**(im-1)*g**(in-1) sumr=sumr+rhat(im,in)*fmgn sumdrf=sumdrf+rhat(im,in)*fmgn*(im-1.d0)/f sumdrg=sumdrg+rhat(im,in)*fmgn*(in-1.d0)/g 30 continue 20 continue fnorm=(1.d0+f)**(M-1)*(1.d0+g)**(N-1) c g32=dexp(1.5d0*dlog(g)) ggp132=dexp(1.5d0*dlog(g+g*g)) c rhohat=f/(1.d0+f)*g**1.5d0*(1.d0+g)**1.5d0*sumr/fnorm c rhohat=f/(1.d0+f)*g32*gp132*sumr/fnorm rhohat=f/(1.d0+f)*ggp132*sumr/fnorm drdfg=rhohat*(1.d0/f/(1.d0+f)+sumdrf/sumr-(m-1.d0)/(1.d0+f)) drdgf=rhohat* 1 (1.5d0/g+1.5d0/(1.d0+g)+sumdrg/sumr-(n-1.d0)/(1.d0+g)) drdfb=(drdfg*dbdgf-drdgf*dbdfg)/dbdgf c compute the correction to f delta=(rhohat-rhoh0)/drdfb c print *,it,f,delta/f f0=f-delta c if the correction is small enough, then u.r. done if (dabs(delta/f).lt.1.d-7) go to 11 10 continue c 11 if (it.ge.20) then print *, ' E.F.F. iterations failed to converge' write (31,*) ' E.F.F. iterations failed to converge' write (31,*) ' log rho/mue=',dlog10(domue),' log T=',dlog10(t) endif c f=f0 c if (f.lt.0.d0) print *,' EFF: f=',f0 c if (f.lt.0.d0) f=1.d-10 g=beta*dsqrt(1.d0+f) ggp132=dexp(1.5d0*dlog(g+g*g)) dbdfg=-beta/(2.d0*(1.d0+f)) dbdgf=beta/g c c now use the value of f to compute Pe and Ue and Psi c sump=0.d0 sumu=0.d0 sumdpf=0.d0 sumdpg=0.d0 sumduf=0.d0 sumdug=0.d0 do 40 im=1,M do 50 in=1,N fmgn=f**(im-1)*g**(in-1) sump=sump+phat(im,in)*fmgn sumdpf=sumdpf+phat(im,in)*fmgn*(im-1.d0)/f sumdpg=sumdpg+phat(im,in)*fmgn*(in-1.d0)/g sumu=sumu+uhat(im,in)*fmgn sumduf=sumduf+uhat(im,in)*fmgn*(im-1.d0)/f sumdug=sumdug+uhat(im,in)*fmgn*(in-1.d0)/g 50 continue 40 continue fnorm=(1.d0+f)**(M-1)*(1.d0+g)**(N-1) c x=dsqrt(1.d0+f) psi=2.d0*x+dlog((x-1.d0)/(x+1.d0)) c g52=dexp(2.5d0*dlog(g)) c gp132=dexp(1.5d0*dlog(g+1.d0)) g52p132=g*ggp132 c PeHat=f/(1.d0+f)*g**2.5d0*(1.d0+g)**1.5d0*sump/fnorm c UeHat=f/(1.d0+f)*g**2.5d0*(1.d0+g)**1.5d0*sumu/fnorm c PeHat=f/(1.d0+f)*g52*gp132*sump/fnorm PeHat=f/(1.d0+f)*g52p132*sump/fnorm c UeHat=f/(1.d0+f)*g52*gp132*sumu/fnorm UeHat=f/(1.d0+f)*g52p132*sumu/fnorm c c compute derivatives w.r.t. f at constant g dpdfg=PeHat*(1.d0/f/(1.d0+f)+sumdpf/sump-(m-1)/(1.d0+f)) dudfg=UeHat*(1.d0/f/(1.d0+f)+sumduf/sumu-(m-1)/(1.d0+f)) c and derivatives w.r.t. g at constant f dpdgf=PeHat*(2.5d0/g+1.5d0/(1.d0+g)+sumdpg/sump-(n-1)/(1.d0+g)) dudgf=UeHat*(2.5d0/g+1.5d0/(1.d0+g)+sumdug/sumu-(n-1)/(1.d0+g)) c compute due/dT at constant rho = duhhat/dbeta*8*pi*k/flam**3 dudbr=(dudfg*drdgf-dudgf*drdfg)/(dbdfg*drdgf-dbdgf*drdfg) dUedTr=dudbr*8.d0*3.14159254d0*1.38d-16/flcomp**3 c now dPe/d(rho/me) at constant beta dPhDrh=(dpdgf*dbdfg-dpdfg*dbdgf)/(dbdfg*drdgf-drdfg*dbdgf) dPedOm=dphdrh*9.109d-28*2.998d10**2/1.6724d-24 c and dPe/dT:rho=(dPh/dbeta)*8*pi*k/flam**3 dphdb=(dpdfg*drdgf-dpdgf*drdfg)/(dbdfg*drdgf-dbdgf*drdfg) dPedT=dPhdb*8.d0*3.141592654d0*1.38d-16/flcomp**3 c Pe=PeHat/fac Ue=UeHat/fac c return end