Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2022 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module pc_eos
21 : use utils_lib
22 : use math_lib
23 : use auto_diff
24 : use const_def, only: dp, PI
25 :
26 : implicit none
27 :
28 : contains
29 :
30 : !* Equation of state for fully ionized electron-ion plasmas (EOS EIP)
31 : ! A.Y.Potekhin & G.Chabrier, Contrib. Plasma Phys., 50 (2010) 82,
32 : ! and references therein
33 : ! Please communicate comments/suggestions to Alexander Potekhin:
34 : ! palex@astro.ioffe.ru
35 : ! Previously distributed versions (obsolete):
36 : ! eos2000, eos2002, eos2004, eos2006, eos2007, eos2009, eos10, eos11,
37 : ! and eos13.
38 : ! Last update: 28.05.15. All updates since 2008 are listed below.
39 : !* L I S T O F S U B R O U T I N E S :
40 : ! MAIN (normally commented-out) - example driving routine.
41 : ! MELANGE9 - for arbitrary ionic mixture, renders total (ion+electron)
42 : ! pressure, internal energy, entropy, heat capacity (all
43 : ! normalized to the ionic ideal-gas values), logarithmic
44 : ! derivatives of pressure over temperature and density.
45 : ! EOSFI8 - nonideal (ion-ion + ion-electron + electron-electron)
46 : ! contributions to the free and internal energies, pressure,
47 : ! entropy, heat capacity, derivatives of pressure over
48 : ! logarithm of temperature and over logarithm of density (all
49 : ! normalized to the ionic ideal-gas values) for one ionic
50 : ! component in a mixture.
51 : ! FITION9 - ion-ion interaction contributions to the free and internal
52 : ! energies, pressure, entropy, heat capacity, derivatives of
53 : ! pressure over logarithms of temperature and density.
54 : ! FSCRliq8 - ion-electron (screening) contributions to the free and
55 : ! internal energies, pressure, entropy, heat capacity,
56 : ! derivatives of pressure over logarithms of temperature and
57 : ! density in the liquid phase for one ionic component in a
58 : ! mixture.
59 : ! FSCRsol8 - ion-electron (screening) contributions to the free and
60 : ! internal energies, pressure, entropy, heat capacity,
61 : ! derivatives of pressure over logarithms of temperature and
62 : ! density for monoionic solid.
63 : ! FHARM12 - harmonic (including static-lattice and zero-point)
64 : ! contributions to the free and internal energies, pressure,
65 : ! entropy, heat capacity, derivatives of pressure over
66 : ! logarithms of temperature and density for solid OCP.
67 : ! HLfit12 - the same as FHARM12, but only for thermal contributions
68 : ! ANHARM8 - anharmonic contributions to the free and internal energies,
69 : ! pressure, entropy, heat capacity, derivatives of pressure
70 : ! over logarithms of temperature and density for solid OCP.
71 : ! CORMIX - correction to the linear mixing rule for the Coulomb
72 : ! contributions to the thermodynamic functions in the liquid.
73 : ! ELECT11 - for an ideal electron gas of arbitrary degeneracy and
74 : ! relativity at given temperature and electron chemical
75 : ! potential, renders number density (in atomic units), free
76 : ! energy, pressure, internal energy, entropy, heat capacity
77 : ! (normalized to the electron ideal-gas values), logarithmic
78 : ! derivatives of pressure over temperature and density.
79 : ! EXCOR7 - electron-electron (exchange-correlation) contributions to
80 : ! the free and internal energies, pressure, entropy, heat
81 : ! capacity, derivatives of pressure over logarithm of
82 : ! temperature and over logarithm of density (all normalized
83 : ! to the classical electron ideal-gas values).
84 : ! FERINV7 - inverse non-relativistic Fermi integrals of orders -1/2,
85 : ! 1/2, 3/2, 5/2, and their first and second derivatives.
86 : ! BLIN9 - relativistic Fermi-Dirac integrals of orders 1/2, 3/2, 5/2,
87 : ! and their first, second, and some third derivatives.
88 : ! CHEMFIT7 - electron chemical potential at given density and
89 : ! temperature, and its first derivatives over density and
90 : ! temperature and the second derivative over temperature.
91 : !* I M P R O V E M E N T S O F 2 0 0 8 - 2 0 0 9 :
92 : ! FHARM8 uses a fit HLfit8 to the thermal free energy of the harmonic
93 : ! Coulomb lattice, which is more accurate than its predecessor FHARM7.
94 : ! Resulting corrections amount up to 20% for the ion heat capacity.
95 : ! Accordingly, S/R D3fit and FthCHA7 deleted (not used anymore).
96 : ! BLIN7 upgraded to BLIN8:
97 : ! - cleaned (a never-reached if-else branch deleted);
98 : ! - Sommerfeld (high-\chi) expansion improved;
99 : ! - some third derivatives added.
100 : ! CORMIX added (and MELANGE7 upgraded to MELANGE8 accordingly).
101 : ! ANHARM7 upgraded to ANHARM8, more consistent with Wigner-Kirkwood.
102 : ! Since the T- and rho-dependences of individual Z values in a mixture
103 : ! are not considered, the corresponding inputs (AYLR, AYLT) are
104 : ! excluded from MELANGE8 (and EOSFI7 changed to EOSFI8 accordingly).
105 : ! ELECT7 upgraded to ELECT9 (high-degeneracy behaviour is improved)
106 : !* P O S T - P U B L I C A T I O N (2 0 1 0 +) IMPROVEMENTS :
107 : ! ELECT9 upgraded (smooth match of two fits at chi >> 1)
108 : ! BLIN8 replaced by BLIN9 - smooth fit interfaces at chi=0.6 and 14.
109 : ! MELANGE8 replaced by MELANGE9 - slightly modified input/output
110 : ! 08.08.11 - corrected mistake (unlikely to have an effect) in CHEMFIT7
111 : ! 16.11.11 - ELECT9 upgraded to ELECT11 (additional output)
112 : ! 20.04.12 - FHARM8 and HLfit8 upgraded to FHARM12 and HLfit12:
113 : ! output of HLfit12 does not include zero-point vibr., but provides U1
114 : ! 22.12.12 - MELANGE9 now includes a correction to the linear mixing
115 : ! rule (LMR) for the Madelung energy in the random bcc multi-ion
116 : ! lattice.
117 : ! 14.05.13 - an accidental error in programming the newly introduced
118 : ! correction to the LMR is fixed.
119 : ! 20.05.13 - calculation of the Wigner-Kirkwood quantum diffraction term
120 : ! for the liquid plasma is moved from EOSFI8 into MELANGE9.
121 : ! 10.12.14 - slight cleaning of the text (no effect on the results)
122 : ! 28.05.15 - an accidental error in Wigner-Kirkwood entropy correction
123 : ! is fixed (it was in the line "Stot=Stot+FWK*DENSI" since 20.05.13).
124 : ! ***********************************************************************
125 : ! MAIN program: Version 02.06.09
126 : ! This driving routine allows one to compile and run this code "as is".
127 : ! In practice, however, one usually needs to link subroutines from this
128 : ! file to another (external) code, therefore the MAIN program is
129 : ! normally commented-out.
130 : ! C%C implicit double precision (A-H), double precision (O-Z)
131 : ! C%C parameter(MAXY=10,UN_T6=.3157746,EPS=1.d-7)
132 : ! C%C dimension AY(MAXY),AZion(MAXY),ACMI(MAXY)
133 : ! C%C write(*,'('' Introduce the chemical composition (up to'',I3,
134 : ! C%C * '' ion species):''/
135 : ! C%C * '' charge number Z_i, atomic weight A_i,'',
136 : ! C%C * '' partial number density x_i, derivatives d x_i / d ln T'',
137 : ! C%C * '' and d x_i / d ln rho''/
138 : ! C%C / '' (non-positive Z, A, or x=1 terminates the input)'')') MAXY
139 : ! C%C NMIX=0
140 : ! C%C 3 continue
141 : ! C%C XSUM=0.
142 : ! C%C do IX=1,MAXY
143 : ! C%C write(*,'(''Z, A ('',I2,''): ''$)') IX
144 : ! C%C read*,AZion(IX),ACMI(IX)
145 : ! C%C if (AZion(IX) <= 0..or.ACMI(IX) <= 0.) GOTO 2
146 : ! C%C write(*,'(''x ('',I2,''): ''$)') IX
147 : ! C%C read*,AY(IX)
148 : ! C%C XSUM=XSUM+AY(IX)
149 : ! C%C if (AY(IX) <= 0.) GOTO 2
150 : ! C%C NMIX=IX
151 : ! C%C if (dabs(XSUM-1.d0) < EPS) GOTO 2
152 : ! C%C end do
153 : ! C%C 2 continue
154 : ! C%C if (NMIX == 0) then
155 : ! C%C print*,'There must be at least one set of positive (x,Z,A).'
156 : ! C%C GOTO 3
157 : ! C%C end if
158 : ! C%C write(*,114)
159 : ! C%C do IX=1,NMIX
160 : ! C%C write(*,113) IX,AZion(IX),ACMI(IX),AY(IX)
161 : ! C%C end do
162 : ! C%C 9 continue
163 : ! C%C write(*,'('' Input T (K) (<0 to stop): ''$)')
164 : ! C%C read*,T
165 : ! C%C if (T <= 0.) stop
166 : ! C%C 10 continue
167 : ! C%C write(*,'('' Input RHO [g/cc] (<0 to new T): ''$)')
168 : ! C%C read*,RHO
169 : ! C%C if (RHO <= 0.) GOTO 9
170 : ! C%C RHOlg=dlog10(RHO)
171 : ! C%C Tlg=dlog10(T)
172 : ! C%C T6=10.d0**(Tlg-6.d0)
173 : ! C%C RHO=10.d0**RHOlg
174 : ! C%C write(*,112)
175 : ! C%C 1 continue
176 : ! C%C TEMP=T6/UN_T6 ! T [au]
177 : ! C%C call MELANGE9(NMIX,AY,AZion,ACMI,RHO,TEMP,150.0,175.0, ! input
178 : ! C%C * PRADnkT, ! additional output - radiative pressure
179 : ! C%C * DENS,Zmean,CMImean,Z2mean,GAMI,CHI,TPT,LIQSOL, ! output param.
180 : ! C%C * PnkT,UNkT,SNk,CV,CHIR,CHIT) ! output dimensionless TD functions
181 : ! C%C Tnk=8.31447d13/CMImean*RHO*T6 ! n_i kT [erg/cc]
182 : ! C%C P=PnkT*Tnk/1.d12 ! P [Mbar]
183 : ! C%C TEGRAD=CHIT/(CHIT**2+CHIR*CV/PnkT) ! from Maxwell relat.
184 : ! -------------------- OUTPUT -------------------------------- *
185 : ! Here in the output we have:
186 : ! RHO - mass density in g/cc
187 : ! P - total pressure in Mbar (i.e. in 1.e12 dyn/cm^2)
188 : ! PnkT=P/nkT, where n is the number density of ions, T temperature
189 : ! CV - heat capacity at constant volume, divided by number of ions, /k
190 : ! CHIT - logarithmic derivative of pressure \chi_T
191 : ! CHIR - logarithmic derivative of pressure \chi_\rho
192 : ! UNkT - internal energy divided by NkT, N being the number of ions
193 : ! SNk - entropy divided by number of ions, /k
194 : ! GAMI - ionic Coulomb coupling parameter
195 : ! TPT=T_p/T, where T_p is the ion plasma temperature
196 : ! CHI - electron chemical potential, divided by kT
197 : ! LIQSOL = 0 in the liquid state, = 1 in the solid state
198 : ! C%C write(*,111) RHO,T6,P,PnkT,CV,CHIT,CHIR,UNkT,SNk,GAMI,TPT,CHI,
199 : ! C%C * LIQSOL
200 : ! C%C GOTO 10
201 : ! C%C 112 format(/
202 : ! C%C * ' rho [g/cc] T6 [K] P [Mbar] P/(n_i kT) Cv/(N k)',
203 : ! C%C * ' chi_T chi_r U/(N k T) S/(N k) Gamma_i',
204 : ! C%C * ' T_p/T chi_e liq/sol')
205 : ! C%C 111 format(1P,12E12.3,I2)
206 : ! C%C 113 format(I3,2F8.3,1PE12.4)
207 : ! C%C 114 format(' Z CMI x_j')
208 : ! C%C end
209 :
210 0 : subroutine MELANGE9(NMIX,AY,AZion,ACMI,RHO,TEMP, &
211 : GAMIlo,GAMIhi,PRADnkT, &
212 : DENS,Zmean,CMImean,Z2mean,GAMImean,CHI,TPT,LIQSOL, &
213 : PnkT,UNkT,SNk,CV,CHIR,CHIT,ierr)
214 : ! Version 20.05.13
215 : ! slight optimization 10.12.14
216 : ! Wigner-Kirkwood correction for the entropy corrected 28.05.15
217 : ! Stems from MELANGE8 v.26.12.09.
218 : ! Difference: output PRADnkT instead of input KRAD
219 : ! + EOS of fully ionized electron-ion plasma mixture.
220 : ! Limitations:
221 : ! (a) inapplicable in the regimes of
222 : ! (1) bound-state formation,
223 : ! (2) quantum liquid,
224 : ! (3) presence of positrons;
225 : ! (b) for the case of a composition gradually depending on RHO or TEMP,
226 : ! second-order functions (CV,CHIR,CHIT in output) should not be trusted
227 : ! Choice of the liquid or solid regime - criterion GAMI [because the
228 : ! choice based on comparison of total (non-OCP) free energies can be
229 : ! sometimes dangerous because of the fit uncertainties ("Local field
230 : ! correction" in solid and quantum effects in liquid are unknown)].
231 : ! Input: NMIX - number of different elements;
232 : ! AY - their partial number densities,
233 : ! AZion and ACMI - their charge and mass numbers,
234 : ! RHO - total mass density [g/cc]
235 : ! TEMP - temperature [in a.u.=2Ryd=3.1577e5 K].
236 : ! GAMIlo - begin mixing liquid and solid solutions
237 : ! GAMIhi - phase transition complete into fully solid
238 : ! NB: instead of RHO, a true input is CHI, defined below
239 : ! Hence, disagreement between RHO and DENS is the fit error (<0.4%)
240 : ! Output:
241 : ! AY - rescaled so that to sum up to 1
242 : ! DENS - electron number density [in a.u.=6.7483346e24 cm^{-3}]
243 : ! Zmean=<Z>, CMImean=<A> - mean ion charge and mass numbers,
244 : ! Z2mean=<Z^2> - mean-square ion charge number
245 : ! GAMImean - effective ion-ion Coulomb coupling constant
246 : ! CHI = mu_e/kT, where mu_e is the electron chem.potential
247 : ! TPT - effective ionic quantum parameter (T_p/T)
248 : ! LIQSOL=0/1 for liquid/solid
249 : ! SNk - dimensionless entropy per 1 ion
250 : ! UNkT - internal energy per kT per ion
251 : ! PnkT - pressure / n_i kT, where n_i is the ion number density
252 : ! PRADnkT - radiative pressure / n_i kT
253 : ! CV - heat capacity per ion, div. by Boltzmann const.
254 : ! CHIR - inverse compressibility -(d ln P / d ln V)_T ("\chi_r")
255 : ! CHIT = (d ln P / d ln T)_V ("\chi_T")
256 : integer, intent(in) :: NMIX
257 : integer, intent(out) :: LIQSOL, ierr
258 : real(dp), intent(in) :: AZion(:), ACMI(:), GAMIlo, GAMIhi
259 : real(dp), intent(inout) :: AY(:)
260 : type(auto_diff_real_2var_order1), intent(in) :: RHO, TEMP
261 : type(auto_diff_real_2var_order1), intent(out) :: PRADnkT, DENS, GAMImean
262 : real(dp), intent(out) :: Zmean, CMImean, Z2mean
263 : type(auto_diff_real_2var_order1), intent(out) :: CHI, TPT, PnkT, UNkT, SNk, CV, CHIR, CHIT
264 :
265 : integer :: IX, I, J
266 0 : real(dp) :: Y, Z13, Z52, Z53, Z73, Z321, Zion, CMI
267 : type(auto_diff_real_2var_order1) :: UINTRAD, PRESSRAD
268 : type(auto_diff_real_2var_order1) :: FEid,PEid,UEid,SEid,CVE,CHITE,CHIRE
269 : type(auto_diff_real_2var_order1) :: DlnDH,DlnDT,DlnDHH,DlnDTT,DlnDHT
270 : type(auto_diff_real_2var_order1) :: DTE, PRESSE, UINTE, RS, RSI
271 : type(auto_diff_real_2var_order1) :: GAME, alfa, beta, TPT2
272 : type(auto_diff_real_2var_order1) :: UINT, PRESS, CVtot, Stot, PDLT, PDLR
273 : type(auto_diff_real_2var_order1) :: DENSI, PRESSI, GAMI, DNI, PRI
274 : type(auto_diff_real_2var_order1) :: FC1,UC1,PC1,SC1,CV1,PDT1,PDR1
275 : type(auto_diff_real_2var_order1) :: FC2,UC2,PC2,SC2,CV2,PDT2,PDR2
276 : type(auto_diff_real_2var_order1) :: FC1_0,UC1_0,PC1_0,SC1_0,CV1_0,PDT1_0,PDR1_0
277 : type(auto_diff_real_2var_order1) :: FC2_0,UC2_0,PC2_0,SC2_0,CV2_0,PDT2_0,PDR2_0
278 : type(auto_diff_real_2var_order1) :: FC1_1,UC1_1,PC1_1,SC1_1,CV1_1,PDT1_1,PDR1_1
279 : type(auto_diff_real_2var_order1) :: FC2_1,UC2_1,PC2_1,SC2_1,CV2_1,PDT2_1,PDR2_1
280 : type(auto_diff_real_2var_order1) :: FMIX,UMIX,PMIX,CVMIX,PDTMIX,PDRMIX
281 :
282 :
283 : type(auto_diff_real_2var_order1) :: CTP, FWK, UWK
284 : type(auto_diff_real_2var_order1) :: X, X1, X2, RZ, DeltaG
285 :
286 : real(dp), parameter :: TINY=1.d-7
287 : real(dp), parameter :: AUM=1822.888d0 ! a.m.u./m_e
288 : real(dp), parameter :: GAMIMELT=175.d0 ! OCP value of Gamma_i for melting (not used, replaced by GAMIlo/hi)
289 : real(dp), parameter :: RSIMELT=140.d0 ! ion density parameter of quantum melting
290 : real(dp), parameter :: RAD=2.5568570411948021d-07 ! Radiation constant (=4\sigma/c) (in a.u.)
291 :
292 :
293 0 : ierr = 0
294 0 : Y=0.d0
295 0 : do IX=1,NMIX
296 0 : Y=Y+AY(IX)
297 : end do
298 0 : if (abs(Y-1.d0)>TINY) then
299 0 : do IX=1,NMIX
300 0 : AY(IX)=AY(IX)/Y
301 : end do
302 : ! print*,'MELANGE9: partial densities (and derivatives)',
303 : ! * ' are rescaled by factor',1./Y
304 : end if
305 0 : Zmean=0d0
306 0 : Z2mean=0d0
307 0 : Z52=0d0
308 0 : Z53=0d0
309 0 : Z73=0d0
310 0 : Z321=0d0 ! corr.26.12.09
311 0 : CMImean=0d0
312 0 : do IX=1,NMIX
313 0 : if (AY(IX) < TINY) cycle
314 0 : Zmean=Zmean+AY(IX)*AZion(IX)
315 0 : Z2mean=Z2mean+AY(IX)*AZion(IX)*AZion(IX)
316 0 : Z13=pow(AZion(IX),1d0/3d0)
317 0 : Z53=Z53+AY(IX)*pow5(Z13)
318 0 : Z73=Z73+AY(IX)*pow7(Z13)
319 0 : Z52=Z52+AY(IX)*pow5(sqrt(AZion(IX)))
320 0 : Z321=Z321+AY(IX)*AZion(IX)*pow3(sqrt(AZion(IX)+1.d0)) ! 26.12.09
321 0 : CMImean=CMImean+AY(IX)*ACMI(IX)
322 : end do
323 : ! (0) Photons:
324 0 : UINTRAD=RAD*TEMP*TEMP*TEMP*TEMP
325 0 : PRESSRAD=UINTRAD/3d0
326 : ! CVRAD=4.*UINTRAD/TEMP
327 : ! (1) ideal electron gas (including relativity and degeneracy) ----- *
328 0 : DENS=RHO/11.20587d0*Zmean/CMImean ! number density of electrons [au]
329 0 : call CHEMFIT(DENS,TEMP,CHI)
330 :
331 : ! NB: CHI can be used as true input instead of RHO or DENS
332 : call ELECT11(TEMP,CHI, &
333 : DENS,FEid,PEid,UEid,SEid,CVE,CHITE,CHIRE, &
334 0 : DlnDH,DlnDT,DlnDHH,DlnDTT,DlnDHT)
335 : ! NB: at this point DENS is redefined (the difference can be ~0.1%)
336 0 : DTE=DENS*TEMP
337 0 : PRESSE=PEid*DTE ! P_e [a.u.]
338 0 : UINTE=UEid*DTE ! U_e / V [a.u.]
339 : ! (2) non-ideal Coulomb EIP ---------------------------------------- *
340 0 : RS=pow(0.75d0/PI/DENS,1d0/3d0) ! r_s - electron density parameter
341 0 : RSI=RS*CMImean*Z73*AUM ! R_S - ion density parameter
342 0 : GAME=1d0/RS/TEMP ! electron Coulomb parameter Gamma_e
343 0 : GAMImean=Z53*GAME ! effective Gamma_i - ion Coulomb parameter
344 0 : if (RSI<RSIMELT) then ! doesn't happen in "typical" wd cases
345 0 : LIQSOL=0 ! liquid regime
346 0 : else if (GAMImean<GAMIlo) then
347 0 : LIQSOL=0 ! liquid regime
348 0 : else if (GAMImean>GAMIhi) then
349 0 : LIQSOL=1 ! solid regime
350 : else ! blend of liquid and solid
351 0 : LIQSOL=-1
352 0 : alfa = (GAMImean - GAMIlo)/(GAMIhi - GAMIlo) ! 1 for solid, 0 for liquid
353 0 : beta = 1d0 - alfa
354 : end if
355 : ! Calculate partial thermodynamic quantities and combine them together:
356 0 : UINT=UINTE
357 0 : PRESS=PRESSE
358 0 : CVtot=CVE*DENS
359 0 : Stot=SEid*DENS
360 0 : PDLT=PRESSE*CHITE ! d P_e[a.u.] / d ln T
361 0 : PDLR=PRESSE*CHIRE ! d P_e[a.u.] / d ln\rho
362 0 : DENSI=DENS/Zmean ! number density of all ions
363 0 : PRESSI=DENSI*TEMP ! ideal-ions total pressure (normalization)
364 0 : TPT2=0d0
365 0 : CTP=4.d0*PI/AUM/(TEMP*TEMP) ! common coefficient for TPT2.10.12.14
366 : ! Add Coulomb+xc nonideal contributions, and ideal free energy:
367 0 : do IX=1,NMIX
368 0 : if (AY(IX)<TINY) cycle ! skip this species
369 0 : Zion=AZion(IX)
370 0 : if (Zion==0d0) cycle ! skip neutrons
371 0 : CMI=ACMI(IX)
372 0 : GAMI=pow(Zion,5d0/3d0)*GAME ! Gamma_i for given ion species
373 0 : DNI=DENSI*AY(IX) ! number density of ions of given type
374 0 : PRI=DNI*TEMP ! = ideal-ions partial pressure (normalization)
375 0 : if (LIQSOL == 0 .or. LIQSOL == 1) then
376 : call EOSFI8(LIQSOL,CMI,Zion,RS,GAMI, &
377 : FC1,UC1,PC1,SC1,CV1,PDT1,PDR1, &
378 0 : FC2,UC2,PC2,SC2,CV2,PDT2,PDR2,ierr)
379 0 : if (ierr /= 0) return
380 : else
381 : call EOSFI8(0,CMI,Zion,RS,GAMI, &
382 : FC1_0,UC1_0,PC1_0,SC1_0,CV1_0,PDT1_0,PDR1_0, &
383 0 : FC2_0,UC2_0,PC2_0,SC2_0,CV2_0,PDT2_0,PDR2_0,ierr)
384 0 : if (ierr /= 0) return
385 : call EOSFI8(1,CMI,Zion,RS,GAMI, &
386 : FC1_1,UC1_1,PC1_1,SC1_1,CV1_1,PDT1_1,PDR1_1, &
387 0 : FC2_1,UC2_1,PC2_1,SC2_1,CV2_1,PDT2_1,PDR2_1,ierr)
388 0 : if (ierr /= 0) return
389 0 : FC1 = alfa*FC1_1 + beta*FC1_0
390 0 : UC1 = alfa*UC1_1 + beta*UC1_0
391 0 : PC1 = alfa*PC1_1 + beta*PC1_0
392 0 : SC1 = alfa*SC1_1 + beta*SC1_0
393 0 : CV1 = alfa*CV1_1 + beta*CV1_0
394 0 : PDT1 = alfa*PDT1_1 + beta*PDT1_0
395 0 : PDR1 = alfa*PDR1_1 + beta*PDR1_0
396 0 : FC2 = alfa*FC2_1 + beta*FC2_0
397 0 : UC2 = alfa*UC2_1 + beta*UC2_0
398 0 : PC2 = alfa*PC2_1 + beta*PC2_0
399 0 : SC2 = alfa*SC2_1 + beta*SC2_0
400 0 : CV2 = alfa*CV2_1 + beta*CV2_0
401 0 : PDT2 = alfa*PDT2_1 + beta*PDT2_0
402 0 : PDR2 = alfa*PDR2_1 + beta*PDR2_0
403 : end if
404 : ! First-order TD functions:
405 0 : UINT=UINT+UC2*PRI ! internal energy density (e+i+Coul.)
406 0 : Stot=Stot+DNI*(SC2-log(AY(IX))) !entropy per unit volume[a.u.]
407 0 : PRESS=PRESS+PC2*PRI ! pressure (e+i+Coul.) [a.u.]
408 : ! Second-order functions (they take into account compositional changes):
409 0 : CVtot=CVtot+DNI*CV2 ! C_V (e+i+Coul.)/ V (optim.10.12.14)
410 0 : PDLT=PDLT+PRI*PDT2 ! d P / d ln T
411 0 : PDLR=PDLR+PRI*PDR2 ! d P / d ln\rho
412 0 : TPT2=TPT2+CTP*DNI/ACMI(IX)*AZion(IX)**2 ! opt.10.12.14
413 : end do ! next IX
414 : ! Wigner-Kirkwood perturbative correction for liquid:
415 0 : TPT=sqrt(TPT2) ! effective T_p/T - ion quantum parameter
416 : ! (in the case of a mixture, this estimate is crude)
417 0 : if (LIQSOL==0) then
418 0 : FWK=TPT2/24.d0 ! Wigner-Kirkwood (quantum diffr.) term
419 : ! MESA doesn't warn/error when this term gets large because
420 : ! it is not clear that we are better off falling back to HELM
421 : !
422 : ! if (FWK > .7d0) then
423 : ! print*,'MELANGE9: strong quantum effects in liquid!'
424 : ! ierr = -1
425 : ! return
426 : ! end if
427 0 : UWK=2.d0*FWK
428 0 : UINT=UINT+UWK*PRESSI
429 0 : Stot=Stot+FWK*DENSI ! corrected 28.05.15
430 0 : PRESS=PRESS+FWK*PRESSI
431 0 : CVtot=CVtot-UWK*DENSI ! corrected by JWS 17.04.20
432 0 : PDLT=PDLT-FWK*PRESSI
433 0 : PDLR=PDLR+UWK*PRESSI
434 : end if
435 : ! Corrections to the linear mixing rule:
436 0 : if (LIQSOL==0) then ! liquid phase
437 : call CORMIX(RS,GAME,Zmean,Z2mean,Z52,Z53,Z321, &
438 0 : FMIX,UMIX,PMIX,CVMIX,PDTMIX,PDRMIX)
439 : else ! solid phase (only Madelung contribution) [22.12.12]
440 0 : FMIX=0d0
441 0 : do I=1,NMIX
442 0 : if (AY(I)<TINY) cycle
443 0 : do J=I+1,NMIX
444 0 : if (AY(J)<TINY) cycle
445 0 : RZ=AZion(J)/AZion(I)
446 0 : X2=AY(J)/(AY(I)+AY(J))
447 0 : X1=dim(1.d0,X2)
448 0 : if (X2<TINY) cycle
449 0 : X=X2/RZ+(1.d0-1.d0/RZ)*pow(X2, RZ)
450 0 : GAMI=pow(AZion(I),5d0/3d0)*GAME ! Gamma_i corrected 14.05.13
451 0 : DeltaG=.012d0*(1.d0-1.d0/pow2(RZ))*(X1+X2*pow(RZ,5d0/3d0))
452 0 : DeltaG=DeltaG*X/X2*dim(1.d0,X)/X1
453 0 : FMIX=FMIX+AY(I)*AY(J)*GAMI*DeltaG
454 : end do
455 : end do
456 0 : UMIX=FMIX
457 0 : PMIX=FMIX/3.d0
458 0 : CVMIX=0d0
459 0 : PDTMIX=0d0
460 0 : PDRMIX=FMIX/2.25d0
461 : end if
462 0 : UINT=UINT+UMIX*PRESSI
463 0 : Stot=Stot+DENSI*(UMIX-FMIX)
464 0 : PRESS=PRESS+PMIX*PRESSI
465 0 : CVtot=CVtot+DENSI*CVMIX
466 0 : PDLT=PDLT+PRESSI*PDTMIX
467 0 : PDLR=PDLR+PRESSI*PDRMIX
468 : ! First-order:
469 0 : PRADnkT=PRESSRAD/PRESSI ! radiative pressure / n_i k T
470 : ! CVtot=CVtot+CVRAD
471 : ! Stot=Stot+CVRAD/3.
472 0 : PnkT=PRESS/PRESSI ! P / n_i k T
473 0 : UNkT=UINT/PRESSI ! U / N_i k T
474 : ! UNkT=UNkT+UINTRAD/PRESSI
475 0 : SNk=Stot/DENSI ! S / N_i k
476 : ! Second-order:
477 0 : CV=CVtot/DENSI ! C_V per ion
478 0 : CHIR=PDLR/PRESS ! d ln P / d ln\rho
479 0 : CHIT=PDLT/PRESS ! d ln P / d ln T
480 : ! CHIT=CHIT+4.*PRESSRAD/PRESS ! d ln P / d ln T
481 0 : return
482 : end subroutine MELANGE9
483 :
484 0 : subroutine EOSFI8(LIQSOL,CMI,Zion,RS,GAMI, &
485 : FC1,UC1,PC1,SC1,CV1,PDT1,PDR1, &
486 : FC2,UC2,PC2,SC2,CV2,PDT2,PDR2,ierr)
487 : ! Version 16.09.08
488 : ! call FHARM8 has been replaced by call FHARM12 27.04.12
489 : ! Wigner-Kirkwood correction excluded 20.05.13
490 : ! slight cleaning 10.12.14
491 : ! Non-ideal parts of thermodynamic functions in the fully ionized plasma
492 : ! Stems from EOSFI5 and EOSFI05 v.04.10.05
493 : ! Input: LIQSOL=0/1(liquid/solid),
494 : ! Zion,CMI - ion charge and mass numbers,
495 : ! RS=r_s (electronic density parameter),
496 : ! GAMI=Gamma_i (ion coupling),
497 : ! Output: FC1 and UC1 - non-ideal "ii+ie+ee" contribution to the
498 : ! free and internal energies (per ion per kT),
499 : ! PC1 - analogous contribution to pressure divided by (n_i kT),
500 : ! CV1 - "ii+ie+ee" heat capacity per ion [units of k]
501 : ! PDT1=(1/n_i kT)*(d P_C/d ln T)_V
502 : ! PDR1=(1/n_i kT)*(d P_C/d ln\rho)_T
503 : ! FC2,UC2,PC2,SC2,CV2 - analogous to FC1,UC1,PC1,SC1,CV1, but including
504 : ! the part corresponding to the ideal ion gas. This is useful for
505 : ! preventing accuracy loss in some cases (e.g., when SC2 << SC1).
506 : ! FC2 does not take into account the entropy of mixing S_{mix}: in a
507 : ! mixture, S_{mix}/(N_i k) has to be added externally (see MELANGE9).
508 : ! FC2 does not take into account the ion spin degeneracy either.
509 : ! When needed, the spin term must be added to the entropy externally.
510 : integer, intent(in) :: LIQSOL
511 : real(dp), intent(in) :: CMI, Zion
512 : type(auto_diff_real_2var_order1), intent(in) :: RS, GAMI
513 : type(auto_diff_real_2var_order1), intent(out) :: FC1,UC1,PC1,SC1,CV1,PDT1,PDR1
514 : type(auto_diff_real_2var_order1), intent(out) :: FC2,UC2,PC2,SC2,CV2,PDT2,PDR2
515 : integer, intent(out) :: ierr
516 :
517 : type(auto_diff_real_2var_order1) :: GAME,FXC,UXC,PXC,CVXC,SXC,PDTXC,PDRXC
518 : type(auto_diff_real_2var_order1) :: TPT, COTPT, FidION
519 : type(auto_diff_real_2var_order1) :: FION, UION, PION, CVii, PDTii, PDRii
520 : type(auto_diff_real_2var_order1) :: FItot, UItot, PItot, CVItot, SCItot
521 : type(auto_diff_real_2var_order1) :: PDTi, PDRi
522 : type(auto_diff_real_2var_order1) :: Fharm,Uharm,Pharm,CVharm,Sharm,PDTharm,PDRharm
523 : type(auto_diff_real_2var_order1) :: Fah,Uah,Pah,CVah,PDTah,PDRah
524 : type(auto_diff_real_2var_order1) :: FSCR,USCR,PSCR,S_SCR,CVSCR,PDTSCR,PDRSCR
525 : type(auto_diff_real_2var_order1) :: FC0, UC0, PC0, SC0, CV0, PDT0, PDR0
526 :
527 : real(dp), parameter :: TINY=1.d-20
528 : real(dp), parameter :: AUM=1822.888d0 ! a.m.u/m_e
529 :
530 0 : ierr = 0
531 0 : if (LIQSOL/=1.and.LIQSOL/=0) then
532 0 : ierr = -1
533 0 : return
534 : !call mesa_error(__FILE__,__LINE__,'EOSFI8: invalid LIQSOL')
535 : end if
536 0 : if (CMI<=.1d0) then
537 0 : ierr = -1
538 0 : return
539 : !call mesa_error(__FILE__,__LINE__,'EOSFI8: too small CMI')
540 : end if
541 0 : if (Zion<=.1d0) then
542 0 : ierr = -1
543 0 : return
544 : !call mesa_error(__FILE__,__LINE__,'EOSFI8: too small Zion')
545 : end if
546 0 : if (RS<=.0d0) then
547 0 : ierr = -1
548 0 : return
549 : !call mesa_error(__FILE__,__LINE__,'EOSFI8: invalid RS')
550 : end if
551 0 : if (GAMI<=.0d0) then
552 0 : ierr = -1
553 0 : return
554 : !call mesa_error(__FILE__,__LINE__,'EOSFI8: invalid GAMI')
555 : end if
556 0 : GAME=GAMI/pow(Zion,5d0/3d0)
557 0 : call EXCOR7(RS,GAME,FXC,UXC,PXC,CVXC,SXC,PDTXC,PDRXC) ! "ee"("xc")
558 : ! Calculate "ii" part:
559 0 : COTPT=sqrt(3d0/AUM/CMI)/pow(Zion,7d0/6d0) ! auxiliary coefficient
560 0 : TPT=GAMI/sqrt(RS)*COTPT ! T_p/T
561 0 : FidION=1.5d0*log(TPT*TPT/GAMI)-1.323515d0
562 : ! 1.3235=1+0.5*ln(6/pi); FidION = F_{id.ion gas}/(N_i kT), but without
563 : ! the term x_i ln x_i = -S_{mix}/(N_i k).
564 0 : if (LIQSOL==0) then ! liquid
565 : call FITION9(GAMI, &
566 0 : FION,UION,PION,CVii,PDTii,PDRii)
567 0 : FItot=FION+FidION
568 0 : UItot=UION+1.5d0
569 0 : PItot=PION+1.0d0
570 0 : CVItot=CVii+1.5d0
571 0 : SCItot=UItot-FItot
572 0 : PDTi=PDTii+1.d0
573 0 : PDRi=PDRii+1.d0
574 : else ! solid
575 : call FHARM12(GAMI,TPT, &
576 0 : Fharm,Uharm,Pharm,CVharm,Sharm,PDTharm,PDRharm) ! harm."ii"
577 0 : call ANHARM8(GAMI,TPT,Fah,Uah,Pah,CVah,PDTah,PDRah) ! anharm.
578 0 : FItot=Fharm+Fah
579 0 : FION=FItot-FidION
580 0 : UItot=Uharm+Uah
581 0 : UION=UItot-1.5d0 ! minus 1.5=ideal-gas, in order to get "ii"
582 0 : PItot=Pharm+Pah
583 0 : PION=PItot-1.d0 ! minus 1=ideal-gas
584 0 : PDTi=PDTharm+PDTah
585 0 : PDRi=PDRharm+PDRah
586 0 : PDTii=PDTi-1.d0 ! minus 1=ideal-gas
587 0 : PDRii=PDRi-1.d0 ! minus 1=ideal-gas
588 0 : CVItot=CVharm+CVah
589 0 : SCItot=Sharm+Uah-Fah
590 0 : CVii=CVItot-1.5d0 ! minus 1.5=ideal-gas
591 : end if
592 : ! Calculate "ie" part:
593 0 : if (LIQSOL==1) then
594 : call FSCRsol8(RS,GAMI,Zion,TPT, &
595 0 : FSCR,USCR,PSCR,S_SCR,CVSCR,PDTSCR,PDRSCR,ierr)
596 0 : if (ierr /= 0) return
597 : else
598 : call FSCRliq8(RS,GAME,Zion, &
599 0 : FSCR,USCR,PSCR,CVSCR,PDTSCR,PDRSCR,ierr)
600 0 : if (ierr /= 0) return
601 0 : S_SCR=USCR-FSCR
602 : end if
603 : ! Total excess quantities ("ii"+"ie"+"ee", per ion):
604 0 : FC0=FSCR+Zion*FXC
605 0 : UC0=USCR+Zion*UXC
606 0 : PC0=PSCR+Zion*PXC
607 0 : SC0=S_SCR+Zion*SXC
608 0 : CV0=CVSCR+Zion*CVXC
609 0 : PDT0=PDTSCR+Zion*PDTXC
610 0 : PDR0=PDRSCR+Zion*PDRXC
611 0 : FC1=FION+FC0
612 0 : UC1=UION+UC0
613 0 : PC1=PION+PC0
614 0 : SC1=(UION-FION)+SC0
615 0 : CV1=CVii+CV0
616 0 : PDT1=PDTii+PDT0
617 0 : PDR1=PDRii+PDR0
618 : ! Total excess + ideal-ion quantities
619 0 : FC2=FItot+FC0
620 0 : UC2=UItot+UC0
621 0 : PC2=PItot+PC0
622 0 : SC2=SCItot+SC0
623 0 : CV2=CVItot+CV0
624 0 : PDT2=PDTi+PDT0
625 0 : PDR2=PDRi+PDR0
626 0 : return
627 : end subroutine EOSFI8
628 :
629 : ! ================== ELECTRON-ION COULOMB LIQUID =================== *
630 0 : subroutine FITION9(GAMI, &
631 : FION,UION,PION,CVii,PDTii,PDRii)
632 : ! Version 11.09.08
633 : ! Non-ideal contributions to thermodynamic functions of classical OCP,
634 : ! corrected at small density for a mixture.
635 : ! Stems from FITION00 v.24.05.00.
636 : ! Input: GAMI - ion coupling parameter
637 : ! Output: FION - ii free energy / N_i kT
638 : ! UION - ii internal energy / N_i kT
639 : ! PION - ii pressure / n_i kT
640 : ! CVii - ii heat capacity / N_i k
641 : ! PDTii = PION + d(PION)/d ln T = (1/N_i kT)*(d P_{ii}/d ln T)
642 : ! PDRii = PION + d(PION)/d ln\rho
643 : ! Parameters adjusted to Caillol (1999).
644 : use pc_support
645 : type(auto_diff_real_2var_order1), intent(in) :: GAMI
646 : type(auto_diff_real_2var_order1), intent(out) :: FION, UION, PION, CVii, PDTii, PDRii
647 : type(auto_diff_real_2var_order1) :: F0, U0
648 :
649 : real(dp), parameter :: A1=-.907347d0
650 : real(dp), parameter :: A2=.62849d0
651 : real(dp) :: A3
652 : real(dp), parameter :: C1=.004500d0
653 : real(dp), parameter :: G1=170.0d0
654 : real(dp), parameter :: C2=-8.4d-5
655 : real(dp), parameter :: G2=.0037d0
656 : real(dp), parameter :: SQ32=.8660254038d0 ! SQ32=sqrt(3)/2
657 :
658 : real(dp) :: &
659 0 : xFION, dFION_dlnGAMI, &
660 0 : xUION, dUION_dlnGAMI, &
661 0 : xPION, dPION_dlnGAMI, &
662 0 : xCVii, dCVii_dlnGAMI, &
663 0 : xPDTii, dPDTii_dlnGAMI, &
664 0 : xPDRii, dPDRii_dlnGAMI
665 0 : real(dp) :: dlnGAMI_dT, dlnGAMI_dRho
666 : integer :: ierr
667 : logical :: skip
668 : logical, parameter :: use_FITION9_table = .false.
669 : logical, parameter :: debug_FITION9_table = .false.
670 :
671 :
672 : if (.not. use_FITION9_table .or. debug_FITION9_table) then
673 :
674 0 : A3=-SQ32-A1/sqrt(A2)
675 : F0=A1*(sqrt(GAMI*(A2+GAMI)) &
676 : - A2*log(sqrt(GAMI/A2)+sqrt(1d0+GAMI/A2))) &
677 0 : + 2d0*A3*(sqrt(GAMI)-atan(sqrt(GAMI)))
678 0 : U0=pow3(sqrt(GAMI))*(A1/sqrt(A2+GAMI)+A3/(1.d0+GAMI))
679 : ! This is the zeroth approximation. Correction:
680 0 : UION=U0+C1*GAMI*GAMI/(G1+GAMI)+C2*GAMI*GAMI/(G2+GAMI*GAMI)
681 : FION=F0+C1*(GAMI-G1*log(1.d0+GAMI/G1)) &
682 0 : + C2/2d0*log(1.d0+GAMI*GAMI/G2)
683 : CVii=-0.5d0*pow3(sqrt(GAMI))*(A1*A2/pow3(sqrt(A2+GAMI)) &
684 : + A3*(1.d0-GAMI)/pow2(1.d0+GAMI)) &
685 0 : - GAMI*GAMI*(C1*G1/pow2(G1+GAMI)+C2*(G2-GAMI*GAMI)/pow2(G2+GAMI*GAMI))
686 0 : PION=UION/3.0d0
687 0 : PDRii=(4.0d0*UION-CVii)/9.0d0 ! p_{ii} + d p_{ii} / d ln\rho
688 0 : PDTii=CVii/3.0d0 ! p_{ii} + d p_{ii} / d ln T
689 :
690 : end if
691 :
692 : if (use_FITION9_table .or. debug_FITION9_table) then
693 : ierr = 0
694 : call get_FITION9(GAMI%val, &
695 : xFION, dFION_dlnGAMI, &
696 : xUION, dUION_dlnGAMI, &
697 : xPION, dPION_dlnGAMI, &
698 : xCVii, dCVii_dlnGAMI, &
699 : xPDTii, dPDTii_dlnGAMI, &
700 : xPDRii, dPDRii_dlnGAMI, &
701 : skip, ierr)
702 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'failed in call get_FITION9')
703 : else
704 0 : skip = .true.
705 : end if
706 :
707 :
708 : if (debug_FITION9_table .and. .not. skip) then
709 :
710 : if (.not. check1(FION, xFION, 'FION')) return
711 : if (.not. check1(UION, xUION, 'UION')) return
712 : if (.not. check1(PION, xPION, 'PION')) return
713 : if (.not. check1(CVii, xCVii, 'CVii')) return
714 : if (.not. check1(PDTii, xPDTii, 'PDTii')) return
715 : if (.not. check1(PDRii, xPDRii, 'PDRii')) return
716 :
717 : end if
718 :
719 :
720 0 : if (use_FITION9_table .and. .not. skip) then
721 :
722 : ! values
723 : FION% val = xFION
724 : UION% val = xUION
725 : PION% val = xPION
726 : CVii% val = xCVii
727 : PDTii% val = xPDTii
728 : PDRii% val = xPDRii
729 :
730 : ! dT (val1) derivatives
731 : ! dlnRs_dT = RS% d1val1 / RS% val = 0
732 : dlnGAMI_dT = GAMI% d1val1 / GAMI% val
733 :
734 : FION% d1val1 = dFION_dlnGAMI * dlnGAMI_dT
735 : UION% d1val1 = dUION_dlnGAMI * dlnGAMI_dT
736 : PION% d1val1 = dPION_dlnGAMI * dlnGAMI_dT
737 : CVii% d1val1 = dCVii_dlnGAMI * dlnGAMI_dT
738 : PDTii% d1val1 = dPDTii_dlnGAMI * dlnGAMI_dT
739 : PDRii% d1val1 = dPDRii_dlnGAMI * dlnGAMI_dT
740 :
741 : ! dRho (val2) derivatives
742 : dlnGAMI_dRho = GAMI% d1val2 / GAMI% val
743 :
744 : FION% d1val2 = dFION_dlnGAMI * dlnGAMI_dRho
745 : UION% d1val2 = dUION_dlnGAMI * dlnGAMI_dRho
746 : PION% d1val2 = dPION_dlnGAMI * dlnGAMI_dRho
747 : CVii% d1val2 = dCVii_dlnGAMI * dlnGAMI_dRho
748 : PDTii% d1val2 = dPDTii_dlnGAMI * dlnGAMI_dRho
749 : PDRii% d1val2 = dPDRii_dlnGAMI * dlnGAMI_dRho
750 :
751 : end if
752 :
753 :
754 : contains
755 :
756 : logical function check1(v, xv, str)
757 : type(auto_diff_real_2var_order1), intent(in) :: v
758 : real(dp), intent(in) :: xv
759 : character (len=*), intent(in) :: str
760 : real(dp) :: val
761 : real(dp), parameter :: atol = 1d-8, rtol = 1d-6
762 : include 'formats'
763 : val = v%val
764 : check1 = .false.
765 : if (is_bad(xv)) then
766 : write(*,*) 'is_bad ' // trim(str), xv, val, GAMI
767 : return
768 : end if
769 : if (abs(val - xv) > atol + rtol*max(abs(val),abs(xv))) then
770 : write(*,*) 'rel mismatch ' // trim(str), &
771 : (val - xv)/max(abs(val),abs(xv),1d-99), &
772 : xv, val, GAMI%val
773 : call mesa_error(__FILE__,__LINE__,'FITION9')
774 : return
775 : end if
776 : check1 = .true.
777 0 : end function check1
778 :
779 : end subroutine FITION9
780 :
781 0 : subroutine FSCRliq8(RS,GAME,Zion, &
782 : FSCR,USCR,PSCR,CVSCR,PDTSCR,PDRSCR,ierr) ! fit to the el.-ion scr.
783 : ! Version 11.09.08
784 : ! cleaned 16.06.09
785 : ! Stems from FSCRliq7 v. 09.06.07. Included a check for RS=0.
786 : ! INPUT: RS - density parameter, GAME - electron Coulomb parameter,
787 : ! Zion - ion charge number,
788 : ! OUTPUT: FSCR - screening (e-i) free energy per kT per 1 ion,
789 : ! USCR - internal energy per kT per 1 ion (screen.contrib.)
790 : ! PSCR - pressure divided by (n_i kT) (screen.contrib.)
791 : ! CVSCR - heat capacity per 1 ion (screen.contrib.)
792 : ! PDTSCR,PDRSCR = PSCR + d PSCR / d ln(T,\rho)
793 : use pc_support
794 : type(auto_diff_real_2var_order1), intent(in) :: RS,GAME
795 : real(dp), intent(in) :: Zion
796 : type(auto_diff_real_2var_order1), intent(out) :: FSCR,USCR,PSCR,CVSCR,PDTSCR,PDRSCR
797 : integer, intent(out) :: ierr
798 :
799 : type(auto_diff_real_2var_order1) :: SQG, SQR, SQZ1, SQZ, CDH0, CDH, ZLN, Z13
800 : type(auto_diff_real_2var_order1) :: X, CTF, P01, P03, PTX
801 : type(auto_diff_real_2var_order1) :: TX, TXDG, TXDGG, TY1, TY1DG, TY1DGG
802 : type(auto_diff_real_2var_order1) :: TY2, TY2DX, TY2DXX, TY, TYX, TYDX, TYDG, P1
803 : type(auto_diff_real_2var_order1) :: COR1, COR1DX, COR1DG, COR1DXX, COR1DGG, COR1DXG
804 : type(auto_diff_real_2var_order1) :: U0, U0DX, U0DG, U0DXX, U0DGG, U0DXG
805 : type(auto_diff_real_2var_order1) :: D0DG, D0, D0DX, D0DXX
806 : type(auto_diff_real_2var_order1) :: COR0, COR0DX, COR0DG, COR0DXX, COR0DGG, COR0DXG
807 : type(auto_diff_real_2var_order1) :: RELE, Q1, Q2, H1U, H1D, H1, H1X, H1DX, H1DXX
808 : type(auto_diff_real_2var_order1) :: UP, UPDX, UPDG, UPDXX, UPDGG, UPDXG
809 : type(auto_diff_real_2var_order1) :: DN1, DN1DX, DN1DG, DN1DXX, DN1DGG, DN1DXG
810 : type(auto_diff_real_2var_order1) :: DN, DNDX, DNDG, DNDXX, DNDGG, DNDXG
811 : type(auto_diff_real_2var_order1) :: FX, FXDG, FDX, FG, FDG, FDGDH, FDXX, FDGG, FDXG
812 :
813 : real(dp), parameter :: XRS=.0140047d0
814 : real(dp), parameter :: TINY=1.d-19
815 :
816 : real(dp) :: &
817 0 : xFSCR, dFSCR_dlnRS, dFSCR_dlnGAME, &
818 0 : xUSCR, dUSCR_dlnRS, dUSCR_dlnGAME, &
819 0 : xPSCR, dPSCR_dlnRS, dPSCR_dlnGAME, &
820 0 : xCVSCR, dCVSCR_dlnRS, dCVSCR_dlnGAME, &
821 0 : xPDTSCR, dPDTSCR_dlnRS, dPDTSCR_dlnGAME, &
822 0 : xPDRSCR, dPDRSCR_dlnRS, dPDRSCR_dlnGAME
823 0 : real(dp) :: dlnRs_dT, dlnRs_dRho, dlnGAME_dT, dlnGAME_dRho
824 : logical :: skip
825 : logical, parameter :: use_FSCRliq8_table = .false.
826 : logical, parameter :: debug_FSCRliq8_table = .false.
827 :
828 0 : ierr = 0
829 0 : if (RS<0d0) then
830 0 : ierr = -1
831 0 : return
832 : !call mesa_error(__FILE__,__LINE__,'FSCRliq8: RS < 0')
833 : end if
834 0 : if (RS<TINY) then
835 0 : FSCR=0.d0
836 0 : USCR=0.d0
837 0 : PSCR=0.d0
838 0 : CVSCR=0.d0
839 0 : PDTSCR=0.d0
840 0 : PDRSCR=0.d0
841 0 : return
842 : end if
843 :
844 : if (use_FSCRliq8_table .or. debug_FSCRliq8_table) then
845 : ierr = 0
846 : call get_FSCRliq8(int(Zion), RS%val, GAME%val, &
847 : xFSCR, dFSCR_dlnRS, dFSCR_dlnGAME, &
848 : xUSCR, dUSCR_dlnRS, dUSCR_dlnGAME, &
849 : xPSCR, dPSCR_dlnRS, dPSCR_dlnGAME, &
850 : xCVSCR, dCVSCR_dlnRS, dCVSCR_dlnGAME, &
851 : xPDTSCR, dPDTSCR_dlnRS, dPDTSCR_dlnGAME, &
852 : xPDRSCR, dPDRSCR_dlnRS, dPDRSCR_dlnGAME, skip, ierr)
853 : if (ierr /= 0) return
854 : else
855 0 : skip = .true.
856 : end if
857 :
858 : if (.not. use_FSCRliq8_table .or. debug_FSCRliq8_table .or. skip) then
859 :
860 0 : SQG=sqrt(GAME)
861 0 : SQR=sqrt(RS)
862 0 : SQZ1=sqrt(1d0+Zion)
863 0 : SQZ=sqrt(Zion)
864 0 : CDH0=Zion/1.73205d0 ! 1.73205=sqrt(3.)
865 0 : CDH=CDH0*(SQZ1*SQZ1*SQZ1-SQZ*SQZ*SQZ-1d0)
866 0 : SQG=sqrt(GAME)
867 0 : ZLN=log(Zion)
868 0 : Z13=exp(ZLN/3.d0) ! Zion**(1./3.)
869 0 : X=XRS/RS ! relativity parameter
870 0 : CTF=Zion*Zion*.2513d0*(Z13-1d0+.2d0/sqrt(Z13))
871 : ! Thomas-Fermi constant; .2513=(18/175)(12/\pi)^{2/3}
872 0 : P01=1.11d0*exp(0.475d0*ZLN)
873 0 : P03=0.2d0+0.078d0*ZLN*ZLN
874 0 : PTX=1.16d0+0.08d0*ZLN
875 0 : TX=pow(GAME,PTX)
876 0 : TXDG=PTX*TX/GAME
877 0 : TXDGG=(PTX-1.d0)*TXDG/GAME
878 0 : TY1=1d0/(1.d-3*Zion*Zion+2d0*GAME)
879 0 : TY1DG=-2d0*TY1*TY1
880 0 : TY1DGG=-4d0*TY1*TY1DG
881 0 : TY2=1d0+6d0*RS*RS
882 0 : TY2DX=-12d0*RS*RS/X
883 0 : TY2DXX=-3d0*TY2DX/X
884 0 : TY=RS*RS*RS/TY2*(1d0+TY1)
885 0 : TYX=3d0/X+TY2DX/TY2
886 0 : TYDX=-TY*TYX
887 0 : TYDG=RS*RS*RS*TY1DG/TY2
888 0 : P1=(Zion-1d0)/9d0
889 0 : COR1=1d0+P1*TY
890 0 : COR1DX=P1*TYDX
891 0 : COR1DG=P1*TYDG
892 0 : COR1DXX=P1*(TY*(3d0/(X*X)+(TY2DX/TY2)*(TY2DX/TY2)-TY2DXX/TY2)-TYDX*TYX)
893 0 : COR1DGG=P1*RS*RS*RS*TY1DGG/TY2
894 0 : COR1DXG=-P1*TYDG*TYX
895 0 : U0=0.78d0*sqrt(GAME/Zion)*RS*RS*RS
896 0 : U0DX=-3d0*U0/X
897 0 : U0DG=.5d0*U0/GAME
898 0 : U0DXX=-4.d0*U0DX/X
899 0 : U0DGG=-.5d0*U0DG/GAME
900 0 : U0DXG=-3.d0*U0DG/X
901 0 : D0DG=Zion*Zion*Zion
902 0 : D0=GAME*D0DG+21d0*RS*RS*RS
903 0 : D0DX=-63d0*RS*RS*RS/X
904 0 : D0DXX=252d0*RS*RS*RS/(X*X)
905 0 : COR0=1d0+U0/D0
906 0 : COR0DX=(U0DX-U0*D0DX/D0)/D0
907 0 : COR0DG=(U0DG-U0*D0DG/D0)/D0
908 0 : COR0DXX=(U0DXX-(2d0*U0DX*D0DX+U0*D0DXX)/D0+2d0*(D0DX/D0)*(D0DX/D0))/D0
909 0 : COR0DGG=(U0DGG-2d0*U0DG*D0DG/D0+2d0*U0*(D0DG/D0)*(D0DG/D0))/D0
910 0 : COR0DXG=(U0DXG-(U0DX*D0DG+U0DG*D0DX)/D0+2d0*U0*D0DX*D0DG/(D0*D0))/D0
911 : ! Relativism:
912 0 : RELE=sqrt(1.d0+X*X)
913 0 : Q1=0.18d0/sqrt(sqrt(Zion))
914 0 : Q2=0.2d0+0.37d0/sqrt(Zion)
915 0 : H1U=1d0+X*X/5.d0
916 0 : H1D=1d0+Q1*X+Q2*X*X
917 0 : H1=H1U/H1D
918 0 : H1X=0.4d0*X/H1U-(Q1+2d0*Q2*X)/H1D
919 0 : H1DX=H1*H1X
920 : H1DXX=H1DX*H1X &
921 0 : + H1*(0.4d0/H1U-(0.4d0*X/H1U)*(0.4d0*X/H1U)-2d0*Q2/H1D+pow2((Q1+2d0*Q2*X)/H1D))
922 0 : UP=CDH*SQG+P01*CTF*TX*COR0*H1
923 0 : UPDX=P01*CTF*TX*(COR0DX*H1+COR0*H1DX)
924 0 : UPDG=.5d0*CDH/SQG+P01*CTF*(TXDG*COR0+TX*COR0DG)*H1
925 0 : UPDXX=P01*CTF*TX*(COR0DXX*H1+2d0*COR0DX*H1DX+COR0*H1DXX)
926 : UPDGG=-.25d0*CDH/(SQG*GAME) &
927 0 : + P01*CTF*(TXDGG*COR0+2d0*TXDG*COR0DG+TX*COR0DGG)*H1
928 : UPDXG=P01*CTF*(TXDG*(COR0DX*H1+COR0*H1DX) &
929 0 : + TX*(COR0DXG*H1+COR0DG*H1DX))
930 0 : DN1=P03*SQG+P01/RS*TX*COR1
931 0 : DN1DX=P01*TX*(COR1/XRS+COR1DX/RS)
932 0 : DN1DG=.5d0*P03/SQG+P01/RS*(TXDG*COR1+TX*COR1DG)
933 0 : DN1DXX=P01*TX/XRS*(2d0*COR1DX+X*COR1DXX)
934 : DN1DGG=-.25d0*P03/(GAME*SQG) &
935 0 : + P01/RS*(TXDGG*COR1+2d0*TXDG*COR1DG+TX*COR1DGG)
936 0 : DN1DXG=P01*(TXDG*(COR1/XRS+COR1DX/RS)+TX*(COR1DG/XRS+COR1DXG/RS))
937 0 : DN=1d0+DN1/RELE
938 0 : DNDX=DN1DX/RELE-X*DN1/(RELE*RELE*RELE)
939 0 : DNDXX=(DN1DXX-((2d0*X*DN1DX+DN1)-3.d0*X*X*DN1/(RELE*RELE))/(RELE*RELE))/RELE
940 0 : DNDG=DN1DG/RELE
941 0 : DNDGG=DN1DGG/RELE
942 0 : DNDXG=DN1DXG/RELE-X*DN1DG/(RELE*RELE*RELE)
943 0 : FSCR=-UP/DN*GAME
944 0 : FX=(UP*DNDX/DN-UPDX)/DN
945 0 : FXDG=((UPDG*DNDX+UPDX*DNDG+UP*DNDXG-2d0*UP*DNDX*DNDG/DN)/DN-UPDXG)/DN
946 0 : FDX=FX*GAME
947 0 : FG=(UP*DNDG/DN-UPDG)/DN
948 0 : FDG=FG*GAME-UP/DN
949 0 : FDGDH=SQG*DNDG/(DN*DN) ! d FDG / d CDH
950 0 : FDXX=((UP*DNDXX+2d0*(UPDX*DNDX-UP*DNDX*DNDX/DN))/DN-UPDXX)/DN*GAME
951 0 : FDGG=2d0*FG+GAME*((2d0*DNDG*(UPDG-UP*DNDG/DN)+UP*DNDGG)/DN-UPDGG)/DN
952 0 : FDXG=FX+GAME*FXDG
953 0 : USCR=GAME*FDG
954 0 : CVSCR=-GAME*GAME*FDGG
955 0 : PSCR=(X*FDX+GAME*FDG)/3.d0
956 0 : PDTSCR=-GAME*GAME*(X*FXDG+FDGG)/3.d0
957 0 : PDRSCR=(12d0*PSCR+X*X*FDXX+2d0*X*GAME*FDXG+GAME*GAME*FDGG)/9d0
958 :
959 : end if
960 :
961 : if (debug_FSCRliq8_table .and. .not. skip) then
962 :
963 : if (.not. check1(FSCR, xFSCR, 'FSCR')) return
964 : if (.not. check1(USCR, xUSCR, 'USCR')) return
965 : if (.not. check1(PSCR, xPSCR, 'PSCR')) return
966 : if (.not. check1(CVSCR, xCVSCR, 'CVSCR')) return
967 : if (.not. check1(PDTSCR, xPDTSCR, 'PDTSCR')) return
968 : if (.not. check1(PDRSCR, xPDRSCR, 'PDRSCR')) return
969 :
970 : end if
971 :
972 0 : if (use_FSCRliq8_table .and. .not. skip) then
973 :
974 : ! values
975 : FSCR% val = xFSCR
976 : USCR% val = xUSCR
977 : PSCR% val = xPSCR
978 : CVSCR% val = xCVSCR
979 : PDTSCR% val = xPDTSCR
980 : PDRSCR% val = xPDRSCR
981 :
982 : ! dT (val1) derivatives
983 : ! dlnRs_dT = RS% d1val1 / RS% val = 0
984 : dlnGAME_dT = GAME% d1val1 / GAME% val
985 :
986 : FSCR% d1val1 = dFSCR_dlnGAME * dlnGAME_dT
987 : USCR% d1val1 = dUSCR_dlnGAME * dlnGAME_dT
988 : PSCR% d1val1 = dPSCR_dlnGAME * dlnGAME_dT
989 : CVSCR% d1val1 = dCVSCR_dlnGAME * dlnGAME_dT
990 : PDTSCR% d1val1 = dPDTSCR_dlnGAME * dlnGAME_dT
991 : PDRSCR% d1val1 = dPDRSCR_dlnGAME * dlnGAME_dT
992 :
993 : ! dRho (val2) derivatives
994 : dlnRs_dRho = RS% d1val2 / RS% val
995 : dlnGAME_dRho = GAME% d1val2 / GAME% val
996 :
997 : FSCR% d1val2 = dFSCR_dlnRS * dlnRS_dRho + dFSCR_dlnGAME * dlnGAME_dRho
998 : USCR% d1val2 = dUSCR_dlnRS * dlnRS_dRho + dUSCR_dlnGAME * dlnGAME_dRho
999 : PSCR% d1val2 = dPSCR_dlnRS * dlnRS_dRho + dPSCR_dlnGAME * dlnGAME_dRho
1000 : CVSCR% d1val2 = dCVSCR_dlnRS * dlnRS_dRho + dCVSCR_dlnGAME * dlnGAME_dRho
1001 : PDTSCR% d1val2 = dPDTSCR_dlnRS * dlnRS_dRho + dPDTSCR_dlnGAME * dlnGAME_dRho
1002 : PDRSCR% d1val2 = dPDRSCR_dlnRS * dlnRS_dRho + dPDRSCR_dlnGAME * dlnGAME_dRho
1003 :
1004 : end if
1005 :
1006 : contains
1007 :
1008 : logical function check1(v, xv, str)
1009 : type(auto_diff_real_2var_order1), intent(in) :: v
1010 : real(dp), intent(in) :: xv
1011 : character (len=*), intent(in) :: str
1012 : real(dp) :: val
1013 : real(dp), parameter :: atol = 1d-8, rtol = 1d-6
1014 : include 'formats'
1015 : val = v%val
1016 : check1 = .false.
1017 : if (is_bad(xv)) then
1018 : write(*,*) 'is_bad ' // trim(str), xv, val, int(Zion), RS, GAME
1019 : return
1020 : end if
1021 : if (abs(val - xv) > atol + rtol*max(abs(val),abs(xv))) then
1022 : write(*,*) 'rel mismatch ' // trim(str), &
1023 : (val - xv)/max(abs(val),abs(xv),1d-99), &
1024 : xv, val, int(Zion), RS%val, GAME%val
1025 : call mesa_error(__FILE__,__LINE__,'FSCRliq8')
1026 : return
1027 : end if
1028 : check1 = .true.
1029 0 : end function check1
1030 :
1031 : end subroutine FSCRliq8
1032 :
1033 : ! ============== SUBROUTINES FOR THE SOLID STATE ================= *
1034 0 : subroutine FSCRsol8(RS,GAMI,Zion,TPT, &
1035 : FSCR,USCR,PSCR,S_SCR,CVSCR,PDTSCR,PDRSCR,ierr)
1036 : ! Version 28.05.08
1037 : ! undefined zero variable Q1DXG is wiped out 21.06.10
1038 : ! cosmetic change 16.05.13
1039 : ! Fit to the el.-ion screening in bcc or fcc Coulomb solid
1040 : ! Stems from FSCRsol8 v.09.06.07. Included a check for RS=0.
1041 : ! INPUT: RS - el. density parameter, GAMI - ion coupling parameter,
1042 : ! Zion - ion charge, TPT=T_p/T - ion quantum parameter
1043 : ! OUTPUT: FSCR - screening (e-i) free energy per kT per 1 ion,
1044 : ! USCR - internal energy per kT per 1 ion (screen.contrib.)
1045 : ! PSCR - pressure divided by (n_i kT) (screen.contrib.)
1046 : ! S_SCR - screening entropy contribution / (N_i k)
1047 : ! CVSCR - heat capacity per 1 ion (screen.contrib.)
1048 : ! PDTSCR,PDRSCR = PSCR + d PSCR / d ln(T,\rho)
1049 : type(auto_diff_real_2var_order1), intent(in) :: RS,GAMI
1050 : real(dp), intent(in) :: Zion
1051 : type(auto_diff_real_2var_order1) :: TPT
1052 : type(auto_diff_real_2var_order1), intent(out) :: FSCR,USCR,PSCR,S_SCR,CVSCR,PDTSCR,PDRSCR
1053 : integer, intent(out) :: ierr
1054 :
1055 : type(auto_diff_real_2var_order1) :: XSR, P1, P2, Finf, FinfX, FinfDX, FinfDXX
1056 0 : real(dp) :: Z13, ZLN, R1, R2, R3
1057 : type(auto_diff_real_2var_order1) :: Q1U, Q1D, Q1, Q1X, Q1XDX, Q1DX, Q1DXX
1058 : type(auto_diff_real_2var_order1) :: Y0, Y0DX, Y0DG, Y0DXX, Y0DGG, Y0DXG
1059 : type(auto_diff_real_2var_order1) :: Y1, Y1DX, Y1DG, Y1DXX, Y1DGG, Y1DXG
1060 : type(auto_diff_real_2var_order1) :: SA, SUPA, SUPADX, SUPADG, SUPADXX, SUPADGG, SUPADXG
1061 : type(auto_diff_real_2var_order1) :: EM2, EM2Y1
1062 : type(auto_diff_real_2var_order1) :: SB, SUPB, SUPBDX, SUPBDG, SUPBDXX, SUPBDGG, SUPBDXG
1063 : type(auto_diff_real_2var_order1) :: SUP, SUPX, SUPDX, SUPG, SUPDG, SUPDXX, SUPDGG, SUPDXG
1064 : type(auto_diff_real_2var_order1) :: GR3, GR3X, GR3DX, GR3G, GR3DG, GR3DXX, GR3DGG, GR3DXG
1065 : type(auto_diff_real_2var_order1) :: W, WDX, WDG, WDXX, WDGG, WDXG
1066 : type(auto_diff_real_2var_order1) :: FDX, FDG, FDXX, FDGG, FDXG
1067 :
1068 0 : real(dp) :: AP(4)
1069 : real(dp) :: ENAT
1070 : real(dp) :: TINY
1071 : real(dp) :: PX
1072 :
1073 :
1074 0 : AP(1) = 1.1866d0
1075 0 : AP(2) = 0.684d0
1076 0 : AP(3) = 17.9d0
1077 0 : AP(4) = 41.5d0
1078 :
1079 0 : ENAT=2.7182818285d0
1080 0 : TINY=1.d-19
1081 0 : PX=0.205d0
1082 :
1083 : ! data AP/1.1857,.663,17.1,40./,PX/.212/ ! for fcc lattice
1084 0 : ierr = 0
1085 0 : if (RS<0d0) then
1086 0 : ierr = -1
1087 0 : return
1088 : !call mesa_error(__FILE__,__LINE__,'FSCRsol8: RS < 0')
1089 : end if
1090 0 : if (RS<TINY) then
1091 0 : FSCR=0.d0
1092 0 : USCR=0.d0
1093 0 : PSCR=0.d0
1094 0 : S_SCR=0.d0
1095 0 : CVSCR=0.d0
1096 0 : PDTSCR=0.d0
1097 0 : PDRSCR=0.d0
1098 0 : return
1099 : end if
1100 0 : XSR=0.0140047d0/RS ! relativity parameter
1101 0 : Z13=pow(Zion,1d0/3d0)
1102 0 : P1=0.00352d0*(1d0-AP(1)/pow(Zion,0.267d0)+0.27d0/Zion)
1103 : P2=1d0+2.25d0/Z13* &
1104 0 : (1d0+AP(2)*pow5(Zion)+0.222d0*pow6(Zion))/(1d0+0.222d0*pow6(Zion))
1105 0 : ZLN=log(Zion)
1106 0 : Finf=sqrt(P2/(XSR*XSR)+1d0)*Z13*Z13*P1 ! The TF limit
1107 0 : FinfX=-P2/((P2+XSR*XSR)*XSR)
1108 0 : FinfDX=Finf*FinfX
1109 0 : FinfDXX=FinfDX*FinfX-FinfDX*(P2+3d0*XSR*XSR)/((P2+XSR*XSR)*XSR)
1110 0 : R1=AP(4)/(1d0+ZLN)
1111 0 : R2=0.395d0*ZLN+0.347d0/Zion/sqrt(Zion)
1112 0 : R3=1d0/(1d0+ZLN*sqrt(ZLN)*0.01d0+0.097d0/(Zion*Zion))
1113 0 : Q1U=R1+AP(3)*XSR*XSR
1114 0 : Q1D=1d0+R2*XSR*XSR
1115 0 : Q1=Q1U/Q1D
1116 0 : Q1X=2d0*XSR*(AP(3)/Q1U-R2/Q1D)
1117 0 : Q1XDX=Q1X/XSR+4d0*XSR*XSR*((R2/Q1D)*(R2/Q1D)-(AP(3)/Q1U)*(AP(3)/Q1U))
1118 0 : Q1DX=Q1*Q1X
1119 0 : Q1DXX=Q1DX*Q1X+Q1*Q1XDX
1120 : ! New quantum factor, in order to suppress CVSCR at TPT >> 1
1121 0 : if (TPT<6d0/PX) then
1122 0 : Y0=(PX*TPT)*(PX*TPT)
1123 0 : Y0DX=Y0/XSR
1124 0 : Y0DG=2d0*Y0/GAMI
1125 0 : Y0DXX=0d0
1126 0 : Y0DGG=Y0DG/GAMI
1127 0 : Y0DXG=Y0DG/XSR
1128 0 : Y1=exp(Y0)
1129 0 : Y1DX=Y1*Y0DX
1130 0 : Y1DG=Y1*Y0DG
1131 0 : Y1DXX=Y1*(Y0DX*Y0DX+Y0DXX)
1132 0 : Y1DGG=Y1*(Y0DG*Y0DG+Y0DGG)
1133 0 : Y1DXG=Y1*(Y0DX*Y0DG+Y0DXG)
1134 0 : SA=1.d0+Y1
1135 0 : SUPA=log(SA)
1136 0 : SUPADX=Y1DX/SA
1137 0 : SUPADG=Y1DG/SA
1138 0 : SUPADXX=(Y1DXX-Y1DX*Y1DX/SA)/SA
1139 0 : SUPADGG=(Y1DGG-Y1DG*Y1DG/SA)/SA
1140 0 : SUPADXG=(Y1DXG-Y1DX*Y1DG/SA)/SA
1141 0 : EM2=ENAT-2.d0
1142 0 : SB=ENAT-EM2/Y1
1143 0 : SUPB=log(SB)
1144 0 : EM2Y1=EM2/(Y1*Y1*SB)
1145 0 : SUPBDX=EM2Y1*Y1DX
1146 0 : SUPBDG=EM2Y1*Y1DG
1147 0 : SUPBDXX=EM2Y1*(Y1DXX-2d0*Y1DX*Y1DX/Y1-Y1DX*SUPBDX)
1148 0 : SUPBDGG=EM2Y1*(Y1DGG-2d0*Y1DG*Y1DG/Y1-Y1DG*SUPBDG)
1149 0 : SUPBDXG=EM2Y1*(Y1DXG-2d0*Y1DX*Y1DG/Y1-Y1DG*SUPBDX)
1150 0 : SUP=sqrt(SUPA/SUPB)
1151 0 : SUPX=0.5d0*(SUPADX/SUPA-SUPBDX/SUPB)
1152 0 : SUPDX=SUP*SUPX
1153 0 : SUPG=0.5d0*(SUPADG/SUPA-SUPBDG/SUPB)
1154 0 : SUPDG=SUP*SUPG
1155 : SUPDXX=SUPDX*SUPX &
1156 : + SUP*0.5d0*(SUPADXX/SUPA-(SUPADX/SUPA)*(SUPADX/SUPA) &
1157 0 : - SUPBDXX/SUPB+(SUPBDX/SUPB)*(SUPBDX/SUPB))
1158 : SUPDGG=SUPDG*SUPG &
1159 : + SUP*0.5d0*(SUPADGG/SUPA-(SUPADG/SUPA)*(SUPADG/SUPA) &
1160 0 : - SUPBDGG/SUPB+(SUPBDG/SUPB)*(SUPBDG/SUPB))
1161 : SUPDXG=SUPDX*SUPG &
1162 : + SUP*0.5d0*((SUPADXG-SUPADX*SUPADG/SUPA)/SUPA &
1163 0 : - (SUPBDXG-SUPBDX*SUPBDG/SUPB)/SUPB)
1164 : else
1165 0 : SUP=PX*TPT
1166 0 : SUPDX=0.5d0*PX*TPT/XSR
1167 0 : SUPDG=PX*TPT/GAMI
1168 0 : SUPDXX=-0.5d0*SUPDX/XSR
1169 0 : SUPDGG=0d0
1170 0 : SUPDXG=SUPDX/GAMI
1171 : end if
1172 0 : GR3=pow(GAMI/SUP,R3)
1173 0 : GR3X=-R3*SUPDX/SUP
1174 0 : GR3DX=GR3*GR3X
1175 0 : GR3DXX=GR3DX*GR3X-R3*GR3*(SUPDXX/SUP-(SUPDX/SUP)*(SUPDX/SUP))
1176 0 : GR3G=R3*(1d0/GAMI-SUPDG/SUP)
1177 0 : GR3DG=GR3*GR3G
1178 0 : GR3DGG=GR3DG*GR3G+GR3*R3*((SUPDG/SUP)*(SUPDG/SUP)-SUPDGG/SUP-1d0/(GAMI*GAMI))
1179 0 : GR3DXG=GR3DG*GR3X+GR3*R3*(SUPDX*SUPDG/(SUP*SUP)-SUPDXG/SUP)
1180 0 : W=1d0+Q1/GR3
1181 0 : WDX=Q1DX/GR3-Q1*GR3DX/(GR3*GR3*GR3)
1182 0 : WDG=-Q1*GR3DG/(GR3*GR3)
1183 0 : WDXX=Q1DXX/GR3-(2d0*Q1DX*GR3DX+Q1*(GR3DXX-2d0*GR3DX*GR3DX/GR3))/(GR3*GR3)
1184 0 : WDGG=Q1*(2d0*GR3DG*GR3DG/GR3-GR3DGG)/(GR3*GR3)
1185 0 : WDXG=-(Q1DX*GR3DG+Q1*(GR3DXG-2d0*GR3DX*GR3DG/GR3))/(GR3*GR3)
1186 0 : FSCR=-GAMI*Finf*W
1187 0 : FDX=-GAMI*(FinfDX*W+Finf*WDX)
1188 0 : FDXX=-GAMI*(FinfDXX*W+2d0*FinfDX*WDX+Finf*WDXX)
1189 0 : FDG=-Finf*W-GAMI*Finf*WDG
1190 0 : FDGG=-2d0*Finf*WDG-GAMI*Finf*WDGG
1191 0 : FDXG=-FinfDX*W-Finf*WDX-GAMI*(FinfDX*WDG+Finf*WDXG)
1192 0 : S_SCR=-GAMI*GAMI*Finf*WDG
1193 0 : USCR=S_SCR+FSCR
1194 0 : CVSCR=-GAMI*GAMI*FDGG
1195 0 : PSCR=(XSR*FDX+GAMI*FDG)/3d0
1196 0 : PDTSCR=GAMI*GAMI*(XSR*Finf*(FinfX*WDG+WDXG)-FDGG)/3d0
1197 : PDRSCR=(12d0*PSCR+XSR*XSR*FDXX+2d0*XSR*GAMI*FDXG &
1198 0 : + GAMI*GAMI*FDGG)/9d0
1199 0 : return
1200 : end subroutine FSCRsol8
1201 :
1202 0 : subroutine ANHARM8(GAMI,TPT,Fah,Uah,Pah,CVah,PDTah,PDRah)
1203 : ! ANHARMONIC free energy Version 27.07.07
1204 : ! cleaned 16.06.09
1205 : ! Stems from ANHARM8b. Difference: AC=0., B1=.12 (.1217 - over accuracy)
1206 : ! Input: GAMI - ionic Gamma, TPT=Tp/T - ionic quantum parameter
1207 : ! Output: anharm.free en. Fah=F_{AH}/(N_i kT), internal energy Uah,
1208 : ! pressure Pah=P_{AH}/(n_i kT), specific heat CVah = C_{V,AH}/(N_i k),
1209 : ! PDTah = Pah + d Pah / d ln T, PDRah = Pah + d Pah / d ln\rho
1210 : type(auto_diff_real_2var_order1), intent(in) :: GAMI,TPT
1211 : type(auto_diff_real_2var_order1), intent(out) :: Fah,Uah,Pah,CVah,PDTah,PDRah
1212 : integer, parameter :: NM=3
1213 :
1214 : integer :: N
1215 0 : real(dp) :: AA(3)
1216 : type(auto_diff_real_2var_order1) :: CK, TPT2, TPT4, TQ, TK2, SUP
1217 : type(auto_diff_real_2var_order1) :: CN, SUPGN, ACN, PN
1218 : type(auto_diff_real_2var_order1) :: B1
1219 :
1220 0 : B1 = 0.12d0 ! coeff.at \eta^2/\Gamma at T=0
1221 :
1222 : ! Farouki & Hamaguchi'93
1223 0 : AA(1) = 10.9d0
1224 0 : AA(2) = 247.0d0
1225 0 : AA(3) = 1.765d5
1226 :
1227 :
1228 0 : CK=B1/AA(1) ! fit coefficient
1229 0 : TPT2=TPT*TPT
1230 0 : TPT4=TPT2*TPT2
1231 0 : TQ=B1*TPT2/GAMI ! quantum dependence
1232 0 : TK2=CK*TPT2
1233 0 : SUP=exp(-TK2) ! suppress.factor of class.anharmonicity
1234 0 : Fah=0.d0
1235 0 : Uah=0.d0
1236 0 : Pah=0.d0
1237 0 : CVah=0.d0
1238 0 : PDTah=0.d0
1239 0 : PDRah=0.d0
1240 0 : SUPGN=SUP
1241 0 : do N=1,NM
1242 0 : CN=1d0*N
1243 0 : SUPGN=SUPGN/GAMI ! SUP/Gamma^n
1244 0 : ACN=AA(N)
1245 0 : Fah=Fah-ACN/CN*SUPGN
1246 0 : Uah=Uah+(ACN*(1d0+2d0*TK2/CN))*SUPGN
1247 0 : PN=AA(N)/3d0+TK2*AA(N)/CN
1248 0 : Pah=Pah+PN*SUPGN
1249 : CVah=CVah+((CN+1d0)*AA(N)+(4d0-2d0/CN)*AA(N)*TK2 &
1250 0 : + 4d0*AA(N)*CK*CK/CN*TPT4)*SUPGN
1251 0 : PDTah=PDTah+(PN*(1d0+CN+2d0*TK2)-2d0/CN*AA(N)*TK2)*SUPGN
1252 0 : PDRah=PDRah+(PN*(1d0-CN/3d0-TK2)+AA(N)/CN*TK2)*SUPGN
1253 : end do
1254 0 : Fah=Fah-TQ
1255 0 : Uah=Uah-TQ
1256 0 : Pah=Pah-TQ/1.5d0
1257 0 : PDRah=PDRah-TQ/4.5d0
1258 0 : return
1259 : end subroutine ANHARM8
1260 :
1261 0 : subroutine FHARM12(GAMI,TPT, &
1262 : Fharm,Uharm,Pharm,CVth,Sharm,PDTharm,PDRharm)
1263 : ! Thermodynamic functions of a harmonic crystal, incl.stat.Coul.lattice
1264 : !
1265 : ! Version 27.04.12
1266 : ! Stems from FHARM8 v.15.02.08
1267 : ! Replaced HLfit8 with HLfit12: rearranged output.
1268 : ! Input: GAMI - ionic Gamma, TPT=T_{p,i}/T
1269 : ! Output: Fharm=F/(N_i T), Uharm=U/(N_i T), Pharm=P/(n_i T),
1270 : ! CVth=C_V/N_i, Sharm=S/N_i
1271 : ! PDTharm = Pharm + d Pharm / d ln T, PDRharm = Pharm + d Pharm/d ln\rho
1272 : type(auto_diff_real_2var_order1), intent(in) :: GAMI,TPT
1273 : type(auto_diff_real_2var_order1), intent(out) :: Fharm,Uharm,Pharm,CVth,Sharm,PDTharm,PDRharm
1274 :
1275 : type(auto_diff_real_2var_order1) :: Fth,Uth,Sth,U0,E0
1276 : type(auto_diff_real_2var_order1) :: F,U,U1
1277 :
1278 : real(dp), parameter :: CM = .895929256d0 ! Madelung
1279 :
1280 0 : call HLfit12(TPT,F,U,CVth,Sth,U1,1)
1281 0 : U0=-CM*GAMI ! perfect lattice
1282 0 : E0=1.5d0*U1*TPT ! zero-point energy
1283 0 : Uth=U+E0
1284 0 : Fth=F+E0
1285 0 : Uharm=U0+Uth
1286 0 : Fharm=U0+Fth
1287 0 : Sharm=Sth
1288 0 : Pharm=U0/3d0+Uth/2d0
1289 0 : PDTharm=0.5d0*CVth
1290 0 : PDRharm=U0/2.25d0+0.75d0*Uth-0.25d0*CVth
1291 0 : return
1292 : end subroutine FHARM12
1293 :
1294 0 : subroutine HLfit12(eta,F,U,CV,S,U1,LATTICE)
1295 : ! Version 24.04.12
1296 : ! Stems from HLfit8 v.03.12.08;
1297 : ! differences: E0 excluded from U and F;
1298 : ! U1 and d(CV)/d\ln(T) are added on the output.
1299 : ! Fit to thermal part of the thermodynamic functions.
1300 : ! Baiko, Potekhin, & Yakovlev (2001).
1301 : ! Zero-point lattice quantum energy 1.5u_1\eta EXCLUDED (unlike HLfit8).
1302 : ! Input: eta=Tp/T, LATTICE=1 for bcc, 2 for fcc
1303 : ! Output: F and U (normalized to NkT) - due to phonon excitations,
1304 : ! CV and S (normalized to Nk) in the HL model,
1305 : ! U1 - the 1st phonon moment,
1306 :
1307 : type(auto_diff_real_2var_order1) :: eta ! can be modified, not sure if this is an intended side-effect
1308 : type(auto_diff_real_2var_order1), intent(out) :: F, U, CV, S, U1
1309 : integer, intent(in) :: LATTICE
1310 :
1311 0 : real(dp) :: CLM, ALPHA, BETA, GAMMA
1312 0 : real(dp) :: A1, A2, A3, A4, A6, A8
1313 0 : real(dp) :: B0, B2, B4, B5, B6, B7, B9, B11
1314 0 : real(dp) :: C9, C11
1315 : type(auto_diff_real_2var_order1) :: UP, DN, EA, EB, EG, UP1, UP2, DN1, DN2, E0
1316 :
1317 : real(dp) :: EPS
1318 : real(dp) :: TINY
1319 :
1320 0 : EPS=1.d-5
1321 0 : TINY=1.d-99
1322 :
1323 0 : if (LATTICE==1) then ! bcc lattice
1324 0 : CLM=-2.49389d0 ! 3*ln<\omega/\omega_p>
1325 0 : U1=0.5113875d0
1326 0 : ALPHA=0.265764d0
1327 0 : BETA=0.334547d0
1328 0 : GAMMA=0.932446d0
1329 0 : A1=0.1839d0
1330 0 : A2=0.593586d0
1331 0 : A3=0.0054814d0
1332 0 : A4=5.01813d-4
1333 0 : A6=3.9247d-7
1334 0 : A8=5.8356d-11
1335 0 : B0=261.66d0
1336 0 : B2=7.07997d0
1337 0 : B4=0.0409484d0
1338 0 : B5=0.000397355d0
1339 0 : B6=5.11148d-5
1340 0 : B7=2.19749d-6
1341 0 : C9=0.004757014d0
1342 0 : C11=0.0047770935d0
1343 0 : elseif (LATTICE==2) then ! fcc lattice
1344 0 : CLM=-2.45373d0
1345 0 : U1=0.513194d0
1346 0 : ALPHA=0.257591d0
1347 0 : BETA=0.365284d0
1348 0 : GAMMA=0.9167070d0
1349 0 : A1=0.0d0
1350 0 : A2=0.532535d0
1351 0 : A3=0.0d0
1352 0 : A4=3.76545d-4
1353 0 : A6=2.63013d-7
1354 0 : A8=6.6318d-11
1355 0 : B0=303.20d0
1356 0 : B2=7.7255d0
1357 0 : B4=0.0439597d0
1358 0 : B5=0.000114295d0
1359 0 : B6=5.63434d-5
1360 0 : B7=1.36488d-6
1361 0 : C9=0.00492387d0
1362 0 : C11=0.00437506d0
1363 : else
1364 0 : call mesa_error(__FILE__,__LINE__,'HLfit: unknown lattice type')
1365 : end if
1366 0 : if (eta>1d0/EPS) then ! asymptote of Eq.(13) of BPY'01
1367 0 : U=3d0/(C11*eta*eta*eta)
1368 0 : F=-U/3d0
1369 0 : CV=4d0*U
1370 0 : else if (eta<EPS) then ! Eq.(17) of BPY'01
1371 0 : if (eta<TINY) eta = TINY !call mesa_error(__FILE__,__LINE__,'HLfit8: eta is too small')
1372 0 : F=3d0*log(eta)+CLM-1.5d0*U1*eta+eta*eta/24.d0
1373 0 : U=3d0-1.5d0*U1*eta+eta*eta/12d0
1374 0 : CV=3d0-eta*eta/12d0
1375 : else
1376 0 : B9=A6*C9
1377 0 : B11=A8*C11
1378 :
1379 : !UP=1d0+A1*eta+A2*eta**2+A3*eta**3+A4*eta**4+A6*eta**6+A8*eta**8
1380 0 : UP=1d0+eta*(A1+eta*(A2+eta*(A3+eta*(A4+eta*eta*(A6+eta*eta*A8)))))
1381 :
1382 : !DN=B0+B2*eta**2+B4*eta**4+B5*eta**5+B6*eta**6+B7*eta**7+B9*eta**9+B11*eta**11
1383 0 : DN=B0+eta*eta*(B2+eta*eta*(B4+eta*(B5+eta*(B6+eta*(B7+eta*eta*(B9+eta*eta*B11))))))
1384 :
1385 0 : EA=exp(-ALPHA*eta)
1386 0 : EB=exp(-BETA*eta)
1387 0 : EG=exp(-GAMMA*eta)
1388 0 : F=log(1.d0-EA)+log(1.d0-EB)+log(1.d0-EG)-UP/DN ! F_{thermal}/NT
1389 :
1390 : !UP1=A1+2d0*A2*eta+3.*A3*eta**2+4.*A4*eta**3+6d0*A6*eta**5+8.*A8*eta**7
1391 0 : UP1=A1+eta*(2d0*A2+eta*(3.d0*A3+eta*(4.d0*A4+eta*eta*(6d0*A6+eta*eta*8d0*A8))))
1392 :
1393 : !UP2=2d0*A2+6d0A3*eta+12d0*A4*eta**2+30.*A6*eta**4+56d0A8*eta**6
1394 0 : UP2=2d0*A2+eta*(6d0*A3+eta*(12d0*A4+eta*eta*(30d0*A6+eta*eta*56d0*A8)))
1395 :
1396 : !DN1=2d0*B2*eta+4.*B4*eta**3+5.*B5*eta**4+6d0*B6*eta**5+7.*B7*eta**6+9.*B9*eta**8+11d0*B11*eta**10.
1397 0 : DN1=eta*(2d0*B2+eta*eta*(4d0*B4+eta*(5d0*B5+eta*(6d0*B6+eta*(7d0*B7+eta*eta*(9d0*B9+eta*eta*11d0*B11))))))
1398 :
1399 : !DN2=2d0*B2+12d0*B4*eta**2+20.*B5*eta**3+30.*B6*eta**4+42d0*B7*eta**5+72d0*B9*eta**7+110.*B11*eta**9
1400 0 : DN2=2d0*B2+eta*eta*(12d0*B4+eta*(20d0*B5+eta*(30d0*B6+eta*(42d0*B7+eta*eta*(72d0*B9+eta*eta*110d0*B11)))))
1401 :
1402 : U=ALPHA*EA/(1.d0-EA)+BETA*EB/(1.d0-EB)+GAMMA*EG/(1.d0-EG) &
1403 0 : - (UP1*DN-DN1*UP)/(DN*DN) ! int.en./NT/eta
1404 : CV=ALPHA*ALPHA*EA/((1.d0-EA)*(1.d0-EA))+BETA*BETA*EB/((1.d0-EB)*(1.d0-EB)) &
1405 : + GAMMA*GAMMA*EG/((1.d0-EG)*(1.d0-EG)) &
1406 0 : + ((UP2*DN-DN2*UP)*DN-2d0*(UP1*DN-DN1*UP)*DN1)/(DN*DN*DN) ! cV/eta^2
1407 0 : U=U*eta
1408 0 : CV=CV*eta*eta
1409 : end if
1410 0 : S=U-F
1411 0 : return
1412 : end subroutine HLfit12
1413 :
1414 0 : subroutine CORMIX(RS,GAME,Zmean,Z2mean,Z52,Z53,Z321, &
1415 : FMIX,UMIX,PMIX,CVMIX,PDTMIX,PDRMIX)
1416 : ! Version 02.07.09
1417 : ! Correction to the linear mixing rule for moderate to small Gamma
1418 : ! Input: RS=r_s (if RS=0, then OCP, otherwise EIP)
1419 : ! GAME=\Gamma_e
1420 : ! Zmean=<Z> (average Z of all ions, without electrons)
1421 : ! Z2mean=<Z^2>, Z52=<Z^2.5>, Z53=<Z^{5/3}>, Z321=<Z(Z+1)^1.5>
1422 : ! Output: FMIX=\Delta f - corr.to the reduced free energy f=F/N_{ion}kT
1423 : ! UMIX=\Delta u - corr.to the reduced internal energy u
1424 : ! PMIX=\Delta u - corr.to the reduced pressure P=P/n_{ion}kT
1425 : ! CVMIX=\Delta c - corr.to the reduced heat capacity c_V
1426 : ! PDTMIX=(1/n_{ion}kT)d\Delta P / d ln T
1427 : ! = \Delta p + d \Delta p / d ln T
1428 : ! PDRMIX=(1/n_{ion}kT)d\Delta P / d ln n_e
1429 : ! (composition is assumed fixed: Zmean,Z2mean,Z52,Z53=constant)
1430 : type(auto_diff_real_2var_order1), intent(in) :: RS,GAME
1431 : real(dp), intent(in) :: Zmean,Z2mean,Z52,Z53,Z321
1432 : type(auto_diff_real_2var_order1), intent(out) :: FMIX,UMIX,PMIX,CVMIX,PDTMIX,PDRMIX
1433 :
1434 : type(auto_diff_real_2var_order1) :: GAMImean, Dif0, DifR, DifFDH, D
1435 : type(auto_diff_real_2var_order1) :: P3, D0, GP, FMIX0, Q, R, GQ, G, GDG, UDG
1436 :
1437 : real(dp) :: TINY
1438 :
1439 0 : TINY=1.d-9
1440 :
1441 0 : GAMImean=GAME*Z53
1442 0 : if (RS<TINY) then ! OCP
1443 0 : Dif0=Z52-sqrt(Z2mean*Z2mean*Z2mean/Zmean)
1444 : else
1445 0 : Dif0=Z321-sqrt((Z2mean+Zmean)*(Z2mean+Zmean)*(Z2mean+Zmean)/Zmean)
1446 : end if
1447 0 : DifR=Dif0/Z52
1448 0 : DifFDH=Dif0*GAME*sqrt(GAME/3d0) ! F_DH - F_LM(DH)
1449 0 : D=Z2mean/(Zmean*Zmean)
1450 0 : if (abs(D-1.d0)<TINY) then ! no correction
1451 0 : FMIX=0d0
1452 0 : UMIX=0d0
1453 0 : PMIX=0d0
1454 0 : CVMIX=0d0
1455 0 : PDTMIX=0d0
1456 0 : PDRMIX=0d0
1457 0 : return
1458 : end if
1459 0 : P3=pow(D,-0.2d0)
1460 0 : D0=(2.6d0*DifR+14d0*DifR*DifR*DifR)/(1.d0-P3)
1461 0 : GP=D0*pow(GAMImean,P3)
1462 0 : FMIX0=DifFDH/(1d0+GP)
1463 0 : Q=D*D*0.0117d0
1464 0 : R=1.5d0/P3-1d0
1465 0 : GQ=Q*GP
1466 0 : FMIX=FMIX0/pow(1d0+GQ,R)
1467 0 : G=1.5d0-P3*GP/(1d0+GP)-R*P3*GQ/(1d0+GQ)
1468 0 : UMIX=FMIX*G
1469 0 : PMIX=UMIX/3.d0
1470 0 : GDG=-P3*P3*(GP/((1.d0+GP)*(1.d0+GP))+R*GQ/((1.d0+GQ)*(1.d0+GQ))) ! d G /d ln Gamma
1471 0 : UDG=UMIX*G+FMIX*GDG ! d u_mix /d ln Gamma
1472 0 : CVMIX=UMIX-UDG
1473 0 : PDTMIX=PMIX-UDG/3d0
1474 0 : PDRMIX=PMIX+UDG/9d0
1475 0 : return
1476 : end subroutine CORMIX
1477 :
1478 : ! =================== IDEAL ELECTRON GAS =========================== *
1479 0 : subroutine ELECT11(TEMP,CHI, &
1480 : DENS,FEid,PEid,UEid,SEid,CVE,CHITE,CHIRE, &
1481 : DlnDH,DlnDT,DlnDHH,DlnDTT,DlnDHT)
1482 : ! Version 17.11.11
1483 : ! ELECT9 v.04.03.09 + smooth match of two fits at chi >> 1 + add.outputs
1484 : ! Compared to ELECTRON v.06.07.00, this S/R is completely rewritten:
1485 : ! numerical differentiation is avoided now.
1486 : ! Compared to ELECT7 v.06.06.07,
1487 : ! - call BLIN7 is changed to call BLIN9,
1488 : ! - Sommerfeld expansion is used at chi >~ 28 i.o. 1.e4
1489 : ! - Sommerfeld expansion is corrected: introduced DeltaEF, D1 and D2.
1490 : ! Ideal electron-gas EOS.
1491 : ! Input: TEMP - T [a.u.], CHI=\mu/T
1492 : ! Output: DENS - electron number density n_e [a.u.],
1493 : ! FEid - free energy / N_e kT, UEid - internal energy / N_e kT,
1494 : ! PEid - pressure (P_e) / n_e kT, SEid - entropy / N_e k,
1495 : ! CVE - heat capacity / N_e k,
1496 : ! CHITE=(d ln P_e/d ln T)_V, CHIRE=(d ln P_e/d ln n_e)_T
1497 : ! DlnDH=(d ln n_e/d CHI)_T = (T/n_e) (d n_e/d\mu)_T
1498 : ! DlnDT=(d ln n_e/d ln T)_CHI
1499 : ! DlnDHH=(d^2 ln n_e/d CHI^2)_T
1500 : ! DlnDTT=(d^2 ln n_e/d (ln T)^2)_CHI
1501 : ! DlnDHT=d^2 ln n_e/d (ln T) d CHI
1502 : type(auto_diff_real_2var_order1), intent(in) :: TEMP,CHI
1503 : type(auto_diff_real_2var_order1), intent(out) :: DENS,FEid,PEid,UEid,SEid,CVE,CHITE,CHIRE,DlnDH,DlnDT,DlnDHH,DlnDTT,DlnDHT
1504 :
1505 : type(auto_diff_real_2var_order1) :: X2, FP, FM
1506 : type(auto_diff_real_2var_order1) :: DENSa,FEida,PEida,UEida,SEida,CVEa,CHITEa,CHIREa,DlnDHa,DlnDTa,DlnDHHa,DlnDTTa,DlnDHTa
1507 : type(auto_diff_real_2var_order1) :: DENSb,FEidb,PEidb,UEidb,SEidb,CVEb,CHITEb,CHIREb,DlnDHb,DlnDTb,DlnDHHb,DlnDTTb,DlnDHTb
1508 :
1509 : type(auto_diff_real_2var_order1) :: CHI1
1510 : type(auto_diff_real_2var_order1) :: CHI2
1511 : type(auto_diff_real_2var_order1) :: XMAX
1512 : type(auto_diff_real_2var_order1) :: DCHI1
1513 : type(auto_diff_real_2var_order1) :: DCHI2
1514 : type(auto_diff_real_2var_order1) :: XSCAL2
1515 :
1516 0 : CHI1=0.6d0
1517 0 : CHI2=28.d0
1518 0 : XMAX=20.d0
1519 0 : DCHI1=0.1d0
1520 0 : DCHI2=CHI2-CHI1-DCHI1
1521 0 : XSCAL2=XMAX/DCHI2
1522 :
1523 0 : X2=(CHI-CHI2)*XSCAL2
1524 0 : if (X2<-XMAX) then
1525 : call ELECT11a(TEMP,CHI, &
1526 : DENS,FEid,PEid,UEid,SEid,CVE,CHITE,CHIRE,&
1527 0 : DlnDH,DlnDT,DlnDHH,DlnDTT,DlnDHT)
1528 0 : elseif (X2>XMAX) then
1529 : call ELECT11b(TEMP,CHI, &
1530 : DENS,FEid,PEid,UEid,SEid,CVE,CHITE,CHIRE,&
1531 0 : DlnDH,DlnDT,DlnDHH,DlnDTT,DlnDHT)
1532 : else
1533 0 : call FERMI10(X2,XMAX,FP,FM)
1534 : call ELECT11a(TEMP,CHI, &
1535 : DENSa,FEida,PEida,UEida,SEida,CVEa,CHITEa,CHIREa, &
1536 0 : DlnDHa,DlnDTa,DlnDHHa,DlnDTTa,DlnDHTa)
1537 : call ELECT11b(TEMP,CHI, &
1538 : DENSb,FEidb,PEidb,UEidb,SEidb,CVEb,CHITEb,CHIREb,&
1539 0 : DlnDHb,DlnDTb,DlnDHHb,DlnDTTb,DlnDHTb)
1540 0 : DENS=DENSa*FP+DENSb*FM
1541 0 : FEid=FEida*FP+FEidb*FM
1542 0 : PEid=PEida*FP+PEidb*FM
1543 0 : UEid=UEida*FP+UEidb*FM
1544 0 : SEid=SEida*FP+SEidb*FM
1545 0 : CVE=CVEa*FP+CVEb*FM
1546 0 : CHITE=CHITEa*FP+CHITEb*FM
1547 0 : CHIRE=CHIREa*FP+CHIREb*FM
1548 0 : DlnDH=DlnDHa*FP+DlnDHb*FM
1549 0 : DlnDT=DlnDTa*FP+DlnDTb*FM
1550 0 : DlnDHH=DlnDHHa*FP+DlnDHHb*FM
1551 0 : DlnDHT=DlnDHTa*FP+DlnDHTb*FM
1552 0 : DlnDTT=DlnDTTa*FP+DlnDTTb*FM
1553 : end if
1554 0 : return
1555 : end subroutine ELECT11
1556 :
1557 0 : subroutine ELECT11a(TEMP,CHI, &
1558 : DENS,FEid,PEid,UEid,SEid,CVE,CHITE,CHIRE, &
1559 : DlnDH,DlnDT,DlnDHH,DlnDTT,DlnDHT)
1560 : ! Version 16.11.11
1561 : ! This is THE FIRST PART of ELECT9 v.04.03.09.
1562 : type(auto_diff_real_2var_order1), intent(in) :: TEMP,CHI
1563 : type(auto_diff_real_2var_order1), intent(out) :: DENS,FEid,PEid,UEid,SEid,CVE,CHITE,CHIRE,DlnDH,DlnDT,DlnDHH,DlnDTT,DlnDHT
1564 :
1565 : type(auto_diff_real_2var_order1) :: TEMR
1566 : type(auto_diff_real_2var_order1) :: W0,W0DX,W0DT,W0DXX,W0DTT,W0DXT
1567 : type(auto_diff_real_2var_order1) :: W1,W1DX,W1DT,W1DXX,W1DTT,W1DXT
1568 : type(auto_diff_real_2var_order1) :: W2,W2DX,W2DT,W2DXX,W2DTT,W2DXT
1569 : type(auto_diff_real_2var_order1) :: W0XXX,W0XTT,W0XXT
1570 : type(auto_diff_real_2var_order1) :: TPI, DENR, PR, U
1571 : type(auto_diff_real_2var_order1) :: dndT, dPdT, dUdT, dndH, dPdH, dUdH
1572 : type(auto_diff_real_2var_order1) :: dndHH, dndHT, dndTT
1573 :
1574 : type(auto_diff_real_2var_order1) :: BOHR
1575 : type(auto_diff_real_2var_order1) :: PI2
1576 : type(auto_diff_real_2var_order1) :: BOHR2
1577 : type(auto_diff_real_2var_order1) :: BOHR3
1578 :
1579 0 : BOHR=137.036d0
1580 0 : PI2=PI*PI
1581 0 : BOHR2=BOHR*BOHR
1582 0 : BOHR3=BOHR2*BOHR !cleaned 15/6
1583 :
1584 0 : TEMR=TEMP/BOHR2 ! T in rel.units (=T/mc^2)
1585 : call BLIN9(TEMR,CHI, &
1586 : W0,W0DX,W0DT,W0DXX,W0DTT,W0DXT, &
1587 : W1,W1DX,W1DT,W1DXX,W1DTT,W1DXT, &
1588 : W2,W2DX,W2DT,W2DXX,W2DTT,W2DXT, &
1589 0 : W0XXX,W0XTT,W0XXT)
1590 0 : TPI=TEMR*sqrt(2.d0*TEMR)/PI2 ! common pre-factor
1591 0 : DENR=TPI*(W1*TEMR+W0)
1592 0 : PR=TEMR*TPI/3.d0*(W2*TEMR+2d0*W1)
1593 0 : U=TEMR*TPI*(W2*TEMR+W1)
1594 : ! (these are density, pressure, and internal energy in the rel.units)
1595 0 : PEid=PR/(DENR*TEMR)
1596 0 : UEid=U/(DENR*TEMR)
1597 0 : FEid=CHI-PEid
1598 0 : DENS=DENR*BOHR3 ! converts from rel.units to a.u.
1599 0 : SEid=UEid-FEid
1600 : ! derivatives over T at constant chi:
1601 0 : dndT=TPI*(1.5d0*W0/TEMR+2.5d0*W1+W0DT+TEMR*W1DT) ! (d n_e/dT)_\chi
1602 0 : dPdT=TPI/3.d0*(5.d0*W1+2d0*TEMR*W1DT+3.5d0*TEMR*W2+TEMR*TEMR*W2DT) !dP/dT
1603 0 : dUdT=TPI*(2.5d0*W1+TEMR*W1DT+3.5d0*TEMR*W2+TEMR*TEMR*W2DT) !dU/dT_\chi
1604 : ! derivatives over chi at constant T:
1605 0 : dndH=TPI*(W0DX+TEMR*W1DX) ! (d n_e/d\chi)_T
1606 0 : dndHH=TPI*(W0DXX+TEMR*W1DXX) ! (d^2 n_e/d\chi)_T
1607 0 : dndTT=TPI*(0.75d0*W0/pow2(TEMR)+3d0*W0DT/TEMR+W0DTT+3.75d0*W1/TEMR+5d0*W1DT+TEMR*W1DTT)
1608 0 : dndHT=TPI*(1.5d0*W0DX/TEMR+W0DXT+2.5d0*W1DX+TEMR*W1DXT)
1609 0 : DlnDH=dndH/DENR ! (d ln n_e/d\chi)_T
1610 0 : DlnDT=dndT*TEMR/DENR ! (d ln n_e/d ln T)_\chi
1611 0 : DlnDHH=dndHH/DENR-pow2(DlnDH) ! (d^2 ln n_e/d\chi^2)_T
1612 0 : DlnDTT=pow2(TEMR)/DENR*dndTT+DlnDT-pow2(DlnDT) ! d^2 ln n_e/d ln T^2
1613 0 : DlnDHT=TEMR/DENR*(dndHT-dndT*DlnDH) ! d^2 ln n_e/d\chi d ln T
1614 0 : dPdH=TPI/3d0*TEMR*(2d0*W1DX+TEMR*W2DX) ! (d P_e/d\chi)_T
1615 0 : dUdH=TPI*TEMR*(W1DX+TEMR*W2DX) ! (d U_e/d\chi)_T
1616 0 : CVE=(dUdT-dUdH*dndT/dndH)/DENR
1617 0 : CHITE=TEMR/PR*(dPdT-dPdH*dndT/dndH)
1618 0 : CHIRE=DENR/PR*dPdH/dndH ! (dndH*TEMR*PEid) ! DENS/PRE*dPdH/dndH
1619 0 : return
1620 : end subroutine ELECT11a
1621 :
1622 0 : subroutine ELECT11b(TEMP,CHI, &
1623 : DENS,FEid,PEid,UEid,SEid,CVE,CHITE,CHIRE, &
1624 : DlnDH,DlnDT,DlnDHH,DlnDTT,DlnDHT)
1625 : ! Version 17.11.11
1626 : ! Stems from ELECT9b v.19.01.10, Diff. - additional output.
1627 : ! Sommerfeld expansion at very large CHI.
1628 : type(auto_diff_real_2var_order1), intent(in) :: TEMP,CHI
1629 : type(auto_diff_real_2var_order1), intent(out) :: DENS,FEid,PEid,UEid,SEid,CVE,CHITE,CHIRE,DlnDH,DlnDT,DlnDHH,DlnDTT,DlnDHT
1630 :
1631 : type(auto_diff_real_2var_order1) :: TEMR, EF, DeltaEF, G, PF, F, DF, P, DelP, S, U
1632 : type(auto_diff_real_2var_order1) :: DENR, DT, D1, D2
1633 : type(auto_diff_real_2var_order1) :: TPI, dndH, dndT, dndHH, dndHT, dndTT
1634 : type(auto_diff_real_2var_order1) :: W0, W0DX, W0DT, W0DXX, W0DTT, W0DXT, &
1635 : W1,W1DX,W1DT,W1DXX,W1DTT,W1DXT,W2,W2DX,W2DT,W2DXX,W2DTT,W2DXT,W0XXX,W0XTT,W0XXT
1636 :
1637 : type(auto_diff_real_2var_order1) :: BOHR
1638 : type(auto_diff_real_2var_order1) :: PI2
1639 : type(auto_diff_real_2var_order1) :: BOHR2
1640 : type(auto_diff_real_2var_order1) :: BOHR3
1641 :
1642 0 : BOHR=137.036d0
1643 0 : PI2=PI*PI
1644 0 : BOHR2=BOHR*BOHR
1645 0 : BOHR3=BOHR2*BOHR !cleaned 15/6
1646 :
1647 0 : TEMR=TEMP/BOHR2 ! T in rel.units (=T/mc^2)
1648 0 : EF=CHI*TEMR ! Fermi energy in mc^2 - zeroth approx. = CMU1
1649 : DeltaEF=PI2*TEMR*TEMR/6.d0*(1.d0+2.d0*EF*(2.d0+EF)) &
1650 0 : /(EF*(1.d0+EF)*(2.d0+EF)) ! corr. [page 125, equiv. Eq.(6) of PC'10]]
1651 0 : EF=EF+DeltaEF ! corrected Fermi energy (14.02.09)
1652 0 : G=1.d0+EF ! electron Lorentz-factor
1653 0 : if (EF>1.d-5) then ! relativistic expansion (Yak.&Shal.'89)
1654 0 : PF=sqrt(G*G-1.d0) ! Fermi momentum [rel.un.=mc]
1655 0 : F=(PF*(1d0+2.d0*PF*PF)*G-PF*PF*PF/0.375d0-log(PF+G))/8.d0/PI2 !F/V
1656 0 : DF=-TEMR*TEMR*PF*G/6.d0 ! thermal correction to F/V
1657 0 : P=(PF*G*(PF*PF/1.5d0-1.d0)+log(PF+G))/8.d0/PI2 ! P(T=0)
1658 0 : DelP=TEMR*TEMR*PF*(PF*PF+2.d0)/G/18.d0 ! thermal correction to P
1659 0 : CVE=PI2*TEMR*G/(PF*PF)
1660 : else ! nonrelativistic limit
1661 0 : PF=sqrt(2.d0*EF)
1662 0 : F=pow5(PF)*0.1d0/PI2
1663 0 : DF=-TEMR*TEMR*PF/6.d0
1664 0 : P=F/1.5d0
1665 0 : DelP=TEMR*TEMR*PF/9.d0
1666 0 : CVE=PI2*TEMR/EF/2.d0
1667 : end if
1668 0 : F=F+DF
1669 0 : P=P+DelP
1670 0 : S=-2.d0*DF ! entropy per unit volume [rel.un.]
1671 0 : U=F+S
1672 0 : CHIRE=pow5(PF)/(9.d0*PI2*P*G)
1673 0 : CHITE=2.d0*DelP/P
1674 0 : DENR=PF*PF*PF/3.d0/PI2 ! n_e [rel.un.=\Compton^{-3}]
1675 0 : DENS=DENR*BOHR3 ! conversion to a.u.(=\Bohr_radius^{-3})
1676 : ! derivatives over chi at constant T and T at constant chi:
1677 0 : TPI=TEMR*sqrt(2.d0*TEMR)/PI2 ! common pre-factor
1678 : call SOMMERF(TEMR,CHI, &
1679 : W0,W0DX,W0DT,W0DXX,W0DTT,W0DXT, &
1680 : W1,W1DX,W1DT,W1DXX,W1DTT,W1DXT, &
1681 : W2,W2DX,W2DT,W2DXX,W2DTT,W2DXT,&
1682 0 : W0XXX,W0XTT,W0XXT)
1683 0 : dndH=TPI*(W0DX+TEMR*W1DX) ! (d n_e/d\chi)_T
1684 0 : dndT=TPI*(1.5d0*W0/TEMR+2.5d0*W1+W0DT+TEMR*W1DT) ! (d n_e/dT)_\chi
1685 0 : dndHH=TPI*(W0DXX+TEMR*W1DXX) ! (d^2 n_e/d\chi)_T
1686 0 : dndTT=TPI*(0.75d0*W0/pow2(TEMR)+3d0*W0DT/TEMR+W0DTT+3.75d0*W1/TEMR+5d0*W1DT+TEMR*W1DTT)
1687 0 : dndHT=TPI*(1.5d0*W0DX/TEMR+W0DXT+2.5d0*W1DX+TEMR*W1DXT)
1688 0 : DlnDH=dndH/DENR ! (d ln n_e/d\chi)_T
1689 0 : DlnDT=dndT*TEMR/DENR ! (d ln n_e/d ln T)_\chi
1690 0 : DlnDHH=dndHH/DENR-pow2(DlnDH) ! (d^2 ln n_e/d\chi^2)_T
1691 0 : DlnDTT=pow2(TEMR)/DENR*dndTT+DlnDT-pow2(DlnDT) ! d^2 ln n_e/d ln T^2
1692 0 : DlnDHT=TEMR/DENR*(dndHT-dndT*DlnDH) ! d^2 ln n_e/d\chi d ln T
1693 0 : DT=DENR*TEMR
1694 0 : PEid=P/DT
1695 0 : UEid=U/DT
1696 0 : FEid=F/DT
1697 0 : SEid=S/DT
1698 : ! Empirical corrections of 16.02.09:
1699 0 : D1=DeltaEF/EF
1700 0 : D2=D1*(4.d0-2.d0*(PF/G))
1701 0 : CVE=CVE/(1.d0+D2)
1702 0 : SEid=SEid/(1.d0+D1)
1703 0 : CHITE=CHITE/(1.d0+D2)
1704 0 : return
1705 0 : end subroutine ELECT11b
1706 :
1707 0 : subroutine SOMMERF(TEMR,CHI, &
1708 : W0,W0DX,W0DT,W0DXX,W0DTT,W0DXT, &
1709 : W1,W1DX,W1DT,W1DXX,W1DTT,W1DXT, &
1710 : W2,W2DX,W2DT,W2DXX,W2DTT,W2DXT, &
1711 : W0XXX,W0XTT,W0XXT)
1712 : ! Version 17.11.11
1713 : ! Sommerfeld expansion for the Fermi-Dirac integrals
1714 : ! Input: TEMR=T/mc^2; CHI=(\mu-mc^2)/T
1715 : ! Output: Wk - Fermi-Dirac integral of the order k+1/2
1716 : ! WkDX=dWk/dCHI, WkDT = dWk/dT, WkDXX=d^2 Wk / d CHI^2,
1717 : ! WkDTT=d^2 Wk / d T^2, WkDXT=d^2 Wk /dCHIdT,
1718 : ! W0XXX=d^3 W0 / d CHI^3, W0XTT=d^3 W0 /(d CHI d^2 T),
1719 : ! W0XXT=d^3 W0 /dCHI^2 dT
1720 : ! [Draft source: yellow book pages 124-127]
1721 :
1722 : type(auto_diff_real_2var_order1), intent(in) :: TEMR,CHI
1723 : type(auto_diff_real_2var_order1), intent(out) :: W0,W0DX,W0DT,W0DXX,W0DTT,W0DXT
1724 : type(auto_diff_real_2var_order1), intent(out) :: W1,W1DX,W1DT,W1DXX,W1DTT,W1DXT
1725 : type(auto_diff_real_2var_order1), intent(out) :: W2,W2DX,W2DT,W2DXX,W2DTT,W2DXT
1726 : type(auto_diff_real_2var_order1), intent(out) :: W0XXX,W0XTT,W0XXT
1727 :
1728 : type(auto_diff_real_2var_order1) :: CMU, CMU1, PIT26, CN0, CN1, CN2
1729 : type(auto_diff_real_2var_order1) :: CJ00,CJ10,CJ20,CJ01,CJ11, &
1730 : CJ21,CJ02,CJ12,CJ22,CJ03,CJ13,CJ23,CJ04,CJ14,CJ24,CJ05
1731 :
1732 0 : if (CHI<.5d0) call mesa_error(__FILE__,__LINE__,'SOMMERF: non-degenerate (small CHI)')
1733 0 : if (TEMR<=0.d0) call mesa_error(__FILE__,__LINE__,'SOMMERF: T < 0')
1734 0 : CMU1=CHI*TEMR ! chemical potential in rel.units
1735 0 : CMU=1.d0+CMU1
1736 : call SUBFERMJ(CMU1, &
1737 : CJ00,CJ10,CJ20, &
1738 : CJ01,CJ11,CJ21, &
1739 : CJ02,CJ12,CJ22, &
1740 : CJ03,CJ13,CJ23, &
1741 0 : CJ04,CJ14,CJ24,CJ05)
1742 0 : PIT26=pow2(PI*TEMR)/6.d0
1743 : !!! PITAU4=pow2(PIT26)*0.7d0
1744 0 : CN0=sqrt(.5d0/TEMR)/TEMR
1745 0 : CN1=CN0/TEMR
1746 0 : CN2=CN1/TEMR
1747 0 : W0=CN0*(CJ00+PIT26*CJ02) ! +CN0*PITAU4*CJ04
1748 0 : W1=CN1*(CJ10+PIT26*CJ12) ! +CN1*PITAU4*CJ14
1749 0 : W2=CN2*(CJ20+PIT26*CJ22) ! +CN2*PITAU4*CJ24
1750 0 : W0DX=CN0*TEMR*(CJ01+PIT26*CJ03) ! +CN0*PITAU4*CJ05
1751 0 : W1DX=CN0*(CJ11+PIT26*CJ13)
1752 0 : W2DX=CN1*(CJ21+PIT26*CJ23)
1753 0 : W0DT=CN1*(CMU1*CJ01-1.5d0*CJ00+PIT26*(CMU1*CJ03+.5d0*CJ02))
1754 : !!! + CN1*PITAU4*(CMU1*CJ05+2.5d0*CJ04)
1755 0 : W1DT=CN2*(CMU1*CJ11-2.5d0*CJ10+PIT26*(CMU1*CJ13-.5d0*CJ12))
1756 0 : W2DT=CN2/TEMR*(CMU1*CJ21-3.5d0*CJ20+PIT26*(CMU1*CJ23-1.5d0*CJ22))
1757 0 : W0DXX=CN0*pow2(TEMR)*(CJ02+PIT26*CJ04)
1758 0 : W1DXX=CN0*TEMR*(CJ12+PIT26*CJ14)
1759 0 : W2DXX=CN0*(CJ22+PIT26*CJ24)
1760 0 : W0DXT=CN0*(CMU1*CJ02-.5d0*CJ01+PIT26*(CMU1*CJ04+1.5d0*CJ03))
1761 0 : W1DXT=CN1*(CMU1*CJ12-1.5d0*CJ11+PIT26*(CMU1*CJ14+.5d0*CJ13))
1762 0 : W2DXT=CN2*(CMU1*CJ22-2.5d0*CJ21+PIT26*(CMU1*CJ24-.5d0*CJ23))
1763 0 : W0DTT=CN2*(3.75d0*CJ00-3.d0*CMU1*CJ01+pow2(CMU1)*CJ02+PIT26*(-.25d0*CJ02+CMU1*CJ03+pow2(CMU1)*CJ04))
1764 0 : W1DTT=CN2/TEMR*(8.75d0*CJ10-5.d0*CMU1*CJ11+pow2(CMU1)*CJ12+PIT26*(.75d0*CJ12-CMU1*CJ13+pow2(CMU1)*CJ14))
1765 0 : W2DTT=CN2/pow2(TEMR)*(15.75d0*CJ20-7.d0*CMU1*CJ21+pow2(CMU1)*CJ22+PIT26*(3.75d0*CJ22-3.d0*CMU1*CJ23+pow2(CMU1)*CJ24))
1766 0 : W0XXX=CN0*pow3(TEMR)*(CJ03+PIT26*CJ05)
1767 0 : W0XXT=CN0*TEMR*(CMU1*CJ03+.5d0*CJ02+PIT26*(CMU1*CJ05+2.5d0*CJ04))
1768 0 : W0XTT=CN1*(.75d0*CJ01-CMU1*CJ02+pow2(CMU1)*CJ03+PIT26*(.75d0*CJ03+3.d0*CMU1*CJ04+pow2(CMU1)*CJ05))
1769 0 : return
1770 : end subroutine SOMMERF
1771 :
1772 0 : subroutine SUBFERMJ(CMU1, &
1773 : CJ00,CJ10,CJ20, &
1774 : CJ01,CJ11,CJ21, &
1775 : CJ02,CJ12,CJ22, &
1776 : CJ03,CJ13,CJ23, &
1777 : CJ04,CJ14,CJ24,CJ05)
1778 : ! Version 17.11.11
1779 : ! Supplement to SOMMERF
1780 :
1781 : type(auto_diff_real_2var_order1), intent(in) :: CMU1
1782 : type(auto_diff_real_2var_order1), intent(out) :: CJ00,CJ10,CJ20,CJ01,CJ11,&
1783 : CJ21,CJ02,CJ12,CJ22,CJ03,CJ13,CJ23,CJ04,CJ14,CJ24,CJ05
1784 :
1785 : type(auto_diff_real_2var_order1) :: X0, X3, X5, CL, CMU
1786 :
1787 0 : if (CMU1<=0.d0) call mesa_error(__FILE__,__LINE__,'SUBFERMJ: small CHI')
1788 0 : CMU=1.d0+CMU1
1789 0 : X0=sqrt(CMU1*(2.d0+CMU1))
1790 0 : X3=pow3(X0)
1791 0 : X5=pow5(X0)
1792 0 : if (X0<1d-4) then
1793 0 : CJ00=X3/3.d0
1794 0 : CJ10=0.1d0*X5
1795 0 : CJ20=pow7(X0)/28.d0
1796 : else
1797 0 : CL=log(X0+CMU)
1798 0 : CJ00=0.5d0*(X0*CMU-CL) ! J_{1/2}^0
1799 0 : CJ10=X3/3d0-CJ00 ! J_{3/2}^0
1800 0 : CJ20=(0.75d0*CMU-2d0)/3d0*X3+1.25d0*CJ00 ! J_{5/2}^0
1801 : end if
1802 0 : CJ01=X0 ! J_{1/2}^1
1803 0 : CJ11=CJ01*CMU1 ! J_{3/2}^1
1804 0 : CJ21=CJ11*CMU1 ! J_{5/2}^1
1805 0 : CJ02=CMU/X0 ! J_{1/2}^2
1806 0 : CJ12=CMU1/X0*(3.d0+2.d0*CMU1) ! J_{3/2}^2
1807 0 : CJ22=pow2(CMU1)/X0*(5.d0+3.d0*CMU1) ! J_{5/2}^2
1808 0 : CJ03=-1.d0/X3 ! J_{1/2}^3
1809 0 : CJ13=CMU1/X3*(2.d0*pow2(CMU1)+6.d0*CMU1+3.d0)
1810 0 : CJ23=pow2(CMU1)/X3*(6.d0*pow2(CMU1)+2.d1*CMU1+1.5d1)
1811 0 : CJ04=3.d0*CMU/X5
1812 0 : CJ14=-3.d0*CMU1/X5
1813 0 : CJ24=pow2(CMU1)/X5*(6.d0*pow3(CMU1)+3.d1*pow2(CMU1)+45.d0*CMU1+15.d0)
1814 0 : CJ05=(-12.d0*pow2(CMU1)-24.d0*CMU1-15.d0)/(X5*pow2(X0))
1815 0 : return
1816 : end subroutine SUBFERMJ
1817 :
1818 0 : subroutine FERMI10(X,XMAX,FP,FM)
1819 : ! Version 20.01.10
1820 : ! Fermi distribution function and its 3 derivatives
1821 : ! Input: X - argument f(x)
1822 : ! XMAX - max|X| where it is assumed that 0 < f(x) < 1.
1823 : ! Output: FP = f(x)
1824 : ! FM = 1-f(x)
1825 : type(auto_diff_real_2var_order1), intent(in) :: X
1826 : type(auto_diff_real_2var_order1) :: XMAX ! not sure if this side-effect is desired
1827 : type(auto_diff_real_2var_order1), intent(out) :: FP, FM
1828 :
1829 0 : if (XMAX<3.d0) XMAX = 3d0 !call mesa_error(__FILE__,__LINE__,'FERMI: XMAX')
1830 0 : if (X>XMAX) then
1831 0 : FP=0.d0
1832 0 : FM=1.d0
1833 0 : elseif (X<-XMAX) then
1834 0 : FP=1.d0
1835 0 : FM=0.d0
1836 : else
1837 0 : FP=1.d0/(exp(X)+1.d0)
1838 0 : FM=1.d0-FP
1839 : end if
1840 0 : return
1841 : end subroutine FERMI10
1842 :
1843 : ! ============== ELECTRON EXCHANGE AND CORRELATION ================ *
1844 0 : subroutine EXCOR7(RS,GAME,FXC,UXC,PXC,CVXC,SXC,PDTXC,PDRXC)
1845 : ! Version 09.06.07
1846 : ! Exchange-correlation contribution for the electron gas
1847 : ! Stems from TANAKA1 v.03.03.96. Added derivatives.
1848 : ! Input: RS - electron density parameter =electron-sphere radius in a.u.
1849 : ! GAME - electron Coulomb coupling parameter
1850 : ! Output: FXC - excess free energy of e-liquid per kT per one electron
1851 : ! according to Tanaka & Ichimaru 85-87 and Ichimaru 93
1852 : ! UXC - internal energy contr.[per 1 electron, kT]
1853 : ! PXC - pressure contribution divided by (n_e kT)
1854 : ! CVXC - heat capacity divided by N_e k
1855 : ! SXC - entropy divided by N_e k
1856 : ! PDTXC,PDRXC = PXC + d ln PXC / d ln(T,\rho)
1857 : use pc_support
1858 : type(auto_diff_real_2var_order1), intent(in) :: RS, GAME
1859 : type(auto_diff_real_2var_order1), intent(out) :: FXC,UXC,PXC,CVXC,SXC,PDTXC,PDRXC
1860 :
1861 : type(auto_diff_real_2var_order1) :: THETA, SQTH, THETA2, THETA3, THETA4, EXP1TH
1862 : type(auto_diff_real_2var_order1) :: CHT1, SHT1, CHT2, SHT2
1863 : type(auto_diff_real_2var_order1) :: T1, T2, T1DH, T1DHH, T2DH, T2DHH
1864 : type(auto_diff_real_2var_order1) :: A0, A0DH, A0DHH, A1, A1DH, A1DHH, A, AH, ADH, ADHH
1865 : type(auto_diff_real_2var_order1) :: B0, B0DH, B0DHH, B1, B1DH, B1DHH, B, BH, BDH, BDHH
1866 : type(auto_diff_real_2var_order1) :: C, CDH, CDHH, C3, C3DH, C3DHH
1867 : type(auto_diff_real_2var_order1) :: D0, D0DH, D0DHH, D1, D1DH, D1DHH, D, DH, DDH, DDHH
1868 : type(auto_diff_real_2var_order1) :: E0, E0DH, E0DHH, E1, E1DH, E1DHH, E, EH, EDH, EDHH
1869 : type(auto_diff_real_2var_order1) :: DISCR, DIDH, DIDHH
1870 : type(auto_diff_real_2var_order1) :: S1, S1H, S1DH, S1DHH, S1DG, S1DHG
1871 : type(auto_diff_real_2var_order1) :: B2, B2DH, B2DHH, SQGE, B3, B3DH, B3DHH
1872 : type(auto_diff_real_2var_order1) :: S2, S2H, S2DH, S2DHH, S2DG, S2DGG, S2DHG
1873 : type(auto_diff_real_2var_order1) :: R3, R3H, R3DH, R3DHH, R3DG, R3DGG, R3DHG
1874 : type(auto_diff_real_2var_order1) :: S3, S3DH, S3DHH, S3DG, S3DGG, S3DHG
1875 : type(auto_diff_real_2var_order1) :: B4, B4DH, B4DHH
1876 : type(auto_diff_real_2var_order1) :: C4, C4DH, C4DHH, C4DG, C4DGG, C4DHG
1877 : type(auto_diff_real_2var_order1) :: S4A, S4AH, S4ADH, S4ADHH, S4ADG, S4ADGG, S4ADHG
1878 : type(auto_diff_real_2var_order1) :: S4B, S4BDH, S4BDHH, UP1, DN1, UP2, DN2
1879 : type(auto_diff_real_2var_order1) :: S4C, S4CDH, S4CDHH, S4CDG, S4CDGG, S4CDHG
1880 : type(auto_diff_real_2var_order1) :: S4, S4DH, S4DHH, S4DG, S4DGG, S4DHG
1881 : type(auto_diff_real_2var_order1) :: FXCDH, FXCDHH, FXCDG, FXCDGG, FXCDHG
1882 : type(auto_diff_real_2var_order1) :: PDLH, PDLG
1883 :
1884 : real(dp) :: &
1885 : xFXC, dFXC_dlnRS, dFXC_dlnGAME, &
1886 : xUXC, dUXC_dlnRS, dUXC_dlnGAME, &
1887 : xPXC, dPXC_dlnRS, dPXC_dlnGAME, &
1888 : xCVXC, dCVXC_dlnRS, dCVXC_dlnGAME, &
1889 : xSXC, dSXC_dlnRS, dSXC_dlnGAME, &
1890 : xPDTXC, dPDTXC_dlnRS, dPDTXC_dlnGAME, &
1891 : xPDRXC, dPDRXC_dlnRS, dPDRXC_dlnGAME
1892 0 : real(dp) :: dlnRs_dT, dlnRs_dRho, dlnGAME_dT, dlnGAME_dRho
1893 : integer :: ierr
1894 : logical :: skip
1895 : logical, parameter :: use_EXCOR7_table = .true.
1896 : logical, parameter :: debug_EXCOR7_table = .false.
1897 :
1898 : if (use_EXCOR7_table .or. debug_EXCOR7_table) then
1899 : ierr = 0
1900 : call get_EXCOR7(RS%val, GAME%val, &
1901 : xFXC, dFXC_dlnRS, dFXC_dlnGAME, &
1902 : xUXC, dUXC_dlnRS, dUXC_dlnGAME, &
1903 : xPXC, dPXC_dlnRS, dPXC_dlnGAME, &
1904 : xCVXC, dCVXC_dlnRS, dCVXC_dlnGAME, &
1905 : xSXC, dSXC_dlnRS, dSXC_dlnGAME, &
1906 : xPDTXC, dPDTXC_dlnRS, dPDTXC_dlnGAME, &
1907 0 : xPDRXC, dPDRXC_dlnRS, dPDRXC_dlnGAME, skip, ierr)
1908 0 : if (ierr /= 0) return
1909 : else
1910 : skip = .true.
1911 : end if
1912 :
1913 0 : if (.not. use_EXCOR7_table .or. debug_EXCOR7_table .or. skip) then
1914 :
1915 0 : THETA=0.543d0*RS/GAME ! non-relativistic degeneracy parameter
1916 0 : SQTH=sqrt(THETA)
1917 0 : THETA2=THETA*THETA
1918 0 : THETA3=THETA2*THETA
1919 0 : THETA4=THETA3*THETA
1920 0 : if (THETA>.04d0) then
1921 0 : CHT1=cosh(1.d0/THETA)
1922 0 : SHT1=sinh(1.d0/THETA)
1923 0 : CHT2=cosh(1.d0/SQTH)
1924 0 : SHT2=sinh(1.d0/SQTH)
1925 0 : T1=SHT1/CHT1 ! dtanh(1.d0/THETA)
1926 0 : T2=SHT2/CHT2 ! dtanh(1./sqrt(THETA))
1927 0 : T1DH=-1.d0/((THETA*CHT1)*(THETA*CHT1)) ! d T1 / d\theta
1928 0 : T1DHH=2.d0/pow3(THETA*CHT1)*(CHT1-SHT1/THETA)
1929 0 : T2DH=-0.5d0*SQTH/((THETA*CHT2)*(THETA*CHT2))
1930 0 : T2DHH=(0.75d0*SQTH*CHT2-0.5d0*SHT2)/pow3(THETA*CHT2)
1931 : else
1932 0 : T1=1.d0
1933 0 : T2=1.d0
1934 0 : T1DH=0.d0
1935 0 : T2DH=0.d0
1936 0 : T1DHH=0.d0
1937 0 : T2DHH=0.d0
1938 : end if
1939 0 : A0=0.75d0+3.04363d0*THETA2-0.09227d0*THETA3+1.7035d0*THETA4
1940 0 : A0DH=6.08726d0*THETA-0.27681d0*THETA2+6.814d0*THETA3
1941 0 : A0DHH=6.08726d0-0.55362d0*THETA+20.442d0*THETA2
1942 0 : A1=1d0+8.31051d0*THETA2+5.1105d0*THETA4
1943 0 : A1DH=16.62102d0*THETA+20.442d0*THETA3
1944 0 : A1DHH=16.62102d0+61.326d0*THETA2
1945 0 : A=0.610887d0*A0/A1*T1 ! HF fit of Perrot and Dharma-wardana
1946 0 : AH=A0DH/A0-A1DH/A1+T1DH/T1
1947 0 : ADH=A*AH
1948 : ADHH=ADH*AH+A*(A0DHH/A0-pow2(A0DH/A0)-A1DHH/A1+pow2(A1DH/A1) &
1949 0 : + T1DHH/T1-pow2(T1DH/T1))
1950 0 : B0=0.341308d0+12.070873d0*THETA2+1.148889d0*THETA4
1951 0 : B0DH=24.141746d0*THETA+4.595556d0*THETA3
1952 0 : B0DHH=24.141746d0+13.786668d0*THETA2
1953 0 : B1=1d0+10.495346d0*THETA2+1.326623d0*THETA4
1954 0 : B1DH=20.990692d0*THETA+5.306492d0*THETA3
1955 0 : B1DHH=20.990692d0+15.919476d0*THETA2
1956 0 : B=SQTH*T2*B0/B1
1957 0 : BH=0.5d0/THETA+T2DH/T2+B0DH/B0-B1DH/B1
1958 0 : BDH=B*BH
1959 : BDHH=BDH*BH+B*(-0.5d0/THETA2+T2DHH/T2-pow2(T2DH/T2) &
1960 0 : + B0DHH/B0-pow2(B0DH/B0)-B1DHH/B1+pow2(B1DH/B1))
1961 0 : D0=0.614925d0+16.996055d0*THETA2+1.489056d0*THETA4
1962 0 : D0DH=33.99211d0*THETA+5.956224d0*THETA3
1963 0 : D0DHH=33.99211d0+17.868672d0*THETA2
1964 0 : D1=1d0+10.10935d0*THETA2+1.22184d0*THETA4
1965 0 : D1DH=20.2187d0*THETA+4.88736d0*THETA3
1966 0 : D1DHH=20.2187d0+14.66208d0*THETA2
1967 0 : D=SQTH*T2*D0/D1
1968 0 : DH=0.5d0/THETA+T2DH/T2+D0DH/D0-D1DH/D1
1969 0 : DDH=D*DH
1970 : DDHH=DDH*DH+D*(-0.5d0/THETA2+T2DHH/T2-pow2(T2DH/T2) &
1971 0 : + D0DHH/D0-pow2(D0DH/D0)-D1DHH/D1+pow2(D1DH/D1))
1972 0 : E0=0.539409d0+2.522206d0*THETA2+0.178484d0*THETA4
1973 0 : E0DH=5.044412d0*THETA+0.713936d0*THETA3
1974 0 : E0DHH=5.044412d0+2.141808d0*THETA2
1975 0 : E1=1d0+2.555501d0*THETA2+0.146319d0*THETA4
1976 0 : E1DH=5.111002d0*THETA+0.585276d0*THETA3
1977 0 : E1DHH=5.111002d0+1.755828d0*THETA2
1978 0 : E=THETA*T1*E0/E1
1979 0 : EH=1.d0/THETA+T1DH/T1+E0DH/E0-E1DH/E1
1980 0 : EDH=E*EH
1981 : EDHH=EDH*EH+E*(T1DHH/T1-pow2(T1DH/T1)+E0DHH/E0-pow2(E0DH/E0) &
1982 0 : - E1DHH/E1+pow2(E1DH/E1)-1.0d0/THETA2)
1983 0 : EXP1TH=exp(-1.d0/THETA)
1984 0 : C=(0.872496d0+0.025248d0*EXP1TH)*E
1985 0 : CDH=0.025248d0*EXP1TH/THETA2*E+C*EDH/E
1986 : CDHH=0.025248d0*EXP1TH/THETA2*(EDH+(1.0d0-2.0d0*THETA)/THETA2*E) &
1987 0 : + CDH*EDH/E+C*EDHH/E-C*pow2(EDH/E)
1988 0 : DISCR=SQRT(4.0d0*E-D*D)
1989 0 : DIDH=0.5d0/DISCR*(4.0d0*EDH-2.0d0*D*DDH)
1990 0 : DIDHH=(-pow2((2.0d0*EDH-D*DDH)/DISCR)+2.0d0*EDHH-DDH*DDH-D*DDHH)/DISCR
1991 0 : S1=-C/E*GAME
1992 0 : S1H=CDH/C-EDH/E
1993 0 : S1DH=S1*S1H
1994 0 : S1DHH=S1DH*S1H+S1*(CDHH/C-pow2(CDH/C)-EDHH/E+pow2(EDH/E))
1995 0 : S1DG=-C/E ! => S1DGG=0
1996 0 : S1DHG=S1DG*(CDH/C-EDH/E)
1997 0 : B2=B-C*D/E
1998 0 : B2DH=BDH-(CDH*D+C*DDH)/E+C*D*EDH/(E*E)
1999 0 : B2DHH=BDHH-(CDHH*D+2d0*CDH*DDH+C*DDHH)/E+(2d0*(CDH*D+C*DDH-C*D*EDH/E)*EDH+C*D*EDHH)/(E*E)
2000 0 : SQGE=SQRT(GAME)
2001 0 : S2=-2.d0/E*B2*SQGE
2002 0 : S2H=B2DH/B2-EDH/E
2003 0 : S2DH=S2*S2H
2004 0 : S2DHH=S2DH*S2H+S2*(B2DHH/B2-pow2(B2DH/B2)-EDHH/E+pow2(EDH/E))
2005 0 : S2DG=0.5d0*S2/GAME
2006 0 : S2DGG=-0.5d0*S2DG/GAME
2007 0 : S2DHG=0.5d0*S2DH/GAME
2008 0 : R3=E*GAME+D*SQGE+1.0d0
2009 0 : R3DH=EDH*GAME+DDH*SQGE
2010 0 : R3DHH=EDHH*GAME+DDHH*SQGE
2011 0 : R3DG=E+0.5d0*D/SQGE
2012 0 : R3DGG=-0.25d0*D/(GAME*SQGE)
2013 0 : R3DHG=EDH+0.5d0*DDH/SQGE
2014 0 : B3=A-C/E
2015 0 : B3DH=ADH-CDH/E+C*EDH/(E*E)
2016 0 : B3DHH=ADHH-CDHH/E+(2.0d0*CDH*EDH+C*EDHH)/(E*E)-2d0*C*EDH*EDH/pow3(E)
2017 0 : C3=(D/E*B2-B3)/E
2018 0 : C3DH=(DDH*B2+D*B2DH+B3*EDH)/(E*E)-2d0*D*B2*EDH/pow3(E)-B3DH/E
2019 : C3DHH=(-B3DHH &
2020 : + (DDHH*B2+2d0*DDH*B2DH+D*B2DHH+B3DH*EDH+B3*EDHH+B3DH*EDH)/E &
2021 : - 2.0d0*((DDH*B2+D*B2DH+B3*EDH+DDH*B2+D*B2DH)*EDH+D*B2*EDHH)/(E*E) &
2022 0 : + 6.0d0*D*B2*EDH*EDH/pow3(E))/E
2023 0 : S3=C3*log(R3)
2024 0 : S3DH=S3*C3DH/C3+C3*R3DH/R3
2025 : S3DHH=(S3DH*C3DH+S3*C3DHH)/C3-S3*pow2(C3DH/C3) &
2026 0 : + (C3DH*R3DH+C3*R3DHH)/R3-C3*pow2(R3DH/R3)
2027 0 : S3DG=C3*R3DG/R3
2028 0 : S3DGG=C3*(R3DGG/R3-pow2(R3DG/R3))
2029 0 : S3DHG=(C3DH*R3DG+C3*R3DHG)/R3-C3*R3DG*R3DH/(R3*R3)
2030 0 : B4=2.d0-D*D/E
2031 0 : B4DH=EDH*(D/E)*(D/E)-2d0*D*DDH/E
2032 : B4DHH=EDHH*(D/E)*(D/E)+2d0*EDH*(D/E)*(D/E)*(DDH/D-EDH/E) &
2033 0 : - 2d0*(DDH*DDH+D*DDHH)/E+2d0*D*DDH*EDH/(E*E)
2034 0 : C4=2d0*E*SQGE+D
2035 0 : C4DH=2d0*EDH*SQGE+DDH
2036 0 : C4DHH=2d0*EDHH*SQGE+DDHH
2037 0 : C4DG=E/SQGE
2038 0 : C4DGG=-0.5d0*E/(GAME*SQGE)
2039 0 : C4DHG=EDH/SQGE
2040 0 : S4A=2.0d0/E/DISCR
2041 0 : S4AH=EDH/E+DIDH/DISCR
2042 0 : S4ADH=-S4A*S4AH
2043 0 : S4ADHH=-S4ADH*S4AH - S4A*(EDHH/E-(EDH/E)*(EDH/E)+DIDHH/DISCR-pow2(DIDH/DISCR))
2044 0 : S4B=D*B3+B4*B2
2045 0 : S4BDH=DDH*B3+D*B3DH+B4DH*B2+B4*B2DH
2046 0 : S4BDHH=DDHH*B3+2d0*DDH*B3DH+D*B3DHH+B4DHH*B2+2d0*B4DH*B2DH+B4*B2DHH
2047 0 : S4C=atan(C4/DISCR)-atan(D/DISCR)
2048 0 : UP1=C4DH*DISCR-C4*DIDH
2049 0 : DN1=DISCR*DISCR+C4*C4
2050 0 : UP2=DDH*DISCR-D*DIDH
2051 0 : DN2=DISCR*DISCR+D*D
2052 0 : S4CDH=UP1/DN1-UP2/DN2
2053 : S4CDHH=(C4DHH*DISCR-C4*DIDHH)/DN1 &
2054 : - UP1*2d0*(DISCR*DIDH+C4*C4DH)/(DN1*DN1) &
2055 0 : - (DDHH*DISCR-D*DIDHH)/DN2+UP2*2d0*(DISCR*DIDH+D*DDH)/(DN2*DN2)
2056 0 : S4CDG=C4DG*DISCR/DN1
2057 0 : S4CDGG=C4DGG*DISCR/DN1-2d0*C4*DISCR*pow2(C4DG/DN1)
2058 0 : S4CDHG=(C4DHG*DISCR+C4DG*DIDH-C4DG*DISCR/DN1*2d0*(DISCR*DIDH+C4*C4DH))/DN1
2059 0 : S4=S4A*S4B*S4C
2060 0 : S4DH=S4ADH*S4B*S4C+S4A*S4BDH*S4C+S4A*S4B*S4CDH
2061 : S4DHH=S4ADHH*S4B*S4C+S4A*S4BDHH*S4C+S4A*S4B*S4CDHH &
2062 0 : + 2d0*(S4ADH*S4BDH*S4C+S4ADH*S4B*S4CDH+S4A*S4BDH*S4CDH)
2063 0 : S4DG=S4A*S4B*S4CDG
2064 0 : S4DGG=S4A*S4B*S4CDGG
2065 0 : S4DHG=S4A*S4B*S4CDHG+S4CDG*(S4ADH*S4B+S4A*S4BDH)
2066 0 : FXC=S1+S2+S3+S4
2067 0 : FXCDH=S1DH+S2DH+S3DH+S4DH
2068 0 : FXCDG=S1DG+S2DG+S3DG+S4DG
2069 0 : FXCDHH=S1DHH+S2DHH+S3DHH+S4DHH
2070 0 : FXCDGG=S2DGG+S3DGG+S4DGG
2071 0 : FXCDHG=S1DHG+S2DHG+S3DHG+S4DHG
2072 0 : PXC=(GAME*FXCDG-2d0*THETA*FXCDH)/3.d0
2073 0 : UXC=GAME*FXCDG-THETA*FXCDH
2074 0 : SXC=(GAME*S2DG-S2+GAME*S3DG-S3+S4A*S4B*(GAME*S4CDG-S4C))-THETA*FXCDH
2075 0 : if (abs(SXC)<1.d-9*abs(THETA*FXCDH)) SXC=0.d0 ! accuracy loss
2076 0 : CVXC=2d0*THETA*(GAME*FXCDHG-FXCDH)-THETA*THETA*FXCDHH-GAME*GAME*FXCDGG
2077 0 : if (abs(CVXC)<1.d-9*abs(GAME*GAME*FXCDGG)) CVXC=0.d0 ! accuracy
2078 0 : PDLH=THETA*(GAME*FXCDHG-2d0*FXCDH-2d0*THETA*FXCDHH)/3.d0
2079 0 : PDLG=GAME*(FXCDG+GAME*FXCDGG-2d0*THETA*FXCDHG)/3.d0
2080 0 : PDRXC=PXC+(PDLG-2d0*PDLH)/3.d0
2081 0 : PDTXC=GAME*(THETA*FXCDHG-GAME*FXCDGG/3.d0)-THETA*(FXCDH/0.75d0+THETA*FXCDHH/1.5d0)
2082 :
2083 : end if
2084 :
2085 :
2086 : if (debug_EXCOR7_table .and. .not. skip) then
2087 :
2088 : if (.not. check1(FXC, xFXC, 'FXC')) return
2089 : if (.not. check1(UXC, xUXC, 'UXC')) return
2090 : if (.not. check1(PXC, xPXC, 'PXC')) return
2091 : if (.not. check1(CVXC, xCVXC, 'CVXC')) return
2092 : if (.not. check1(SXC, xSXC, 'SXC')) return
2093 : if (.not. check1(PDTXC, xPDTXC, 'PDTXC')) return
2094 : if (.not. check1(PDRXC, xPDRXC, 'PDRXC')) return
2095 :
2096 : end if
2097 :
2098 0 : if (use_EXCOR7_table .and. .not. skip) then
2099 :
2100 : ! values
2101 0 : FXC = xFXC
2102 0 : UXC = xUXC
2103 0 : PXC = xPXC
2104 0 : CVXC = xCVXC
2105 0 : SXC = xSXC
2106 0 : PDTXC = xPDTXC
2107 0 : PDRXC = xPDRXC
2108 :
2109 : ! dT (val1) derivatives
2110 : ! numerically, RS% d1val1 / RS% val is not quite 0
2111 0 : dlnRs_dT = RS% d1val1 / RS% val
2112 0 : dlnGAME_dT = GAME% d1val1 / GAME% val
2113 :
2114 0 : FXC% d1val1 = dFXC_dlnRS * dlnRS_dT + dFXC_dlnGAME * dlnGAME_dT
2115 0 : UXC% d1val1 = dUXC_dlnRS * dlnRS_dT + dUXC_dlnGAME * dlnGAME_dT
2116 0 : PXC% d1val1 = dPXC_dlnRS * dlnRS_dT + dPXC_dlnGAME * dlnGAME_dT
2117 0 : CVXC% d1val1 = dCVXC_dlnRS * dlnRS_dT + dCVXC_dlnGAME * dlnGAME_dT
2118 0 : SXC% d1val1 = dSXC_dlnRS * dlnRS_dT + dSXC_dlnGAME * dlnGAME_dT
2119 0 : PDTXC% d1val1 = dPDTXC_dlnRS * dlnRS_dT + dPDTXC_dlnGAME * dlnGAME_dT
2120 0 : PDRXC% d1val1 = dPDRXC_dlnRS * dlnRS_dT + dPDRXC_dlnGAME * dlnGAME_dT
2121 :
2122 : ! dRho (val2) derivatives
2123 0 : dlnRs_dRho = RS% d1val2 / RS% val
2124 0 : dlnGAME_dRho = GAME% d1val2 / GAME% val
2125 :
2126 0 : FXC% d1val2 = dFXC_dlnRS * dlnRS_dRho + dFXC_dlnGAME * dlnGAME_dRho
2127 0 : UXC% d1val2 = dUXC_dlnRS * dlnRS_dRho + dUXC_dlnGAME * dlnGAME_dRho
2128 0 : PXC% d1val2 = dPXC_dlnRS * dlnRS_dRho + dPXC_dlnGAME * dlnGAME_dRho
2129 0 : CVXC% d1val2 = dCVXC_dlnRS * dlnRS_dRho + dCVXC_dlnGAME * dlnGAME_dRho
2130 0 : SXC% d1val2 = dSXC_dlnRS * dlnRS_dRho + dSXC_dlnGAME * dlnGAME_dRho
2131 0 : PDTXC% d1val2 = dPDTXC_dlnRS * dlnRS_dRho + dPDTXC_dlnGAME * dlnGAME_dRho
2132 0 : PDRXC% d1val2 = dPDRXC_dlnRS * dlnRS_dRho + dPDRXC_dlnGAME * dlnGAME_dRho
2133 :
2134 : end if
2135 :
2136 : contains
2137 :
2138 : logical function check1(v, xv, str)
2139 : type(auto_diff_real_2var_order1), intent(in) :: v
2140 : real(dp), intent(in) :: xv
2141 : character (len=*), intent(in) :: str
2142 : real(dp) :: val
2143 : real(dp), parameter :: atol = 1d-8, rtol = 1d-6
2144 : include 'formats'
2145 : val = v%val
2146 : check1 = .false.
2147 : if (is_bad(xv)) then
2148 : write(*,*) 'is_bad ' // trim(str), xv, val, RS, GAME
2149 : return
2150 : end if
2151 : if (abs(val - xv) > atol + rtol*max(abs(val),abs(xv))) then
2152 : write(*,*) 'rel mismatch ' // trim(str), &
2153 : (val - xv)/max(abs(val),abs(xv),1d-99), &
2154 : xv, val, RS%val, GAME%val
2155 : call mesa_error(__FILE__,__LINE__,'EXCOR7')
2156 : return
2157 : end if
2158 : check1 = .true.
2159 0 : end function check1
2160 :
2161 :
2162 : end subroutine EXCOR7
2163 :
2164 : ! ====================== AUXILIARY SUBROUTINES ==================== *
2165 0 : subroutine FERINV7(F,N,X,XDF,XDFF) ! Inverse Fermi integrals
2166 : ! Version 24.05.07
2167 : ! X_q(f)=F^{-1}_q(f) : H.M.Antia 93 ApJS 84, 101
2168 : ! q=N-1/2=-1/2,1/2,3/2,5/2 (N=0,1,2,3)
2169 : ! Input: F - argument, N=q+1/2
2170 : ! Output: X=X_q, XDF=dX/df, XDFF=d^2 X / df^2
2171 : ! Relative error: N = 0 1 2 3
2172 : ! for X: 3.e-9, 4.2e-9, 2.3e-9, 6.2e-9
2173 : ! jump at f=4:
2174 : ! for XDF: 6.e-7, 5.4e-7, 9.6e-8, 3.1e-7
2175 : ! for XDFF: 4.7e-5, 4.8e-5, 2.3e-6, 1.5e-6
2176 : type(auto_diff_real_2var_order1) :: F ! can be modified, not sure if this is an intended side-effect
2177 : integer, intent(in) :: N
2178 : type(auto_diff_real_2var_order1), intent(out) :: X, XDF, XDFF
2179 : integer :: I
2180 : type(auto_diff_real_2var_order1) :: P, T, T1, T2, UP, UP1, UP2, DOWN, DOWN1, DOWN2
2181 : type(auto_diff_real_2var_order1) :: R, R1, R2, RT
2182 :
2183 : ! The next four are really parameters but there isn't a clean way to initialize them
2184 : ! at declaration time. - Adam Jermyn 4/2/2020
2185 0 : real(dp) :: A(0:5,0:3) ! read only after initialization
2186 0 : real(dp) :: B(0:6,0:3) ! read only after initialization
2187 0 : real(dp) :: C(0:6,0:3) ! read only after initialization
2188 0 : real(dp) :: D(0:6,0:3) ! read only after initialization
2189 :
2190 : integer, parameter :: LA(0:3) = [5,4,3,2]
2191 : integer, parameter :: LB(0:3) = [6,3,4,3]
2192 : integer, parameter :: LD(0:3) = [6,5,5,6]
2193 :
2194 : A(0:5,0) = [ &
2195 : -1.570044577033d4,1.001958278442d4,-2.805343454951d3, &
2196 0 : 4.121170498099d2,-3.174780572961d1,1.d0] ! X_{-1/2}
2197 : A(0:5,1) = [ &
2198 : 1.999266880833d4,5.702479099336d3,6.610132843877d2, &
2199 0 : 3.818838129486d1,1.d0,0d0] ! X_{1/2}
2200 : A(0:5,2) = [ &
2201 : 1.715627994191d2,1.125926232897d2,2.056296753055d1, &
2202 0 : 1.d0,0d0,0d0]
2203 : A(0:5,3) = [ &
2204 0 : 2.138969250409d2,3.539903493971d1,1.d0,0d0,0d0,0d0] ! X_{5/2}
2205 : B(0:6,0) = [ &
2206 : -2.782831558471d4,2.886114034012d4,-1.274243093149d4, &
2207 : 3.063252215963d3,-4.225615045074d2,3.168918168284d1, &
2208 0 : -1.008561571363d0] ! X_{-1/2}
2209 : B(0:6,1) = [ &
2210 : 1.771804140488d4,-2.014785161019d3,9.130355392717d1, &
2211 0 : -1.670718177489d0,0d0,0d0,0d0] ! X_{1/2}
2212 : B(0:6,2) = [ &
2213 : 2.280653583157d2,1.193456203021d2,1.16774311354d1, &
2214 0 : -3.226808804038d-1,3.519268762788d-3,0d0,0d0] ! X_{3/2}
2215 : B(0:6,3) = [ &
2216 : 7.10854551271d2,9.873746988121d1,1.067755522895d0, &
2217 0 : -1.182798726503d-2,0d0,0d0,0d0] ! X_{5/2}
2218 : C(0:6,0) = [ &
2219 : 2.206779160034d-8,-1.437701234283d-6,6.103116850636d-5, &
2220 : -1.169411057416d-3,1.814141021608d-2,-9.588603457639d-2, &
2221 0 : 1.d0]
2222 : C(0:6,1) = [ &
2223 : -1.277060388085d-2,7.187946804945d-2,-4.262314235106d-1, &
2224 : 4.997559426872d-1,-1.285579118012d0,-3.930805454272d-1, &
2225 0 : 1.d0]
2226 : C(0:6,2) = [ &
2227 : -6.321828169799d-3,-2.183147266896d-2,-1.05756279932d-1, &
2228 : -4.657944387545d-1,-5.951932864088d-1,3.6844711771d-1, &
2229 0 : 1.d0]
2230 : C(0:6,3) = [ &
2231 : -3.312041011227d-2,1.315763372315d-1,-4.820942898296d-1, &
2232 : 5.099038074944d-1,5.49561349863d-1,-1.498867562255d0, &
2233 0 : 1.d0]
2234 : D(0:6,0) = [ &
2235 : 8.827116613576d-8,-5.750804196059d-6,2.429627688357d-4, &
2236 : -4.601959491394d-3,6.932122275919d-2,-3.217372489776d-1, &
2237 0 : 3.124344749296d0] ! X_{-1/2}
2238 : D(0:6,1) = [ &
2239 : -9.745794806288d-3,5.485432756838d-2,-3.29946624326d-1, &
2240 : 4.077841975923d-1,-1.145531476975d0,-6.067091689181d-2, &
2241 0 : 0d0]
2242 : D(0:6,2) = [ &
2243 : -4.381942605018d-3,-1.5132365041d-2,-7.850001283886d-2, &
2244 : -3.407561772612d-1,-5.074812565486d-1,-1.387107009074d-1, &
2245 0 : 0d0]
2246 : D(0:6,3) = [ &
2247 : -2.315515517515d-2,9.198776585252d-2,-3.835879295548d-1, &
2248 : 5.415026856351d-1,-3.847241692193d-1,3.739781456585d-2, &
2249 0 : -3.008504449098d-2] ! X_{5/2}
2250 :
2251 0 : if (N<0d0 .or.N>3d0) call mesa_error(__FILE__,__LINE__,'FERINV7: Invalid subscript')
2252 0 : if (F<=0.d0) F = 1d-99 !call mesa_error(__FILE__,__LINE__,'FERINV7: Non-positive argument')
2253 0 : if (F<4.d0) then
2254 0 : T=F
2255 0 : UP=0.d0
2256 0 : UP1=0.d0
2257 0 : UP2=0.d0
2258 0 : DOWN=0.d0
2259 0 : DOWN1=0.d0
2260 0 : DOWN2=0.d0
2261 0 : do I=LA(N),0,-1
2262 0 : UP=UP*T+A(I,N)
2263 0 : if (I>=1) UP1=UP1*T+A(I,N)*I
2264 0 : if (I>=2) UP2=UP2*T+A(I,N)*I*(I-1)
2265 : end do
2266 0 : do I=LB(N),0,-1
2267 0 : DOWN=DOWN*T+B(I,N)
2268 0 : if (I>=1) DOWN1=DOWN1*T+B(I,N)*I
2269 0 : if (I>=2) DOWN2=DOWN2*T+B(I,N)*I*(I-1)
2270 : end do
2271 0 : X=log(T*UP/DOWN)
2272 0 : XDF=1.d0/T+UP1/UP-DOWN1/DOWN
2273 0 : XDFF=-1.d0/(T*T)+UP2/UP-pow2(UP1/UP)-DOWN2/DOWN+pow2(DOWN1/DOWN)
2274 : else
2275 0 : P=-1.d0/(.5d0+N) ! = -1/(1+\nu) = power index
2276 0 : T=pow(F,P) ! t - argument of the rational fraction
2277 0 : T1=P*T/F ! dt/df
2278 0 : T2=P*(P-1.d0)*T/(F*F) ! d^2 t / df^2
2279 0 : UP=0.d0
2280 0 : UP1=0.d0
2281 0 : UP2=0.d0
2282 0 : DOWN=0.d0
2283 0 : DOWN1=0.d0
2284 0 : DOWN2=0.d0
2285 0 : do I=6,0,-1
2286 0 : UP=UP*T+C(I,N)
2287 0 : if (I>=1) UP1=UP1*T+C(I,N)*I
2288 0 : if (I>=2) UP2=UP2*T+C(I,N)*I*(I-1)
2289 : end do
2290 0 : do I=LD(N),0,-1
2291 0 : DOWN=DOWN*T+D(I,N)
2292 0 : if (I>=1) DOWN1=DOWN1*T+D(I,N)*I
2293 0 : if (I>=2) DOWN2=DOWN2*T+D(I,N)*I*(I-1)
2294 : end do
2295 0 : R=UP/DOWN
2296 0 : R1=(UP1-UP*DOWN1/DOWN)/DOWN ! dR/dt
2297 0 : R2=(UP2-(2.0d0*UP1*DOWN1+UP*DOWN2)/DOWN+2.d0*UP*pow2(DOWN1/DOWN))/DOWN
2298 0 : X=R/T
2299 0 : RT=(R1-R/T)/T
2300 0 : XDF=T1*RT
2301 0 : XDFF=T2*RT+T1*T1*(R2-2d0*RT)/T
2302 : end if
2303 0 : return
2304 : end subroutine FERINV7
2305 :
2306 0 : subroutine BLIN9(TEMP,CHI, &
2307 : W0,W0DX,W0DT,W0DXX,W0DTT,W0DXT, &
2308 : W1,W1DX,W1DT,W1DXX,W1DTT,W1DXT, &
2309 : W2,W2DX,W2DT,W2DXX,W2DTT,W2DXT, &
2310 : W0XXX,W0XTT,W0XXT)
2311 : ! Version 21.01.10
2312 : ! Stems from BLIN8 v.24.12.08
2313 : ! Difference - smooth matching of different CHI ranges
2314 : ! Input: TEMP=T/mc^2; CHI=(\mu-mc^2)/T
2315 : ! Output: Wk - Fermi-Dirac integral of the order k+1/2
2316 : ! WkDX=dWk/dCHI, WkDT = dWk/dT, WkDXX=d^2 Wk / d CHI^2,
2317 : ! WkDTT=d^2 Wk / d T^2, WkDXT=d^2 Wk /dCHIdT,
2318 : ! W0XXX=d^3 W0 / d CHI^3, W0XTT=d^3 W0 /(d CHI d^2 T),
2319 : ! W0XXT=d^3 W0 /dCHI^2 dT
2320 : type(auto_diff_real_2var_order1), intent(in) :: TEMP,CHI
2321 : type(auto_diff_real_2var_order1), intent(out) :: W0,W0DX,W0DT,W0DXX,W0DTT,W0DXT
2322 : type(auto_diff_real_2var_order1), intent(out) :: W1,W1DX,W1DT,W1DXX,W1DTT,W1DXT
2323 : type(auto_diff_real_2var_order1), intent(out) :: W2,W2DX,W2DT,W2DXX,W2DTT,W2DXT
2324 : type(auto_diff_real_2var_order1), intent(out) :: W0XXX,W0XTT,W0XXT
2325 :
2326 : type(auto_diff_real_2var_order1) :: X1, X2, FP, FM
2327 : type(auto_diff_real_2var_order1) :: W0a,W0DXa,W0DTa,W0DXXa,W0DTTa,W0DXTa
2328 : type(auto_diff_real_2var_order1) :: W1a,W1DXa,W1DTa,W1DXXa,W1DTTa,W1DXTa
2329 : type(auto_diff_real_2var_order1) :: W2a,W2DXa,W2DTa,W2DXXa,W2DTTa,W2DXTa
2330 : type(auto_diff_real_2var_order1) :: W0XXXa,W0XTTa,W0XXTa
2331 : type(auto_diff_real_2var_order1) :: W0b,W0DXb,W0DTb,W0DXXb,W0DTTb,W0DXTb
2332 : type(auto_diff_real_2var_order1) :: W1b,W1DXb,W1DTb,W1DXXb,W1DTTb,W1DXTb
2333 : type(auto_diff_real_2var_order1) :: W2b,W2DXb,W2DTb,W2DXXb,W2DTTb,W2DXTb
2334 : type(auto_diff_real_2var_order1) :: W0XXXb,W0XTTb,W0XXTb
2335 :
2336 : type(auto_diff_real_2var_order1) :: CHI1
2337 : type(auto_diff_real_2var_order1) :: CHI2
2338 : type(auto_diff_real_2var_order1) :: XMAX
2339 : type(auto_diff_real_2var_order1) :: DCHI1
2340 : type(auto_diff_real_2var_order1) :: DCHI2
2341 : type(auto_diff_real_2var_order1) :: XSCAL1
2342 : type(auto_diff_real_2var_order1) :: XSCAL2
2343 :
2344 0 : CHI1=0.6d0
2345 0 : CHI2=14.d0
2346 0 : XMAX=30.d0
2347 0 : DCHI1=0.1d0
2348 0 : DCHI2=CHI2-CHI1-DCHI1
2349 0 : XSCAL1=XMAX/DCHI1
2350 0 : XSCAL2=XMAX/DCHI2
2351 :
2352 0 : X1=(CHI-CHI1)*XSCAL1
2353 0 : X2=(CHI-CHI2)*XSCAL2
2354 0 : if (X1<-XMAX) then
2355 : call BLIN9a(TEMP,CHI, &
2356 : W0,W0DX,W0DT,W0DXX,W0DTT,W0DXT, &
2357 : W1,W1DX,W1DT,W1DXX,W1DTT,W1DXT, &
2358 : W2,W2DX,W2DT,W2DXX,W2DTT,W2DXT, &
2359 0 : W0XXX,W0XTT,W0XXT)
2360 0 : elseif (X2<XMAX) then ! match two fits
2361 0 : if (X1<XMAX) then ! match fits "a" and "b"
2362 0 : call FERMI10(X1,XMAX,FP,FM)
2363 : call BLIN9a(TEMP,CHI, &
2364 : W0a,W0DXa,W0DTa,W0DXXa,W0DTTa,W0DXTa, &
2365 : W1a,W1DXa,W1DTa,W1DXXa,W1DTTa,W1DXTa, &
2366 : W2a,W2DXa,W2DTa,W2DXXa,W2DTTa,W2DXTa, &
2367 0 : W0XXXa,W0XTTa,W0XXTa)
2368 : call BLIN9b(TEMP,CHI, &
2369 : W0b,W0DXb,W0DTb,W0DXXb,W0DTTb,W0DXTb, &
2370 : W1b,W1DXb,W1DTb,W1DXXb,W1DTTb,W1DXTb, &
2371 : W2b,W2DXb,W2DTb,W2DXXb,W2DTTb,W2DXTb, &
2372 0 : W0XXXb,W0XTTb,W0XXTb)
2373 : else ! match fits "b" and "c"
2374 0 : call FERMI10(X2,XMAX,FP,FM)
2375 : call BLIN9b(TEMP,CHI, &
2376 : W0a,W0DXa,W0DTa,W0DXXa,W0DTTa,W0DXTa, &
2377 : W1a,W1DXa,W1DTa,W1DXXa,W1DTTa,W1DXTa, &
2378 : W2a,W2DXa,W2DTa,W2DXXa,W2DTTa,W2DXTa, &
2379 0 : W0XXXa,W0XTTa,W0XXTa)
2380 : call BLIN9c(TEMP,CHI, &
2381 : W0b,W0DXb,W0DTb,W0DXXb,W0DTTb,W0DXTb, &
2382 : W1b,W1DXb,W1DTb,W1DXXb,W1DTTb,W1DXTb, &
2383 : W2b,W2DXb,W2DTb,W2DXXb,W2DTTb,W2DXTb, &
2384 0 : W0XXXb,W0XTTb,W0XXTb)
2385 : end if
2386 0 : W0=W0a*FP+W0b*FM
2387 0 : W0DX=W0DXa*FP+W0DXb*FM !! +(W0a-W0b)*F1
2388 0 : W0DT=W0DTa*FP+W0DTb*FM
2389 0 : W0DXX=W0DXXa*FP+W0DXXb*FM !! +2.d0*(W0DXa-W0DXb)*F1+(W0a-W0b)*F2
2390 0 : W0DTT=W0DTTa*FP+W0DTTb*FM
2391 0 : W0DXT=W0DXTa*FP+W0DXTb*FM !! +(W0DTa-W0DTb)*F1
2392 0 : W0XXX=W0XXXa*FP+W0XXXb*FM !! +3.d0*(W0DXXa-W0DXXb)*F1+3.d0*(W0DXa-W0DXb)*F2+(W0a-W0b)*F3
2393 0 : W0XTT=W0XTTa*FP+W0XTTb*FM !! +(W0DTTa-W0DTTb)*F1
2394 0 : W0XXT=W0XXTa*FP+W0XXTb*FM !! +2.d0*(W0DXTa-W0DXTb)*F1+(W0DTa-W0DTb)*F2
2395 0 : W1=W1a*FP+W1b*FM
2396 0 : W1DX=W1DXa*FP+W1DXb*FM !! +(W1a-W1b)*F1
2397 0 : W1DT=W1DTa*FP+W1DTb*FM
2398 0 : W1DXX=W1DXXa*FP+W1DXXb*FM !! +2.d0*(W1DXa-W1DXb)*F1+(W1a-W1b)*F2
2399 0 : W1DTT=W1DTTa*FP+W1DTTb*FM
2400 0 : W1DXT=W1DXTa*FP+W1DXTb*FM !! +(W1DTa-W1DTb)*F1
2401 0 : W2=W2a*FP+W2b*FM
2402 0 : W2DX=W2DXa*FP+W2DXb*FM !! +(W2a-W2b)*F1
2403 0 : W2DT=W2DTa*FP+W2DTb*FM
2404 0 : W2DXX=W2DXXa*FP+W2DXXb*FM !! +2.d0*(W2DXa-W2DXb)*F1+(W2a-W2b)*F2
2405 0 : W2DTT=W2DTTa*FP+W2DTTb*FM
2406 0 : W2DXT=W2DXTa*FP+W2DXTb*FM !!
2407 : else
2408 : call BLIN9c(TEMP,CHI, &
2409 : W0,W0DX,W0DT,W0DXX,W0DTT,W0DXT, &
2410 : W1,W1DX,W1DT,W1DXX,W1DTT,W1DXT, &
2411 : W2,W2DX,W2DT,W2DXX,W2DTT,W2DXT, &
2412 0 : W0XXX,W0XTT,W0XXT)
2413 : end if
2414 0 : return
2415 : end subroutine BLIN9
2416 :
2417 0 : subroutine BLIN9a(TEMP,CHI, &
2418 : W0,W0DX,W0DT,W0DXX,W0DTT,W0DXT, &
2419 : W1,W1DX,W1DT,W1DXX,W1DTT,W1DXT, &
2420 : W2,W2DX,W2DT,W2DXX,W2DTT,W2DXT, &
2421 : W0XXX,W0XTT,W0XXT)
2422 : ! Version 19.01.10
2423 : ! First part of BILN9: small CHI. Stems from BLIN9 v.24.12.08
2424 : type(auto_diff_real_2var_order1), intent(in) :: TEMP,CHI
2425 : type(auto_diff_real_2var_order1), intent(out) :: W0,W0DX,W0DT,W0DXX,W0DTT,W0DXT
2426 : type(auto_diff_real_2var_order1), intent(out) :: W1,W1DX,W1DT,W1DXX,W1DTT,W1DXT
2427 : type(auto_diff_real_2var_order1), intent(out) :: W2,W2DX,W2DT,W2DXX,W2DTT,W2DXT
2428 : type(auto_diff_real_2var_order1), intent(out) :: W0XXX,W0XTT,W0XXT
2429 :
2430 : type(auto_diff_real_2var_order1) :: W,WDX,WDT,WDXX,WDTT,WDXT,WDXXX,WDXTT,WDXXT
2431 : type(auto_diff_real_2var_order1) :: ECHI, SQ, DN
2432 :
2433 : integer :: I, J, K
2434 :
2435 : ! The next three are really parameters but there isn't a clean way to initialize them
2436 : ! at declaration time. - Adam Jermyn 4/2/2020
2437 0 : real(dp) :: AC(5,0:2) ! read only after initialization
2438 0 : real(dp) :: AU(5,0:2) ! read only after initialization
2439 0 : real(dp) :: AA(5,0:2) ! read only after initialization
2440 :
2441 : AC(1:5,0) = [ &
2442 : 0.37045057d0, .41258437d0, &
2443 0 : 9.777982d-2, 5.3734153d-3, 3.8746281d-5] ! c_i^0
2444 : AC(1:5,1) = [ &
2445 : .39603109d0, .69468795d0, &
2446 0 : .22322760d0, 1.5262934d-2, 1.3081939d-4] ! c_i^1
2447 : AC(1:5,2) = [ &
2448 : .76934619d0, 1.7891437d0, &
2449 0 : .70754974d0, 5.6755672d-2, 5.5571480d-4] ! c_i^2
2450 : AU(1:5,0) = [ &
2451 : 0.43139881d0, 1.7597537d0, &
2452 0 : 4.1044654d0, 7.7467038d0, 13.457678d0] ! \chi_i^0
2453 : AU(1:5,1) = [ &
2454 : .81763176d0, 2.4723339d0, &
2455 0 : 5.1160061d0, 9.0441465d0, 15.049882d0] ! \chi_i^1
2456 : AU(1:5,2) = [ &
2457 : 1.2558461d0, 3.2070406d0, &
2458 0 : 6.1239082d0, 10.316126d0, 16.597079d0] ! \chi_i^2
2459 :
2460 0 : do J=0,2
2461 0 : do I=1,5
2462 0 : AA(I,J)=exp(-AU(I,J))
2463 : end do
2464 : end do
2465 :
2466 0 : do K=0,2
2467 0 : W=0.d0
2468 0 : WDX=0.d0
2469 0 : WDT=0.d0
2470 0 : WDXX=0.d0
2471 0 : WDTT=0.d0
2472 0 : WDXT=0.d0
2473 0 : WDXXX=0.d0
2474 0 : WDXTT=0.d0
2475 0 : WDXXT=0.d0
2476 0 : ECHI=exp(-CHI)
2477 0 : do I=1,5
2478 0 : SQ=sqrt(1.d0+AU(I,K)*TEMP/2.d0)
2479 0 : DN=AA(I,K)+ECHI ! e^{-\chi_i}+e^{-\chi})
2480 0 : W=W+AC(I,K)*SQ/DN
2481 0 : WDX=WDX+AC(I,K)*SQ/(DN*DN)
2482 0 : WDT=WDT+AC(I,K)*AU(I,K)/(SQ*DN)
2483 0 : WDXX=WDXX+AC(I,K)*SQ*(ECHI-AA(I,K))/pow3(DN)
2484 0 : WDTT=WDTT-AC(I,K)*AU(I,K)*AU(I,K)/(DN*SQ*SQ*SQ)
2485 0 : WDXT=WDXT+AC(I,K)*AU(I,K)/(SQ*DN*DN)
2486 : WDXXX=WDXXX+AC(I,K)*SQ* &
2487 0 : (ECHI*ECHI-4.d0*ECHI*AA(I,K)+AA(I,K)*AA(I,K))/(DN*DN*DN*DN)
2488 0 : WDXTT=WDXTT-AC(I,K)*AU(I,K)*AU(I,K)/(DN*DN*SQ*SQ*SQ)
2489 0 : WDXXT=WDXXT+AC(I,K)*AU(I,K)*(ECHI-AA(I,K))/(SQ*DN*DN*DN)
2490 : end do
2491 0 : WDX=WDX*ECHI
2492 0 : WDT=0.25d0*WDT
2493 0 : WDXX=WDXX*ECHI
2494 0 : WDTT=0.0625d0*WDTT
2495 0 : WDXT=0.25d0*WDXT*ECHI
2496 0 : WDXXX=WDXXX*ECHI
2497 0 : WDXTT=0.0625d0*WDXTT*ECHI
2498 0 : WDXXT=0.25d0*WDXXT*ECHI
2499 0 : if (K==0) then
2500 0 : W0=W
2501 0 : W0DX=WDX
2502 0 : W0DT=WDT
2503 0 : W0DXX=WDXX
2504 0 : W0DTT=WDTT
2505 0 : W0DXT=WDXT
2506 0 : W0XXX=WDXXX
2507 0 : W0XTT=WDXTT
2508 0 : W0XXT=WDXXT
2509 0 : elseif (K==1) then
2510 0 : W1=W
2511 0 : W1DX=WDX
2512 0 : W1DT=WDT
2513 0 : W1DXX=WDXX
2514 0 : W1DTT=WDTT
2515 0 : W1DXT=WDXT
2516 : else
2517 0 : W2=W
2518 0 : W2DX=WDX
2519 0 : W2DT=WDT
2520 0 : W2DXX=WDXX
2521 0 : W2DTT=WDTT
2522 0 : W2DXT=WDXT
2523 : end if
2524 : end do ! next K
2525 0 : return
2526 : end subroutine BLIN9a
2527 :
2528 0 : subroutine BLIN9b(TEMP,CHI, &
2529 : W0,W0DX,W0DT,W0DXX,W0DTT,W0DXT, &
2530 : W1,W1DX,W1DT,W1DXX,W1DTT,W1DXT, &
2531 : W2,W2DX,W2DT,W2DXX,W2DTT,W2DXT, &
2532 : W0XXX,W0XTT,W0XXT)
2533 : ! Version 19.01.10
2534 : ! Second part of BILN9: intermediate CHI. Stems from BLIN8 v.24.12.08
2535 : type(auto_diff_real_2var_order1), intent(in) :: TEMP
2536 : type(auto_diff_real_2var_order1) :: CHI ! can be modified, not sure if this is an intended side-effect
2537 : type(auto_diff_real_2var_order1), intent(out) :: W0,W0DX,W0DT,W0DXX,W0DTT,W0DXT
2538 : type(auto_diff_real_2var_order1), intent(out) :: W1,W1DX,W1DT,W1DXX,W1DTT,W1DXT
2539 : type(auto_diff_real_2var_order1), intent(out) :: W2,W2DX,W2DT,W2DXX,W2DTT,W2DXT
2540 : type(auto_diff_real_2var_order1), intent(out) :: W0XXX,W0XTT,W0XXT
2541 :
2542 : type(auto_diff_real_2var_order1) :: W,WDX,WDT,WDXX,WDTT,WDXT,WDXXX,WDXTT,WDXXT
2543 : type(auto_diff_real_2var_order1) :: SQCHI, CE, ECHI, DE, D, XICHI, DXI
2544 : type(auto_diff_real_2var_order1) :: H, HX, HDX, HXX, HDXX, HT, HDT, HDTT
2545 : type(auto_diff_real_2var_order1) :: HTX, HDXT, HDXXT, HDXTT, HXXX, HDXXX
2546 : type(auto_diff_real_2var_order1) :: V, VX, VDX, VT, VDT, VXX, VDXX, VDXXX
2547 : type(auto_diff_real_2var_order1) :: VXXT, VDTT, VXT, VDXT, VDXXT, VDXTT
2548 :
2549 : integer :: I, J, K
2550 : real(dp), parameter :: AX(5) = &
2551 : [7.265351d-2, .2694608d0, &
2552 : .533122d0, .7868801d0, .9569313d0] ! x_i
2553 : real(dp), parameter :: AXI(5) = &
2554 : [.26356032d0, 1.4134031d0, &
2555 : 3.5964258d0, 7.0858100d0, 12.640801d0] ! \xi_i
2556 : real(dp), parameter :: AH(5) = &
2557 : [3.818735d-2, .1256732d0, &
2558 : .1986308d0, .1976334d0, .1065420d0] ! H_i
2559 : real(dp), parameter :: AV(5) = &
2560 : [.29505869d0, .32064856d0, 7.3915570d-2, &
2561 : 3.6087389d-3, 2.3369894d-5] ! \bar{V}_i
2562 : real(dp), parameter :: EPS=1.d-3
2563 :
2564 0 : if (CHI<EPS) CHI = EPS !call mesa_error(__FILE__,__LINE__,'BLIN9b: CHI is too small')
2565 0 : do K=0,2
2566 0 : W=0.d0
2567 0 : WDX=0.d0
2568 0 : WDT=0.d0
2569 0 : WDXX=0.d0
2570 0 : WDTT=0.d0
2571 0 : WDXT=0.d0
2572 0 : WDXXX=0.d0
2573 0 : WDXTT=0.d0
2574 0 : WDXXT=0.d0
2575 0 : SQCHI=sqrt(CHI)
2576 0 : do I=1,5
2577 0 : CE=AX(I)-1.d0
2578 0 : ECHI=exp(CE*CHI)
2579 0 : DE=1.d0+ECHI
2580 0 : D=1.d0+AX(I)*CHI*TEMP/2.d0
2581 0 : H=pow(CHI,K+1d0)*SQCHI*sqrt(D)/DE
2582 0 : HX=(K+1.5d0)/CHI+.25d0*AX(I)*TEMP/D-ECHI*CE/DE
2583 0 : HDX=H*HX
2584 0 : HXX=(K+1.5d0)/(CHI*CHI)+.125d0*pow2(AX(I)*TEMP/D)+ECHI*pow2(CE/DE)
2585 0 : HDXX=HDX*HX-H*HXX
2586 0 : HT=.25d0*AX(I)*CHI/D
2587 0 : HDT=H*HT
2588 0 : HDTT=-H*HT*HT
2589 0 : HTX=1.d0/CHI-.5d0*AX(I)*TEMP/D
2590 0 : HDXT=HDX*HT+HDT*HTX
2591 : HDXXT=HDXX*HT+HDX*HT*HTX+HDXT*HTX &
2592 0 : + HDT*(.25d0*pow2(AX(I)*TEMP/D)-1.d0/(CHI*CHI))
2593 : HDXTT=HDXT*HT-HDX*.125d0*pow2(AX(I)*CHI/D)+HDTT*HTX &
2594 0 : + HDT*.5d0*AX(I)*(TEMP*.5d0*AX(I)*CHI/(D*D)-1.d0/D)
2595 : HXXX=(2d0*K+3d0)/pow3(CHI)+.125d0*pow3(AX(I)*TEMP/D) &
2596 0 : - ECHI*(1.d0-ECHI)*pow3(CE/DE)
2597 0 : HDXXX=HDXX*HX-2.d0*HDX*HXX+H*HXXX
2598 0 : XICHI=AXI(I)+CHI
2599 0 : DXI=1.d0+XICHI*TEMP/2.d0
2600 0 : V=pow(XICHI,1d0*K)*sqrt(XICHI*DXI)
2601 0 : VX=(K+.5d0)/XICHI+.25d0*TEMP/DXI
2602 0 : VDX=V*VX
2603 0 : VT=.25d0*XICHI/DXI
2604 0 : VDT=V*VT
2605 0 : VXX=(K+0.5d0)/(XICHI*XICHI)+.125d0*pow2(TEMP/DXI)
2606 0 : VDXX=VDX*VX-V*VXX
2607 : VDXXX=VDXX*VX-2d0*VDX*VXX &
2608 0 : + V*((2d0*K+1d0)/pow3(XICHI)+.125d0*pow3(TEMP/DXI))
2609 0 : VXXT=(1.d0-0.5d0*TEMP*XICHI/DXI)/DXI
2610 0 : VDTT=-V*VT*VT
2611 0 : VXT=1.d0/XICHI-0.5d0*TEMP/DXI
2612 0 : VDXT=VDT*VXT+VDX*VT
2613 0 : VDXXT=VDXT*VX+VDX*.25d0*VXXT-VDT*VXX-V*.25d0*TEMP/DXI*VXXT
2614 : VDXTT=VDTT*VXT-VDT*0.5d0*VXXT+VDXT*VT &
2615 0 : - VDX*.125d0*pow2(XICHI/DXI)
2616 0 : W=W+AH(I)*pow(AX(I),K)*H+AV(I)*V
2617 0 : WDX=WDX+AH(I)*pow(AX(I),K)*HDX+AV(I)*VDX
2618 0 : WDT=WDT+AH(I)*pow(AX(I),K)*HDT+AV(I)*VDT
2619 0 : WDXX=WDXX+AH(I)*pow(AX(I),K)*HDXX+AV(I)*VDXX
2620 0 : WDTT=WDTT+AH(I)*pow(AX(I),K)*HDTT+AV(I)*VDTT
2621 0 : WDXT=WDXT+AH(I)*pow(AX(I),K)*HDXT+AV(I)*VDXT
2622 0 : WDXXX=WDXXX+AH(I)*pow(AX(I),K)*HDXXX+AV(I)*VDXXX
2623 0 : WDXTT=WDXTT+AH(I)*pow(AX(I),K)*HDXTT+AV(I)*VDXTT
2624 0 : WDXXT=WDXXT+AH(I)*pow(AX(I),K)*HDXXT+AV(I)*VDXXT
2625 : end do
2626 0 : if (K==0) then
2627 0 : W0=W
2628 0 : W0DX=WDX
2629 0 : W0DT=WDT
2630 0 : W0DXX=WDXX
2631 0 : W0DTT=WDTT
2632 0 : W0DXT=WDXT
2633 0 : W0XXX=WDXXX
2634 0 : W0XTT=WDXTT
2635 0 : W0XXT=WDXXT
2636 0 : elseif (K==1) then
2637 0 : W1=W
2638 0 : W1DX=WDX
2639 0 : W1DT=WDT
2640 0 : W1DXX=WDXX
2641 0 : W1DTT=WDTT
2642 0 : W1DXT=WDXT
2643 : else
2644 0 : W2=W
2645 0 : W2DX=WDX
2646 0 : W2DT=WDT
2647 0 : W2DXX=WDXX
2648 0 : W2DTT=WDTT
2649 0 : W2DXT=WDXT
2650 : end if
2651 : end do ! next K
2652 0 : return
2653 : end subroutine BLIN9b
2654 :
2655 0 : subroutine BLIN9c(TEMP,CHI, &
2656 : W0,W0DX,W0DT,W0DXX,W0DTT,W0DXT, &
2657 : W1,W1DX,W1DT,W1DXX,W1DTT,W1DXT, &
2658 : W2,W2DX,W2DT,W2DXX,W2DTT,W2DXT, &
2659 : W0XXX,W0XTT,W0XXT)
2660 : ! Version 19.01.10
2661 : ! Third part of BILN9: large CHI. Stems from BLIN8 v.24.12.08
2662 : type(auto_diff_real_2var_order1), intent(in) :: CHI, TEMP
2663 : ! type(auto_diff_real_2var_order1) :: CHI ! can be modified, not sure if this is an intended side-effect
2664 : type(auto_diff_real_2var_order1), intent(out) :: W0,W0DX,W0DT,W0DXX,W0DTT,W0DXT
2665 : type(auto_diff_real_2var_order1), intent(out) :: W1,W1DX,W1DT,W1DXX,W1DTT,W1DXT
2666 : type(auto_diff_real_2var_order1), intent(out) :: W2,W2DX,W2DT,W2DXX,W2DTT,W2DXT
2667 : type(auto_diff_real_2var_order1), intent(out) :: W0XXX,W0XTT,W0XXT
2668 :
2669 : type(auto_diff_real_2var_order1), dimension(0:2) :: AM,AMDX,AMDT,AMDXX,AMDTT,AMDXT
2670 : type(auto_diff_real_2var_order1) :: CNU, CHINU, F, FDX, FDXX, FDXXX, C, D, DM
2671 : type(auto_diff_real_2var_order1) :: W, WDX, WDT, WDXX, WDXT, WDTT, WDXXX, WDXXT, WDXTT
2672 : type(auto_diff_real_2var_order1) :: R, RX, RDX, RDT, RXX, RDXX, RDTT, RXT, RDXT
2673 : type(auto_diff_real_2var_order1) :: RXXX, RDXXX, RXTT, RDXTT, RXXT, RDXXT
2674 : type(auto_diff_real_2var_order1) :: FMX1, FMX2, FMX, CKM
2675 : type(auto_diff_real_2var_order1) :: FMT1, FMT2, FMT, FMXX, FMXXX, FMTT
2676 : type(auto_diff_real_2var_order1) :: AMDXXX, FMT1DX, FMT2DX, FMXT, FMTTX
2677 : type(auto_diff_real_2var_order1) :: AMDXTT, FMX1DT, FMX2DT, FMXXT, AMDXXT
2678 : type(auto_diff_real_2var_order1) :: SQ2T, SQ2T3
2679 : type(auto_diff_real_2var_order1) :: A, ADX, ADT, ADXX, ADTT, ADXT, ADXTT, ADXXT
2680 : type(auto_diff_real_2var_order1) :: XT1, Aln, FJ0, ASQ3, ASQ3DX, BXT, BXXT
2681 : type(auto_diff_real_2var_order1) :: FJ0DX, FJ0DT, FJ0DXX, FJ0DTT, FJ0DXT
2682 : type(auto_diff_real_2var_order1) :: FJ0XXX, FJ0XXT, FJ0XTT
2683 : type(auto_diff_real_2var_order1) :: FJ1, FJ1DX, FJ1DT, FJ1DXX, FJ1DXT, FJ1DTT
2684 : type(auto_diff_real_2var_order1) :: FJ2, FJ2DX, FJ2DT, FJ2DXX, FJ2DXT, FJ2DTT
2685 :
2686 : integer :: J, K
2687 : real(dp), parameter :: PI26=PI*PI/6.d0
2688 :
2689 0 : if (CHI*TEMP<.1d0) then
2690 0 : do K=0,2
2691 0 : W=0.d0
2692 0 : WDX=0.d0
2693 0 : WDT=0.d0
2694 0 : WDXX=0.d0
2695 0 : WDTT=0.d0
2696 0 : WDXT=0.d0
2697 0 : WDXXX=0.d0
2698 0 : WDXTT=0.d0
2699 0 : WDXXT=0.d0
2700 0 : do J=0,4 ! for nonrel.Fermi integrals from k+1/2 to k+4.5
2701 0 : CNU=K+J+0.5d0 ! nonrelativistic Fermi integral index \nu
2702 0 : CHINU=pow(CHI,1d0*(K+J))*sqrt(CHI) ! \chi^\nu
2703 : F=CHINU*(CHI/(CNU+1.d0)+PI26*CNU/CHI & ! nonrel.Fermi
2704 0 : + .7d0*PI26*PI26*CNU*(CNU-1.d0)*(CNU-2.d0)/pow3(CHI))
2705 : FDX=CHINU*(1d0+PI26*CNU*(CNU-1.d0)/pow2(CHI) &
2706 0 : +.7d0*PI26*PI26*CNU*(CNU-1.d0)*(CNU-2.d0)*(CNU-3.d0)/pow4(CHI))
2707 : FDXX=CHINU/CHI*CNU*(1.d0+PI26*(CNU-1.d0)*(CNU-2.d0)/pow2(CHI) &
2708 0 : +.7d0*PI26*PI26*(CNU-1.d0)*(CNU-2.d0)*(CNU-3.d0)*(CNU-4.d0)/pow4(CHI))
2709 : FDXXX=CHINU/pow2(CHI)*CNU*(CNU-1.d0)* &
2710 : (1.d0+PI26*(CNU-2.d0)*(CNU-3.d0)/pow2(CHI) &
2711 0 : +.7d0*PI26*PI26*(CNU-2.d0)*(CNU-3.d0)*(CNU-4.d0)*(CNU-5.d0)/pow4(CHI))
2712 0 : if (J==0) then
2713 0 : W=F
2714 0 : WDX=FDX
2715 0 : WDXX=FDXX
2716 0 : WDXXX=FDXXX
2717 0 : elseif (J==1) then
2718 0 : C=.25d0*TEMP
2719 0 : W=W+C*F ! Fermi-Dirac, expressed through Fermi
2720 0 : WDX=WDX+C*FDX
2721 0 : WDXX=WDXX+C*FDXX
2722 0 : WDT=F/4.d0
2723 0 : WDXT=FDX/4.d0
2724 0 : WDTT=0.d0
2725 0 : WDXXX=WDXXX+C*FDXXX
2726 0 : WDXXT=FDXX/4.d0
2727 0 : WDXTT=0.d0
2728 : else
2729 0 : C=-C/(1d0*J)*(2d0*J-3d0)/4d0*TEMP
2730 0 : W=W+C*F
2731 0 : WDX=WDX+C*FDX
2732 0 : WDT=WDT+C*(1d0*J)/TEMP*F
2733 0 : WDXX=WDXX+C*FDXX
2734 0 : WDTT=WDTT+C*(1d0*J)*(1d0*(J-1))/pow2(TEMP)*F
2735 0 : WDXT=WDXT+C*(1d0*J)/TEMP*FDX
2736 0 : WDXXX=WDXXX+C*FDXXX
2737 0 : WDXTT=WDXTT+C*(1d0*J)*(1d0*(J-1))/pow2(TEMP)*FDX
2738 0 : WDXXT=WDXXT+C*(1d0*J)/TEMP*FDXX
2739 : end if
2740 : end do ! next J
2741 0 : if (K==0) then
2742 0 : W0=W
2743 0 : W0DX=WDX
2744 0 : W0DT=WDT
2745 0 : W0DXX=WDXX
2746 0 : W0DTT=WDTT
2747 0 : W0DXT=WDXT
2748 0 : W0XXX=WDXXX
2749 0 : W0XTT=WDXTT
2750 0 : W0XXT=WDXXT
2751 0 : elseif (K==1) then
2752 0 : W1=W
2753 0 : W1DX=WDX
2754 0 : W1DT=WDT
2755 0 : W1DXX=WDXX
2756 0 : W1DTT=WDTT
2757 0 : W1DXT=WDXT
2758 : else
2759 0 : W2=W
2760 0 : W2DX=WDX
2761 0 : W2DT=WDT
2762 0 : W2DXX=WDXX
2763 0 : W2DTT=WDTT
2764 0 : W2DXT=WDXT
2765 : end if
2766 : end do ! next K
2767 : ! ---------------------------------------------------------------- *
2768 : else ! CHI > 14, CHI*TEMP > 0.1: general high-\chi expansion
2769 0 : D=1.d0+CHI*TEMP/2.d0
2770 0 : R=sqrt(CHI*D)
2771 0 : RX=.5d0/CHI+.25d0*TEMP/D
2772 0 : RDX=R*RX
2773 0 : RDT=.25d0*pow2(CHI)/R
2774 0 : RXX=-.5d0/pow2(CHI)-.125d0*pow2(TEMP/D)
2775 0 : RDXX=RDX*RX+R*RXX
2776 0 : RDTT=-.25d0*RDT*CHI/D
2777 0 : RXT=.25d0/D-.125d0*CHI*TEMP/(D*D)
2778 0 : RDXT=RDT*RX+R*RXT
2779 0 : RXXX=1.d0/pow3(CHI)+.125d0*pow3(TEMP/D)
2780 0 : RDXXX=RDXX*RX+2.d0*RDX*RXX+R*RXXX
2781 0 : RXTT=-.25d0/(D*D)*CHI+.125d0*pow2(CHI)*TEMP/pow3(D)
2782 0 : RDXTT=RDTT*RX+2.d0*RDT*RXT+R*RXTT
2783 0 : RXXT=-RXT*TEMP/D
2784 0 : RDXXT=RDXT*RX+RDX*RXT+RDT*RXX+R*RXXT
2785 0 : do K=0,2
2786 0 : DM=K+.5d0+(K+1.d0)*CHI*TEMP/2.d0
2787 0 : AM(K)=pow(CHI,1d0*K)*DM/R
2788 0 : FMX1=.5d0*(K+1.d0)*TEMP/DM
2789 0 : FMX2=.25d0*TEMP/D
2790 0 : FMX=(K-.5d0)/CHI+FMX1-FMX2
2791 0 : AMDX(K)=AM(K)*FMX
2792 0 : CKM=.5d0*(K+1.d0)/DM
2793 0 : FMT1=CKM*CHI
2794 0 : FMT2=.25d0*CHI/D
2795 0 : FMT=FMT1-FMT2
2796 0 : AMDT(K)=AM(K)*FMT
2797 0 : FMXX=-(K-.5d0)/pow2(CHI)-FMX1*FMX1+2.d0*FMX2*FMX2
2798 0 : AMDXX(K)=AMDX(K)*FMX+AM(K)*FMXX
2799 0 : FMTT=2.d0*FMT2*FMT2-FMT1*FMT1
2800 0 : AMDTT(K)=AMDT(K)*FMT+AM(K)*FMTT
2801 : AMDXT(K)=AMDX(K)*FMT+AM(K)*(CKM*(1.d0-CKM*CHI*TEMP) &
2802 0 : - .25d0/D+.125d0*CHI*TEMP/(D*D))
2803 0 : if (K==0) then
2804 0 : FMXXX=(2d0*K-1d0)/pow3(CHI)+2.d0*pow3(FMX1)-8.d0*pow3(FMX2)
2805 0 : AMDXXX=AMDXX(K)*FMX+2.d0*AMDX(K)*FMXX+AM(K)*FMXXX
2806 0 : FMT1DX=CKM-TEMP*CHI*CKM*CKM
2807 0 : FMT2DX=(.25d0-CHI*TEMP*.125d0/D)/D
2808 0 : FMXT=FMT1DX-FMT2DX
2809 0 : FMTTX=4.d0*FMT2*FMT2DX-2.d0*FMT1*FMT1DX
2810 0 : AMDXTT=AMDXT(K)*FMT+AMDT(K)*FMXT+AMDX(K)*FMTT+AM(K)*FMTTX
2811 0 : FMX1DT=CKM-CHI*TEMP*CKM*CKM
2812 0 : FMX2DT=.25d0/D*(1.d0-.5d0*CHI*TEMP/D)
2813 0 : FMXXT=4.d0*FMX2*FMX2DT-2.d0*FMX1*FMX1DT
2814 0 : AMDXXT=AMDXT(K)*FMX+AMDX(K)*FMXT+AMDT(K)*FMXX+AM(K)*FMXXT
2815 : end if
2816 : end do
2817 0 : SQ2T=sqrt(2.d0*TEMP)
2818 0 : SQ2T3=SQ2T*SQ2T*SQ2T
2819 0 : A=1.d0+CHI*TEMP+SQ2T*R
2820 0 : ADX=TEMP+SQ2T*RDX
2821 0 : ADT=CHI+R/SQ2T+SQ2T*RDT
2822 0 : ADXX=SQ2T*RDXX
2823 0 : ADTT=-R/SQ2T3+2.d0/SQ2T*RDT+SQ2T*RDTT
2824 0 : ADXT=1.d0+RDX/SQ2T+SQ2T*RDXT
2825 0 : ADXTT=-RDX/SQ2T3+2.d0/SQ2T*RDXT+SQ2T*RDXTT
2826 0 : ADXXT=RDXX/SQ2T+SQ2T*RDXXT
2827 0 : XT1=CHI+1.d0/TEMP
2828 0 : Aln=log(A)
2829 0 : FJ0=.5d0*XT1*R-Aln/SQ2T3
2830 0 : ASQ3=A*SQ2T3
2831 0 : ASQ3DX=ADX*SQ2T3
2832 0 : FJ0DX=.5d0*(R+XT1*RDX)-ADX/ASQ3
2833 : FJ0DT=.5d0*(XT1*RDT-R/pow2(TEMP))-ADT/ASQ3 &
2834 0 : + .75d0/(TEMP*TEMP*SQ2T)*Aln
2835 0 : FJ0DXX=RDX+.5d0*XT1*RDXX+pow2(ADX/A)/SQ2T3-ADXX/ASQ3
2836 : FJ0DTT=R/pow3(TEMP)-RDT/pow2(TEMP)+.5d0*XT1*RDTT &
2837 : + 3.d0/(ASQ3*TEMP)*ADT &
2838 0 : + pow2(ADT/A)/SQ2T3-ADTT/ASQ3-1.875d0/(TEMP*TEMP*TEMP*SQ2T)*Aln
2839 0 : BXT=1.5d0/TEMP*ADX+ADX*ADT/A-ADXT
2840 : BXXT=1.5d0/TEMP*ADXX+(ADXX*ADT+ADX*ADXT)/A &
2841 0 : - pow2(ADX/A)*ADT-ADXXT
2842 0 : FJ0DXT=.5d0*(RDT-RDX/pow2(TEMP)+XT1*RDXT)+BXT/ASQ3
2843 : FJ0XXX=RDXX*1.5d0+.5d0*XT1*RDXXX &
2844 : +(2.d0*ADX*(ADXX/A-pow2(ADX/A)) &
2845 0 : - SQ2T*RDXXX+ADXX/ASQ3*ASQ3DX)/ASQ3
2846 : FJ0XTT=RDX/pow3(TEMP)-RDXT/pow2(TEMP)+.5d0*(RDTT+XT1*RDXTT) &
2847 : + 3.d0/TEMP*(ADXT-ADT/ASQ3*ASQ3DX)/ASQ3 &
2848 : + (2.d0*ADT*(ADXT/A-ADT*ADX/(A*A)) &
2849 0 : - ADXTT+ADTT*ASQ3DX/ASQ3)/ASQ3-1.875d0/(TEMP*TEMP*TEMP*SQ2T)*ADX/A
2850 : FJ0XXT=.5d0*(RDXT-RDXX/pow2(TEMP)+RDXT+XT1*RDXXT) &
2851 0 : +(BXXT-BXT*ASQ3DX/ASQ3)/ASQ3
2852 0 : W0=FJ0+PI26*AM(0)
2853 0 : W0DX=FJ0DX+PI26*AMDX(0)
2854 0 : W0DT=FJ0DT+PI26*AMDT(0)
2855 0 : W0DXX=FJ0DXX+PI26*AMDXX(0)
2856 0 : W0DTT=FJ0DTT+PI26*AMDTT(0)
2857 0 : W0DXT=FJ0DXT+PI26*AMDXT(0)
2858 0 : W0XXX=FJ0XXX+PI26*AMDXXX
2859 0 : W0XTT=FJ0XTT+PI26*AMDXTT
2860 0 : W0XXT=FJ0XXT+PI26*AMDXXT
2861 0 : FJ1=(R*R*R/1.5d0-FJ0)/TEMP
2862 0 : FJ1DX=(2.d0*R*R*RDX-FJ0DX)/TEMP
2863 0 : FJ1DT=(2.d0*R*R*RDT-FJ0DT-FJ1)/TEMP
2864 0 : FJ1DXX=(4.d0*R*RDX*RDX+2.d0*R*R*RDXX-FJ0DXX)/TEMP
2865 0 : FJ1DTT=(4.d0*R*RDT*RDT+2.d0*R*R*RDTT-FJ0DTT-2.d0*FJ1DT)/TEMP
2866 0 : FJ1DXT=(4.d0*R*RDX*RDT+2.d0*R*R*RDXT-FJ0DXT-FJ1DX)/TEMP
2867 0 : W1=FJ1+PI26*AM(1)
2868 0 : W1DX=FJ1DX+PI26*AMDX(1)
2869 0 : W1DT=FJ1DT+PI26*AMDT(1)
2870 0 : W1DXX=FJ1DXX+PI26*AMDXX(1)
2871 0 : W1DTT=FJ1DTT+PI26*AMDTT(1)
2872 0 : W1DXT=FJ1DXT+PI26*AMDXT(1)
2873 0 : FJ2=(.5d0*CHI*R*R*R-1.25d0*FJ1)/TEMP
2874 0 : FJ2DX=(.5d0*R*R*R+1.5d0*CHI*R*R*RDX-1.25d0*FJ1DX)/TEMP
2875 0 : FJ2DT=(1.5d0*CHI*R*R*RDT-1.25d0*FJ1DT-FJ2)/TEMP
2876 0 : FJ2DXX=(3.d0*R*RDX*(R+CHI*RDX)+1.5d0*CHI*R*R*RDXX-1.25d0*FJ1DXX)/TEMP
2877 0 : FJ2DTT=(3.d0*CHI*R*(RDT*RDT+.5d0*R*RDTT)-1.25d0*FJ1DTT-2.d0*FJ2DT)/TEMP
2878 0 : FJ2DXT=(1.5d0*R*RDT*(R+2.d0*CHI*RDX)+1.5d0*CHI*R*R*RDXT-1.25d0*FJ1DXT-FJ2DX)/TEMP
2879 0 : W2=FJ2+PI26*AM(2)
2880 0 : W2DX=FJ2DX+PI26*AMDX(2)
2881 0 : W2DT=FJ2DT+PI26*AMDT(2)
2882 0 : W2DXX=FJ2DXX+PI26*AMDXX(2)
2883 0 : W2DTT=FJ2DTT+PI26*AMDTT(2)
2884 0 : W2DXT=FJ2DXT+PI26*AMDXT(2)
2885 : end if
2886 0 : return
2887 : end subroutine BLIN9c
2888 :
2889 0 : subroutine CHEMFIT(DENS,TEMP,CHI)
2890 : ! Version 07.06.07
2891 : ! This is merely an interface to CHEMFIT7 for compatibility purposes.
2892 : ! Input: DENS - electron density [a.u.=6.7483346e24 cm^{-3}],
2893 : ! TEMP - temperature [a.u.=2Ryd=3.1577e5 K]
2894 : ! Output: CHI=\mu/TEMP, where \mu - electron chem.pot.w/o rest-energy
2895 : type(auto_diff_real_2var_order1), intent(in) :: DENS, TEMP
2896 : type(auto_diff_real_2var_order1), intent(out) :: CHI
2897 :
2898 : type(auto_diff_real_2var_order1) :: DENR,TEMR,CMU1,CMUDENR,CMUDT,CMUDTT
2899 :
2900 0 : DENR=DENS/2.5733806d6 ! n_e in rel.un.=\lambda_{Compton}^{-3}
2901 0 : TEMR=TEMP/1.8778865d4 ! T in rel.un.=(mc^2/k)=5.93e9 K
2902 0 : call CHEMFIT7(DENR,TEMR,CHI,CMU1,0,CMUDENR,CMUDT,CMUDTT)
2903 0 : return
2904 : end subroutine CHEMFIT
2905 :
2906 0 : subroutine CHEMFIT7(DENR,TEMR,CHI,CMU1,KDERIV, &
2907 : CMUDENR,CMUDT,CMUDTT)
2908 : ! Version 28.05.07
2909 : ! corrected 08.08.11
2910 : ! cosmetic change 16.05.13
2911 : ! Fit to the chemical potential of free electron gas described in:
2912 : ! G.Chabrier & A.Y.Potekhin, Phys.Rev.E 58, 4941 (1998)
2913 : ! Stems from CHEMFIT v.10.10.96. The main difference - derivatives.
2914 : ! All quantities are by default in relativistic units
2915 : ! Input: DENR - electron density, TEMR - temperature
2916 : ! KDERIV=0 if the derivatives are not required
2917 : ! Output: CHI=CMU1/TEMR, where CMU1 = \mu-1 - chem.pot.w/o rest-energy
2918 : ! CMUDENR = (d\mu/d n_e)_T
2919 : ! CMUDT = (d\mu/dT)_V
2920 : ! CMUDTT = (d^2\mu/dT^2)_V
2921 : ! CMUDENR,CMUDT, and CMUDTT =0 on output, if KREDIV=0
2922 : type(auto_diff_real_2var_order1), intent(in) :: DENR,TEMR
2923 : integer, intent(in) :: KDERIV
2924 : type(auto_diff_real_2var_order1), intent(out) :: CHI,CMU1,CMUDENR,CMUDT,CMUDTT
2925 :
2926 : type(auto_diff_real_2var_order1) :: PF0, TF, THETA, THETA32, Q2, T1, U3, THETAC, THETAG
2927 : type(auto_diff_real_2var_order1) :: D3, Q3, Q1, SQT, G, H, CT, F, X, XDF, XDFF
2928 : type(auto_diff_real_2var_order1) :: THETA52, CHIDY, CHIDYY
2929 : type(auto_diff_real_2var_order1) :: Q1D, Q1DD, Q2D, Q2DD, U3D, D3D, D3DD, Q3D, Q3DD
2930 : type(auto_diff_real_2var_order1) :: GDY, GDT, GDYY, GDTT, GDYT, GH
2931 : type(auto_diff_real_2var_order1) :: HDY, HDT, HDYY, HDTT, HDYT
2932 : type(auto_diff_real_2var_order1) :: CTY, CTT, CTDY, CTDT, CTDYY, CTDTT, CTDYT
2933 : type(auto_diff_real_2var_order1) :: CHIDT, CHIDTT, CHIDYT
2934 :
2935 : real(dp), parameter :: PARA=1.612d0
2936 : real(dp), parameter :: PARB=6.192d0
2937 : real(dp), parameter :: PARC=0.0944d0
2938 : real(dp), parameter :: PARF=5.535d0
2939 : real(dp), parameter :: PARG=0.698d0
2940 :
2941 0 : PF0=pow(29.6088132d0*DENR,1d0/3d0) ! Classical Fermi momentum
2942 0 : if (PF0>1.d-4) then
2943 0 : TF=sqrt(1.d0+PF0*PF0)-1.d0 ! Fermi temperature
2944 : else
2945 0 : TF=.5d0*PF0*PF0
2946 : end if
2947 0 : THETA=TEMR/TF
2948 0 : THETA32=THETA*sqrt(THETA)
2949 0 : Q2=12.d0+8.d0/THETA32
2950 0 : T1=exp(-THETA) ! former ('96) 1/T
2951 0 : U3=T1*T1+PARA
2952 0 : THETAC=pow(THETA,PARC)
2953 0 : THETAG=pow(THETA,PARG)
2954 0 : D3=PARB*THETAC*T1*T1+PARF*THETAG
2955 0 : Q3=1.365568127d0-U3/D3 ! 1.365...=2/\pi^{1/3}
2956 0 : if (THETA>1.d-5) then
2957 0 : Q1=1.5d0*T1/(1.d0-T1)
2958 : else
2959 0 : Q1=1.5d0/THETA
2960 : end if
2961 0 : SQT=sqrt(TEMR)
2962 0 : G=(1.d0+Q2*TEMR*Q3+Q1*SQT)*TEMR
2963 0 : H=(1.d0+.5d0*TEMR/THETA)*(1.d0+Q2*TEMR)
2964 0 : CT=1.d0+G/H
2965 0 : F=(2d0/3d0)/THETA32
2966 0 : call FERINV7(F,1,X,XDF,XDFF)
2967 : CHI=X & ! non-relativistic result
2968 0 : - 1.5d0*log(CT) ! Relativistic fit
2969 0 : CMU1=TEMR*CHI ! Fit to chemical potential w/o mc^2
2970 0 : if (KDERIV==0) then ! DISMISS DERIVATIVES
2971 0 : CMUDENR=0.d0
2972 0 : CMUDT=0.d0
2973 0 : CMUDTT=0.d0
2974 0 : return
2975 : end if
2976 : ! CALCULATE DERIVATIVES:
2977 : ! 1: derivatives of CHI over THETA and T
2978 : ! (a): Non-relativistic result:
2979 0 : THETA52=THETA32*THETA
2980 0 : CHIDY=-XDF/THETA52 ! d\chi/d\theta
2981 0 : CHIDYY=(XDFF/pow4(THETA)-2.5d0*CHIDY)/THETA ! d^2\chi/d\theta^2
2982 : ! (b): Relativistic corrections:
2983 0 : if (THETA>1.d-5) then
2984 0 : Q1D=-Q1/(1.d0-T1)
2985 0 : Q1DD=-Q1D*(1.d0+T1)/(1.d0-T1)
2986 : else
2987 0 : Q1D=-1.5d0/pow2(THETA)
2988 0 : Q1DD=-2.d0*Q1D/THETA
2989 : end if
2990 0 : Q2D=-12.d0/THETA52 ! d q_2 / d \theta
2991 0 : Q2DD=30.d0/(THETA52*THETA) ! d^2 q_2 / d \theta^2
2992 0 : U3D=-2.d0*T1*T1
2993 0 : D3D=PARF*PARG*THETAG/THETA+PARB*T1*T1*THETAC*(PARC/THETA-2.d0)
2994 : D3DD=PARF*PARG*(PARG-1.d0)*THETAG/pow2(THETA) &
2995 0 : + PARB*T1*T1*THETAC*(PARC*(PARC-1.d0)/pow2(THETA)-4.d0*PARC/THETA+4.d0)
2996 0 : Q3D=(D3D*U3/D3-U3D)/D3
2997 0 : Q3DD=(2.d0*U3D+(2.d0*U3D*D3D+U3*D3DD)/D3-2.d0*U3*pow2(D3D/D3))/D3
2998 0 : GDY=TEMR*(Q1D*SQT+(Q2D*Q3+Q2*Q3D)*TEMR) ! dG/d\theta
2999 0 : GDT=1.d0+1.5d0*Q1*SQT+2.d0*Q2*Q3*TEMR
3000 0 : GDYY=TEMR*(Q1DD*SQT+(Q2DD*Q3+2.d0*Q2D*Q3D+Q2*Q3DD)*TEMR)
3001 0 : GDTT=.75d0*Q1/SQT+2.d0*Q2*Q3
3002 0 : GDYT=1.5d0*Q1D*SQT+2.d0*(Q2D*Q3+Q2*Q3D)*TEMR
3003 0 : HDY=(-.5d0/pow2(THETA)+Q2D+.5d0*(Q2D-Q2/THETA)/THETA*TEMR)*TEMR
3004 0 : HDT=(.5d0+Q2*TEMR)/THETA+Q2
3005 : HDYY=TEMR/pow3(THETA)+Q2DD*TEMR &
3006 0 : + TEMR*TEMR*(.5d0*Q2DD-Q2D/THETA+Q2/pow2(THETA))/THETA
3007 0 : HDTT=Q2/THETA
3008 0 : HDYT=Q2D*(1.d0+TEMR/THETA)-(.5d0+Q2*TEMR)/pow2(THETA)
3009 0 : CTY=GDY/G-HDY/H
3010 0 : CTT=GDT/G-HDT/H
3011 0 : GH=G/H
3012 0 : CTDY=GH*CTY
3013 0 : CTDT=GH*CTT
3014 0 : CTDYY=CTDY*CTY+GH*(GDYY/G-pow2(GDY/G)-HDYY/H+pow2(HDY/H))
3015 0 : CTDTT=CTDT*CTT+GH*(GDTT/G-pow2(GDT/G)-HDTT/H+pow2(HDT/H))
3016 0 : CTDYT=CTDT*CTY+GH*(GDYT/G-GDY*GDT/pow2(G)-HDYT/H+HDY*HDT/pow2(H))
3017 0 : CHIDY=CHIDY-1.5d0*CTDY/CT
3018 0 : CHIDT=-1.5d0*CTDT/CT
3019 0 : CHIDYY=CHIDYY+1.5d0*(pow2(CTDY/CT)-CTDYY/CT)
3020 0 : CHIDTT=1.5d0*(pow2(CTDT/CT)-CTDTT/CT)
3021 0 : CHIDYT=1.5d0*(CTDY*CTDT/pow2(CT)-CTDYT/CT)
3022 0 : CMUDENR=-pow2(THETA*PF0)/(3.d0*DENR*(1.d0+TF))*CHIDY
3023 0 : CMUDT=CHI+THETA*CHIDY+TEMR*CHIDT
3024 0 : CMUDTT=2.d0*(CHIDY/TF+CHIDT+THETA*CHIDYT)+THETA/TF*CHIDYY+TEMR*CHIDTT
3025 0 : return
3026 : end subroutine CHEMFIT7
3027 :
3028 : end module pc_eos
|
