Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010 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 skye_ideal
21 :
22 : use math_lib
23 : use auto_diff
24 : use const_def, only: dp, pi, amu, planck_h, avo, crad, kerg
25 :
26 : implicit none
27 :
28 : private
29 : public :: compute_F_rad
30 : public :: compute_F_ideal_ion
31 : public :: compute_xne
32 : public :: compute_ideal_ele
33 :
34 : real(dp), parameter :: sifac = planck_h * planck_h * planck_h / (2d0 * pi * amu * sqrt(2d0 * pi * amu))
35 :
36 : logical, parameter :: dbg = .false.
37 :
38 : contains
39 :
40 91040 : type(auto_diff_real_2var_order3) function compute_F_rad(temp, den) result(Frad)
41 : type(auto_diff_real_2var_order3), intent(in) :: temp, den
42 :
43 : ! F = - P / rho
44 45520 : Frad = -crad * pow4(temp) / (3d0 * den)
45 :
46 45520 : end function compute_F_rad
47 :
48 45520 : type(auto_diff_real_2var_order3) function compute_F_ideal_ion(temp, den, abar, species, weights, ya) result(F_ideal_ion)
49 : type(auto_diff_real_2var_order3), intent(in) :: temp, den
50 : integer, intent(in) :: species
51 : real(dp), intent(in) :: weights(species), ya(species), abar
52 :
53 : integer :: j
54 : type(auto_diff_real_2var_order3) :: n, nj, nQ, nQj
55 : type(auto_diff_real_2var_order3) :: y, yy, z
56 : type(auto_diff_real_2var_order3) :: kt, xni, etaion
57 :
58 : ! Helper quantity
59 45520 : kt = kerg * temp
60 45520 : n = den / (amu * abar) ! Number density of ions
61 45520 : nQ = pow(kT, 1.5d0) / sifac ! Quantum density for a 1-amu species
62 :
63 45520 : F_ideal_ion = 0d0
64 356434 : do j=1,species
65 310914 : nj = ya(j) * n
66 310914 : nQj = nQ * pow(weights(j), 1.5d0)
67 356434 : F_ideal_ion = F_ideal_ion + ya(j) * (log(nj / nQj) - 1d0)
68 : end do
69 45520 : F_ideal_ion = F_ideal_ion * kt / (amu * abar)
70 :
71 45520 : end function compute_F_ideal_ion
72 :
73 45520 : type(auto_diff_real_2var_order3) function compute_xne(den, ytot1, zbar) result(xne)
74 : type(auto_diff_real_2var_order3), intent(in) :: den
75 : real(dp), intent(in) :: ytot1, zbar
76 : type(auto_diff_real_2var_order3) :: xni
77 :
78 : ! xne is the electron density due to matter, and so is proportional to density and independent of temperature
79 : ! for a fully ionized system.
80 45520 : xni = avo * ytot1 * den
81 45520 : xne = xni * zbar
82 :
83 45520 : end function compute_xne
84 :
85 45520 : subroutine compute_ideal_ele(temp_in, den_in, din_in, logtemp_in, logden_in, &
86 : zbar, ytot1, ye, ht, F, adr_etaele, adr_xnefer, ierr)
87 : use helm_polynomials
88 : use eos_def
89 : real(dp), intent(in) :: temp_in, den_in, din_in, logtemp_in, logden_in, zbar, ytot1, ye
90 : type (Helm_Table), pointer, intent(in) :: ht
91 :
92 : ! Intermediates
93 45520 : real(dp) :: x, y, z, ww, zz, xni
94 91040 : real(dp) :: temp, den, din, logtemp, logden
95 : integer :: k
96 :
97 : !..for the interpolations
98 : integer :: iat, jat
99 : real(dp) :: dth, dt2, dti, dt2i, dt3i, dd, &
100 1684240 : xt, xd, mxt, mxd, fi(36), &
101 45520 : dindd, dinda, dindz, dindda, dinddz, dindaa, &
102 45520 : dindaz, dindzz, dinddaa, dinddaz, &
103 : w0t, w1t, w2t, w0mt, w1mt, w2mt, &
104 : w0d, w1d, w2d, w0md, w1md, w2md, &
105 : dpepdd_in, dpepddd_in, dpepddt_in
106 :
107 45520 : real(dp) :: si0t, si1t, si2t, si0mt, si1mt, si2mt, &
108 45520 : si0d, si1d, si2d, si0md, si1md, si2md, &
109 45520 : dsi0t, dsi1t, dsi2t, dsi0mt, dsi1mt, dsi2mt, &
110 45520 : dsi0d, dsi1d, dsi2d, dsi0md, dsi1md, dsi2md, &
111 45520 : ddsi0t, ddsi1t, ddsi2t, ddsi0mt, ddsi1mt, ddsi2mt, &
112 45520 : ddsi0d, ddsi1d, ddsi2d, ddsi0md, ddsi1md, ddsi2md, &
113 45520 : dddsi0t, dddsi1t, dddsi2t, &
114 45520 : dddsi0mt, dddsi1mt, dddsi2mt, &
115 45520 : dddsi0d, dddsi1d, dddsi2d, &
116 45520 : dddsi0md, dddsi1md, dddsi2md
117 :
118 : ! For some results we don't need, but which helm_electron_positron.dek computes
119 : integer :: elemult
120 :
121 45520 : real(dp) :: detada,detadz
122 45520 : real(dp) :: detadda,detaddz
123 45520 : real(dp) :: detadta,detadtz,detadaa,detadaz,detadzz
124 45520 : real(dp) :: detadddd, detadddt, detaddtt, detadttt
125 :
126 : real(dp) :: dpepda,dpepdz
127 : real(dp) :: dpepdda,dpepddz
128 : real(dp) :: dpepdta,dpepdtz,dpepdaa
129 : real(dp) :: dpepdaz,dpepdzz
130 : real(dp) :: deepda,deepdz
131 : real(dp) :: deepdda,deepddz
132 : real(dp) :: deepdta,deepdtz,deepdaa
133 : real(dp) :: deepdaz,deepdzz
134 : real(dp) :: dsepda,dsepdz
135 : real(dp) :: dsepdda,dsepddz
136 : real(dp) :: dsepdta,dsepdtz,dsepdaa
137 : real(dp) :: dsepdaz,dsepdzz
138 :
139 : real(dp) :: etapos,zeff
140 :
141 : real(dp) :: dpeledd,dpeledt,dpeleda,dpeledz
142 : real(dp) :: dpeleddd,dpeleddt,dpeledda,dpeleddz
143 : real(dp) :: dpeledtt,dpeledta,dpeledtz,dpeledaa
144 : real(dp) :: dpeledaz,dpeledzz
145 : real(dp) :: deeledd,deeledt,deeleda,deeledz
146 : real(dp) :: deeleddd,deeleddt,deeledda,deeleddz
147 : real(dp) :: deeledtt,deeledta,deeledtz,deeledaa
148 : real(dp) :: deeledaz,deeledzz
149 : real(dp) :: dseledd,dseledt,dseleda,dseledz
150 : real(dp) :: dseleddd,dseleddt,dseledda,dseleddz
151 : real(dp) :: dseledtt,dseledta,dseledtz,dseledaa
152 : real(dp) :: dseledaz,dseledzz
153 :
154 : real(dp) :: ppos,dpposdd,dpposdt,dpposda,dpposdz
155 : real(dp) :: dpposddd,dpposddt,dpposdda,dpposddz
156 : real(dp) :: dpposdtt,dpposdta,dpposdtz,dpposdaa
157 : real(dp) :: dpposdaz,dpposdzz
158 : real(dp) :: epos,deposdd,deposdt,deposda,deposdz
159 : real(dp) :: deposddd,deposddt,deposdda,deposddz
160 : real(dp) :: deposdtt,deposdta,deposdtz,deposdaa
161 : real(dp) :: deposdaz,deposdzz
162 : real(dp) :: spos,dsposdd,dsposdt,dsposda,dsposdz
163 : real(dp) :: dsposddd,dsposddt,dsposdda,dsposddz
164 : real(dp) :: dsposdtt,dsposdta,dsposdtz,dsposdaa
165 : real(dp) :: dsposdaz,dsposdzz
166 :
167 : real(dp) :: xne
168 45520 : real(dp) :: dxnedd,dxnedt,dxneda,dxnedz
169 45520 : real(dp) :: dxneddd,dxneddt,dxnedda,dxneddz
170 45520 : real(dp) :: dxnedtt,dxnedta,dxnedtz,dxnedaa
171 45520 : real(dp) :: dxnedaz,dxnedzz
172 :
173 45520 : real(dp) :: xnefer,dxneferdd,dxneferdt,dxneferda,dxneferdz
174 : real(dp) :: dxneferddd,dxneferddt,dxneferdda,dxneferddz
175 : real(dp) :: dxneferdtt,dxneferdta,dxneferdtz,dxneferdaa
176 : real(dp) :: dxneferdaz,dxneferdzz
177 : real(dp) :: xnpfer,dxnpferdd,dxnpferdt,dxnpferda,dxnpferdz
178 : real(dp) :: dxnpferddd,dxnpferddt,dxnpferdda,dxnpferddz
179 : real(dp) :: dxnpferdtt,dxnpferdta,dxnpferdtz,dxnpferdaa
180 : real(dp) :: dxnpferdaz,dxnpferdzz
181 45520 : real(dp) :: dxneferdddd, dxneferdddt, dxneferddtt, dxneferdttt
182 45520 : real(dp) :: dxneferddda, dxneferdtta, dxneferddta
183 45520 : real(dp) :: dxneferdddz, dxneferdttz, dxneferddtz
184 45520 : real(dp) :: detaddda, detadtta, detaddta
185 45520 : real(dp) :: detadddz, detadttz, detaddtz
186 :
187 : ! Results we do need
188 45520 : real(dp) :: free, df_d, df_t, df_dd, df_tt, df_dt, &
189 45520 : df_ttt, df_dtt, df_ddt, df_ddd
190 45520 : real(dp) :: etaele,detadd,detadt
191 45520 : real(dp) :: detaddd,detaddt,detadtt
192 : type(auto_diff_real_2var_order3), intent(out) :: adr_etaele, adr_xnefer, F
193 :
194 : integer, intent(out) :: ierr
195 :
196 :
197 45520 : temp = temp_in
198 45520 : logtemp = logtemp_in
199 45520 : den = den_in
200 45520 : logden = logden_in
201 45520 : din = din_in
202 :
203 :
204 : !..density derivatives
205 45520 : dindd = ye
206 45520 : dinda = -din*ytot1
207 45520 : dindz = den*ytot1
208 45520 : dindda = -ye*ytot1
209 45520 : dinddz = ytot1
210 45520 : ww = ytot1*ytot1
211 45520 : dindaa = 2.0d0*din*ww
212 45520 : dindaz = -den*ww
213 45520 : dinddaa = 2.0d0*ye*ww
214 45520 : dinddaz = -ww
215 45520 : dindzz = 0.0d0
216 :
217 :
218 : !..hash locate this temperature and density
219 45520 : jat = int((logtemp - ht% logtlo)*ht% logtstpi) + 1
220 45520 : jat = max(1,min(jat,ht% jmax-1))
221 45520 : iat = int((log10(din) - ht% logdlo)*ht% logdstpi) + 1
222 45520 : iat = max(1,min(iat,ht% imax-1))
223 :
224 :
225 : !..access the table locations only once
226 45520 : fi(1) = ht% f(iat,jat)
227 45520 : fi(2) = ht% f(iat+1,jat)
228 45520 : fi(3) = ht% f(iat,jat+1)
229 45520 : fi(4) = ht% f(iat+1,jat+1)
230 45520 : fi(5) = ht% ft(iat,jat)
231 45520 : fi(6) = ht% ft(iat+1,jat)
232 45520 : fi(7) = ht% ft(iat,jat+1)
233 45520 : fi(8) = ht% ft(iat+1,jat+1)
234 45520 : fi(9) = ht% ftt(iat,jat)
235 45520 : fi(10) = ht% ftt(iat+1,jat)
236 45520 : fi(11) = ht% ftt(iat,jat+1)
237 45520 : fi(12) = ht% ftt(iat+1,jat+1)
238 45520 : fi(13) = ht% fd(iat,jat)
239 45520 : fi(14) = ht% fd(iat+1,jat)
240 45520 : fi(15) = ht% fd(iat,jat+1)
241 45520 : fi(16) = ht% fd(iat+1,jat+1)
242 45520 : fi(17) = ht% fdd(iat,jat)
243 45520 : fi(18) = ht% fdd(iat+1,jat)
244 45520 : fi(19) = ht% fdd(iat,jat+1)
245 45520 : fi(20) = ht% fdd(iat+1,jat+1)
246 45520 : fi(21) = ht% fdt(iat,jat)
247 45520 : fi(22) = ht% fdt(iat+1,jat)
248 45520 : fi(23) = ht% fdt(iat,jat+1)
249 45520 : fi(24) = ht% fdt(iat+1,jat+1)
250 45520 : fi(25) = ht% fddt(iat,jat)
251 45520 : fi(26) = ht% fddt(iat+1,jat)
252 45520 : fi(27) = ht% fddt(iat,jat+1)
253 45520 : fi(28) = ht% fddt(iat+1,jat+1)
254 45520 : fi(29) = ht% fdtt(iat,jat)
255 45520 : fi(30) = ht% fdtt(iat+1,jat)
256 45520 : fi(31) = ht% fdtt(iat,jat+1)
257 45520 : fi(32) = ht% fdtt(iat+1,jat+1)
258 45520 : fi(33) = ht% fddtt(iat,jat)
259 45520 : fi(34) = ht% fddtt(iat+1,jat)
260 45520 : fi(35) = ht% fddtt(iat,jat+1)
261 45520 : fi(36) = ht% fddtt(iat+1,jat+1)
262 :
263 : !..various differences
264 45520 : xt = max( (temp - ht% t(jat)) * ht% dti_sav(jat), 0.0d0)
265 45520 : xd = max( (din - ht% d(iat)) * ht% ddi_sav(iat), 0.0d0)
266 45520 : mxt = 1.0d0 - xt
267 45520 : mxd = 1.0d0 - xd
268 :
269 : !..the six density and six temperature basis functions
270 45520 : si0t = psi0(xt)
271 45520 : si1t = psi1(xt) * ht% dt_sav(jat)
272 45520 : si2t = psi2(xt) * ht% dt2_sav(jat)
273 :
274 45520 : si0mt = psi0(mxt)
275 45520 : si1mt = -psi1(mxt) * ht% dt_sav(jat)
276 45520 : si2mt = psi2(mxt) * ht% dt2_sav(jat)
277 :
278 45520 : si0d = psi0(xd)
279 45520 : si1d = psi1(xd) * ht% dd_sav(iat)
280 45520 : si2d = psi2(xd) * ht% dd2_sav(iat)
281 :
282 45520 : si0md = psi0(mxd)
283 45520 : si1md = -psi1(mxd) * ht% dd_sav(iat)
284 45520 : si2md = psi2(mxd) * ht% dd2_sav(iat)
285 :
286 : !..first derivatives of the weight functions
287 45520 : dsi0t = dpsi0(xt) * ht% dti_sav(jat)
288 45520 : dsi1t = dpsi1(xt)
289 45520 : dsi2t = dpsi2(xt) * ht% dt_sav(jat)
290 :
291 45520 : dsi0mt = -dpsi0(mxt) * ht% dti_sav(jat)
292 45520 : dsi1mt = dpsi1(mxt)
293 45520 : dsi2mt = -dpsi2(mxt) * ht% dt_sav(jat)
294 :
295 45520 : dsi0d = dpsi0(xd) * ht% ddi_sav(iat)
296 45520 : dsi1d = dpsi1(xd)
297 45520 : dsi2d = dpsi2(xd) * ht% dd_sav(iat)
298 :
299 45520 : dsi0md = -dpsi0(mxd) * ht% ddi_sav(iat)
300 45520 : dsi1md = dpsi1(mxd)
301 45520 : dsi2md = -dpsi2(mxd) * ht% dd_sav(iat)
302 :
303 : !..second derivatives of the weight functions
304 45520 : ddsi0t = ddpsi0(xt) * ht% dt2i_sav(jat)
305 45520 : ddsi1t = ddpsi1(xt) * ht% dti_sav(jat)
306 45520 : ddsi2t = ddpsi2(xt)
307 :
308 45520 : ddsi0mt = ddpsi0(mxt) * ht% dt2i_sav(jat)
309 45520 : ddsi1mt = -ddpsi1(mxt) * ht% dti_sav(jat)
310 45520 : ddsi2mt = ddpsi2(mxt)
311 :
312 45520 : ddsi0d = ddpsi0(xd) * ht% dd2i_sav(iat)
313 45520 : ddsi1d = ddpsi1(xd) * ht% ddi_sav(iat)
314 45520 : ddsi2d = ddpsi2(xd)
315 :
316 45520 : ddsi0md = ddpsi0(mxd) * ht% dd2i_sav(iat)
317 45520 : ddsi1md = -ddpsi1(mxd) * ht% ddi_sav(iat)
318 45520 : ddsi2md = ddpsi2(mxd)
319 :
320 : !..third derivatives of the weight functions
321 45520 : dddsi0t = dddpsi0(xt) * ht% dt3i_sav(jat)
322 45520 : dddsi1t = dddpsi1(xt) * ht% dt2i_sav(jat)
323 45520 : dddsi2t = dddpsi2(xt) * ht% dti_sav(jat)
324 :
325 45520 : dddsi0mt = -dddpsi0(mxt) * ht% dt3i_sav(jat)
326 45520 : dddsi1mt = dddpsi1(mxt) * ht% dt2i_sav(jat)
327 45520 : dddsi2mt = -dddpsi2(mxt) * ht% dti_sav(jat)
328 :
329 45520 : dddsi0d = dddpsi0(xd) * ht% dd3i_sav(iat)
330 45520 : dddsi1d = dddpsi1(xd) * ht% dd2i_sav(iat)
331 45520 : dddsi2d = dddpsi2(xd) * ht% ddi_sav(iat)
332 :
333 45520 : dddsi0md = -dddpsi0(mxd) * ht% dd3i_sav(iat)
334 45520 : dddsi1md = dddpsi1(mxd) * ht% dd2i_sav(iat)
335 45520 : dddsi2md = -dddpsi2(mxd) * ht% ddi_sav(iat)
336 :
337 :
338 : !..the free energy
339 : free = h5(iat,jat,fi, &
340 : si0t, si1t, si2t, si0mt, si1mt, si2mt, &
341 45520 : si0d, si1d, si2d, si0md, si1md, si2md)
342 :
343 : !..first derivative with respect to density
344 : df_d = h5(iat,jat,fi, &
345 : si0t, si1t, si2t, si0mt, si1mt, si2mt, &
346 45520 : dsi0d, dsi1d, dsi2d, dsi0md, dsi1md, dsi2md)
347 :
348 : !..first derivative with respect to temperature
349 : df_t = h5(iat,jat,fi, &
350 : dsi0t, dsi1t, dsi2t, dsi0mt, dsi1mt, dsi2mt, &
351 45520 : si0d, si1d, si2d, si0md, si1md, si2md)
352 :
353 : !..second derivative with respect to density**2
354 : df_dd = h5(iat,jat,fi, &
355 : si0t, si1t, si2t, si0mt, si1mt, si2mt, &
356 45520 : ddsi0d, ddsi1d, ddsi2d, ddsi0md, ddsi1md, ddsi2md)
357 :
358 : !..second derivative with respect to temperature**2
359 : df_tt = h5(iat,jat,fi, &
360 : ddsi0t, ddsi1t, ddsi2t, ddsi0mt, ddsi1mt, ddsi2mt, &
361 45520 : si0d, si1d, si2d, si0md, si1md, si2md)
362 :
363 : !..second derivative with respect to temperature and density
364 : df_dt = h5(iat,jat,fi, &
365 : dsi0t, dsi1t, dsi2t, dsi0mt, dsi1mt, dsi2mt, &
366 45520 : dsi0d, dsi1d, dsi2d, dsi0md, dsi1md, dsi2md)
367 :
368 : !..third derivative with respect to temperature**3
369 : df_ttt = h5(iat,jat,fi, &
370 : dddsi0t, dddsi1t, dddsi2t, dddsi0mt, dddsi1mt, dddsi2mt, &
371 45520 : si0d, si1d, si2d, si0md, si1md, si2md)
372 :
373 : !..third derivative with respect to density temperature**2
374 : df_dtt = h5(iat,jat,fi, &
375 : ddsi0t, ddsi1t, ddsi2t, ddsi0mt, ddsi1mt, ddsi2mt, &
376 45520 : dsi0d, dsi1d, dsi2d, dsi0md, dsi1md, dsi2md)
377 :
378 : !..third derivative with respect to density**2 temperature
379 : df_ddt = h5(iat,jat,fi, &
380 : dsi0t, dsi1t, dsi2t, dsi0mt, dsi1mt, dsi2mt, &
381 45520 : ddsi0d, ddsi1d, ddsi2d, ddsi0md, ddsi1md, ddsi2md)
382 :
383 : !..third derivative with respect to density**3
384 : df_ddd = h5(iat,jat,fi, &
385 : si0t, si1t, si2t, si0mt, si1mt, si2mt, &
386 45520 : dddsi0d, dddsi1d, dddsi2d, dddsi0md, dddsi1md, dddsi2md)
387 :
388 : !..now get the pressure derivative with density, chemical potential, and
389 : !..electron positron number densities
390 : !..get the cubic interpolation weight functions
391 :
392 45520 : si0t = xpsi0(xt)
393 45520 : si1t = xpsi1(xt) * ht% dt_sav(jat)
394 45520 : si0mt = xpsi0(mxt)
395 45520 : si1mt = -xpsi1(mxt) * ht% dt_sav(jat)
396 :
397 45520 : si0d = xpsi0(xd)
398 45520 : si1d = xpsi1(xd) * ht% dd_sav(iat)
399 45520 : si0md = xpsi0(mxd)
400 45520 : si1md = -xpsi1(mxd) * ht% dd_sav(iat)
401 :
402 : !..first derivatives of weight functions
403 45520 : dsi0t = xdpsi0(xt) * ht% dti_sav(jat)
404 45520 : dsi1t = xdpsi1(xt)
405 45520 : dsi0mt = -xdpsi0(mxt) * ht% dti_sav(jat)
406 45520 : dsi1mt = xdpsi1(mxt)
407 :
408 45520 : dsi0d = xdpsi0(xd) * ht% ddi_sav(iat)
409 45520 : dsi1d = xdpsi1(xd)
410 45520 : dsi0md = -xdpsi0(mxd) * ht% ddi_sav(iat)
411 45520 : dsi1md = xdpsi1(mxd)
412 :
413 : !..second derivatives of weight functions
414 45520 : ddsi0t = xddpsi0(xt) * ht% dt2i_sav(jat)
415 45520 : ddsi1t = xddpsi1(xt) * ht% dti_sav(jat)
416 45520 : ddsi0mt = xddpsi0(mxt) * ht% dt2i_sav(jat)
417 45520 : ddsi1mt = -xddpsi1(mxt) * ht% dti_sav(jat)
418 :
419 45520 : ddsi0d = xddpsi0(xd) * ht% dd2i_sav(iat)
420 45520 : ddsi1d = xddpsi1(xd) * ht% ddi_sav(iat)
421 45520 : ddsi0md = xddpsi0(mxd) * ht% dd2i_sav(iat)
422 45520 : ddsi1md = -xddpsi1(mxd) * ht% ddi_sav(iat)
423 :
424 : !..third derivatives of weight functions
425 45520 : dddsi0t = xdddpsi0(xt) * ht% dt3i_sav(jat)
426 45520 : dddsi1t = xdddpsi1(xt) * ht% dt2i_sav(jat)
427 45520 : dddsi0mt = -xdddpsi0(mxt) * ht% dt3i_sav(jat)
428 45520 : dddsi1mt = xdddpsi1(mxt) * ht% dt2i_sav(jat)
429 :
430 45520 : dddsi0d = xdddpsi0(xd) * ht% dd3i_sav(iat)
431 45520 : dddsi1d = xdddpsi1(xd) * ht% dd2i_sav(iat)
432 45520 : dddsi0md = -xdddpsi0(mxd) * ht% dd3i_sav(iat)
433 45520 : dddsi1md = xdddpsi1(mxd) * ht% dd2i_sav(iat)
434 :
435 : !..look in the electron chemical potential table only once
436 45520 : fi(1) = ht% ef(iat,jat)
437 45520 : fi(2) = ht% ef(iat+1,jat)
438 45520 : fi(3) = ht% ef(iat,jat+1)
439 45520 : fi(4) = ht% ef(iat+1,jat+1)
440 45520 : fi(5) = ht% eft(iat,jat)
441 45520 : fi(6) = ht% eft(iat+1,jat)
442 45520 : fi(7) = ht% eft(iat,jat+1)
443 45520 : fi(8) = ht% eft(iat+1,jat+1)
444 45520 : fi(9) = ht% efd(iat,jat)
445 45520 : fi(10) = ht% efd(iat+1,jat)
446 45520 : fi(11) = ht% efd(iat,jat+1)
447 45520 : fi(12) = ht% efd(iat+1,jat+1)
448 45520 : fi(13) = ht% efdt(iat,jat)
449 45520 : fi(14) = ht% efdt(iat+1,jat)
450 45520 : fi(15) = ht% efdt(iat,jat+1)
451 45520 : fi(16) = ht% efdt(iat+1,jat+1)
452 :
453 :
454 : !..electron chemical potential etaele
455 : etaele = h3(iat,jat,fi, &
456 : si0t, si1t, si0mt, si1mt, &
457 45520 : si0d, si1d, si0md, si1md)
458 :
459 : !..first derivatives
460 : x = h3(iat,jat,fi, &
461 : si0t, si1t, si0mt, si1mt, &
462 45520 : dsi0d, dsi1d, dsi0md, dsi1md)
463 45520 : detadd = ye * x
464 : detadt = h3(iat,jat,fi, &
465 : dsi0t, dsi1t, dsi0mt, dsi1mt, &
466 45520 : si0d, si1d, si0md, si1md)
467 45520 : detada = -x * din * ytot1
468 45520 : detadz = x * den * ytot1
469 :
470 : !..second derivatives
471 : y = h3(iat,jat,fi, &
472 : si0t, si1t, si0mt, si1mt, &
473 45520 : ddsi0d, ddsi1d, ddsi0md, ddsi1md)
474 45520 : detaddd = ye * ye * y
475 : z = h3(iat,jat,fi, &
476 : dsi0t, dsi1t, dsi0mt, dsi1mt, &
477 45520 : dsi0d, dsi1d, dsi0md, dsi1md)
478 45520 : detaddt = ye * z
479 45520 : detadda = -ye*ytot1*x + ye*y*dinda
480 45520 : detaddz = ytot1*x + ye*y*dindz
481 : detadtt = h3(iat,jat,fi, &
482 : ddsi0t, ddsi1t, ddsi0mt, ddsi1mt, &
483 45520 : si0d, si1d, si0md, si1md)
484 45520 : detadta = z * dinda
485 45520 : detadtz = z * dindz
486 45520 : detadaa = (y*din + 2.0d0*x)*din*ww
487 45520 : detadaz = -(y*dindz*din + x*den*ytot1)*ytot1
488 45520 : detadzz = y*dindz*den*ytot1
489 :
490 : !..third derivatives
491 : y = h3(iat,jat,fi, &
492 : si0t, si1t, si0mt, si1mt, &
493 45520 : dddsi0d, dddsi1d, dddsi0md, dddsi1md)
494 45520 : detadddd = ye * ye * ye * y ! Actual interpolation variable is ye * rho, so we multiply by ye to get d/d(density)
495 : ! ! d/drho^3
496 :
497 : detadttt = h3(iat,jat,fi, &
498 : dddsi0t, dddsi1t, dddsi0mt, dddsi1mt, &
499 45520 : si0d, si1d, si0md, si1md) ! d/dT^3
500 :
501 : y = h3(iat,jat,fi, &
502 : dsi0t, dsi1t, dsi0mt, dsi1mt, &
503 45520 : ddsi0d, ddsi1d, ddsi0md, ddsi1md)
504 45520 : detadddt = ye * ye * y ! d/drho^2 d/dT
505 :
506 : y = h3(iat,jat,fi, &
507 : ddsi0t, ddsi1t, ddsi0mt, ddsi1mt, &
508 45520 : dsi0d, dsi1d, dsi0md, dsi1md)
509 :
510 45520 : detaddtt = ye * y ! d/drho d/dT^2
511 :
512 : ! dg/da = dg/d(din) d(din)/da = dg/d(din) d(ye*rho)/da
513 : ! = dg/d(din) d(z rho / a)/da = -(z rho / a^2) dg/d(din)
514 : ! = -(z rho / a^2) dg/d(Rho) d(Rho)/d(din)
515 : ! = -(z rho / a^2) dg/d(Rho) (1 / ye)
516 : ! = -(rho / a) dg/d(Rho) = - ytot1 * rho * dg/d(Rho)
517 : ! = -ytot1 * den * dg/d(Rho)
518 45520 : detaddda = -den * ytot1 * detadddd
519 45520 : detadtta = -den * ytot1 * detaddtt
520 45520 : detaddta = -den * ytot1 * detadddt
521 :
522 : ! dg/dz = dg/d(din) d(din)/dz = dg/d(din) d(ye*rho)/dz
523 : ! = dg/d(din) d(z rho / a)/dz = (rho / a) dg/d(din)
524 : ! = (rho / a) dg/d(Rho) d(Rho)/d(din)
525 : ! = (rho / a) dg/d(Rho) (1 / ye)
526 : ! = (rho / z) dg/d(Rho)
527 : ! = ytot1 * din * dg/d(Rho)
528 45520 : detadddz = din * ytot1 * detadddd
529 45520 : detadttz = din * ytot1 * detaddtt
530 45520 : detaddtz = din * ytot1 * detadddt
531 :
532 : !..look in the number density table only once
533 45520 : fi(1) = ht% xf(iat,jat)
534 45520 : fi(2) = ht% xf(iat+1,jat)
535 45520 : fi(3) = ht% xf(iat,jat+1)
536 45520 : fi(4) = ht% xf(iat+1,jat+1)
537 45520 : fi(5) = ht% xft(iat,jat)
538 45520 : fi(6) = ht% xft(iat+1,jat)
539 45520 : fi(7) = ht% xft(iat,jat+1)
540 45520 : fi(8) = ht% xft(iat+1,jat+1)
541 45520 : fi(9) = ht% xfd(iat,jat)
542 45520 : fi(10) = ht% xfd(iat+1,jat)
543 45520 : fi(11) = ht% xfd(iat,jat+1)
544 45520 : fi(12) = ht% xfd(iat+1,jat+1)
545 45520 : fi(13) = ht% xfdt(iat,jat)
546 45520 : fi(14) = ht% xfdt(iat+1,jat)
547 45520 : fi(15) = ht% xfdt(iat,jat+1)
548 45520 : fi(16) = ht% xfdt(iat+1,jat+1)
549 :
550 : !..electron + positron number densities
551 : xnefer = h3(iat,jat,fi, &
552 : si0t, si1t, si0mt, si1mt, &
553 45520 : si0d, si1d, si0md, si1md)
554 :
555 : !..first derivatives
556 : x = h3(iat,jat,fi, &
557 : si0t, si1t, si0mt, si1mt, &
558 45520 : dsi0d, dsi1d, dsi0md, dsi1md)
559 45520 : x = max(x,1.0d-30)
560 45520 : dxnedd = ye * x
561 : dxnedt = h3(iat,jat,fi, &
562 : dsi0t, dsi1t, dsi0mt, dsi1mt, &
563 45520 : si0d, si1d, si0md, si1md)
564 45520 : dxneda = -x * din * ytot1
565 45520 : dxnedz = x * den * ytot1
566 :
567 : !..second derivatives
568 : y = h3(iat,jat,fi, &
569 : si0t, si1t, si0mt, si1mt, &
570 45520 : ddsi0d, ddsi1d, ddsi0md, ddsi1md)
571 45520 : dxneddd = ye * ye * y
572 : z = h3(iat,jat,fi, &
573 : dsi0t, dsi1t, dsi0mt, dsi1mt, &
574 45520 : dsi0d, dsi1d, dsi0md, dsi1md)
575 45520 : dxneddt = ye * z
576 45520 : dxnedda = -ye*ytot1*x + ye*y*dinda
577 45520 : dxneddz = ytot1*x + ye*y*dindz
578 : dxnedtt = h3(iat,jat,fi, &
579 : ddsi0t, ddsi1t, ddsi0mt, ddsi1mt, &
580 45520 : si0d, si1d, si0md, si1md)
581 45520 : dxnedta = z * dinda
582 45520 : dxnedtz = z * dindz
583 45520 : dxnedaa = (y*din + 2.0d0*x)*din*ytot1*ytot1
584 45520 : dxnedaz = -(y*dindz*din + x*den*ytot1)*ytot1
585 45520 : dxnedzz = y*dindz*den*ytot1
586 :
587 : !..third derivatives
588 : y = h3(iat,jat,fi, &
589 : si0t, si1t, si0mt, si1mt, &
590 45520 : dddsi0d, dddsi1d, dddsi0md, dddsi1md)
591 45520 : dxneferdddd = ye * ye * ye * y ! Actual interpolation variable is ye * rho, so we multiply by ye to get d/d(density)
592 : ! ! d/drho^3
593 :
594 : dxneferdttt = h3(iat,jat,fi, &
595 : dddsi0t, dddsi1t, dddsi0mt, dddsi1mt, &
596 45520 : si0d, si1d, si0md, si1md) ! d/dT^3
597 :
598 :
599 : y = h3(iat,jat,fi, &
600 : dsi0t, dsi1t, dsi0mt, dsi1mt, &
601 45520 : ddsi0d, ddsi1d, ddsi0md, ddsi1md)
602 45520 : dxneferdddt = ye * ye * y ! d/drho^2 d/dT
603 :
604 : y = h3(iat,jat,fi, &
605 : ddsi0t, ddsi1t, ddsi0mt, ddsi1mt, &
606 45520 : dsi0d, dsi1d, dsi0md, dsi1md)
607 :
608 45520 : dxneferddtt = ye * y ! d/drho d/dT^2
609 :
610 :
611 : ! dg/da = dg/d(din) d(din)/da = dg/d(din) d(ye*rho)/da
612 : ! = dg/d(din) d(z rho / a)/da = -(z rho / a^2) dg/d(din)
613 : ! = -(z rho / a^2) dg/d(Rho) d(Rho)/d(din)
614 : ! = -(z rho / a^2) dg/d(Rho) (1 / ye)
615 : ! = -(rho / a) dg/d(Rho) = - ytot1 * rho * dg/d(Rho)
616 : ! = -ytot1 * den * dg/d(Rho)
617 45520 : dxneferddda = -den * ytot1 * dxneferdddd
618 45520 : dxneferdtta = -den * ytot1 * dxneferddtt
619 45520 : dxneferddta = -den * ytot1 * dxneferdddt
620 :
621 : ! dg/dz = dg/d(din) d(din)/dz = dg/d(din) d(ye*rho)/dz
622 : ! = dg/d(din) d(z rho / a)/dz = (rho / a) dg/d(din)
623 : ! = (rho / a) dg/d(Rho) d(Rho)/d(din)
624 : ! = (rho / a) dg/d(Rho) (1 / ye)
625 : ! = (rho / z) dg/d(Rho)
626 : ! = ytot1 * din * dg/d(Rho)
627 45520 : dxneferdddz = din * ytot1 * dxneferdddd
628 45520 : dxneferdttz = din * ytot1 * dxneferddtt
629 45520 : dxneferddtz = din * ytot1 * dxneferdddt
630 :
631 :
632 : ! Pack quantities into auto_diff types
633 :
634 : ! Free energy
635 45520 : F%val = free * ye
636 :
637 45520 : F%d1val1 = df_t * ye
638 45520 : F%d1val2 = df_d * pow2(ye)
639 :
640 45520 : F%d2val1 = df_tt * ye
641 45520 : F%d1val1_d1val2 = df_dt * pow2(ye)
642 45520 : F%d2val2 = df_dd * pow3(ye)
643 :
644 45520 : F%d3val1 = df_ttt * ye
645 45520 : F%d2val1_d1val2 = df_dtt * pow2(ye)
646 45520 : F%d1val1_d2val2 = df_ddt * pow3(ye)
647 45520 : F%d3val2 = df_ddd * pow4(ye)
648 :
649 :
650 : ! Electron chemical potential
651 45520 : adr_etaele = etaele
652 :
653 45520 : adr_etaele%d1val1 = detadt
654 45520 : adr_etaele%d1val2 = detadd
655 : ! adr_etaele%d1val3 = detada ! d/dabar
656 : ! adr_etaele%d1val4 = detadz ! d/dzbar
657 :
658 45520 : adr_etaele%d2val1 = detadtt
659 45520 : adr_etaele%d2val2 = detaddd
660 45520 : adr_etaele%d1val1_d1val2 = detaddt
661 : ! adr_etaele%d1val1_d1val3 = detadta ! d/dabar
662 : ! adr_etaele%d1val1_d1val4 = detadtz ! d/dzbar
663 : ! adr_etaele%d1val2_d1val3 = detadda ! d/dabar
664 : ! adr_etaele%d1val2_d1val4 = detaddz ! d/dzbar
665 :
666 45520 : adr_etaele%d3val1 = detadttt
667 45520 : adr_etaele%d2val1_d1val2 = detaddtt
668 45520 : adr_etaele%d1val1_d2val2 = detadddt
669 45520 : adr_etaele%d3val2 = detadddd
670 :
671 : ! adr_etaele%d2val1_d1val3 = detadtta
672 : ! adr_etaele%d2val2_d1val3 = detaddda
673 : ! adr_etaele%d1val1_d1val2_d1val3 = detaddta
674 : ! adr_etaele%d2val1_d1val4 = detadttz
675 : ! adr_etaele%d2val2_d1val4 = detadddz
676 : ! adr_etaele%d1val1_d1val2_d1val4 = detaddtz
677 :
678 : ! Electron density
679 45520 : adr_xnefer = xnefer
680 :
681 45520 : adr_xnefer%d1val1 = dxnedt
682 45520 : adr_xnefer%d1val2 = dxnedd
683 : ! adr_xnefer%d1val3 = dxneda ! d/dabar
684 : ! adr_xnefer%d1val4 = dxnedz ! d/dzbar
685 :
686 45520 : adr_xnefer%d2val1 = dxnedtt
687 45520 : adr_xnefer%d2val2 = dxneddd
688 45520 : adr_xnefer%d1val1_d1val2 = dxneddt
689 : ! adr_xnefer%d1val1_d1val3 = dxnedta ! d/dabar
690 : ! adr_xnefer%d1val1_d1val4 = dxnedtz ! d/dzbar
691 : ! adr_xnefer%d1val2_d1val3 = dxnedda ! d/dabar
692 : ! adr_xnefer%d1val2_d1val4 = dxneddz ! d/dzbar
693 :
694 45520 : adr_xnefer%d3val1 = dxneferdttt
695 45520 : adr_xnefer%d2val1_d1val2 = dxneferddtt
696 45520 : adr_xnefer%d1val1_d2val2 = dxneferdddt
697 45520 : adr_xnefer%d3val2 = dxneferdddd
698 :
699 : ! adr_xnefer%d2val1_d1val3 = dxneferdtta
700 : ! adr_xnefer%d2val2_d1val3 = dxneferddda
701 : ! adr_xnefer%d1val1_d1val2_d1val3 = dxneferddta
702 : ! adr_xnefer%d2val1_d1val4 = dxneferdttz
703 : ! adr_xnefer%d2val2_d1val4 = dxneferdddz
704 : ! adr_xnefer%d1val1_d1val2_d1val4 = dxneferddtz
705 :
706 45520 : end subroutine compute_ideal_ele
707 :
708 : end module skye_ideal
|
