Note
Go to the end to download the full example code.
Equations of state functions#
Machida06 EoS#
11 import matplotlib.pyplot as plt
12 import numpy as np
13
14 import shamrock
15
16 cs = 190.0
17 rho_c1 = 1.92e-13 * 1000 # g/cm^3 -> kg/m^3
18 rho_c2 = 3.84e-8 * 1000 # g/cm^3 -> kg/m^3
19 rho_c3 = 1.92e-3 * 1000 # g/cm^3 -> kg/m^3
20
21
22 si = shamrock.UnitSystem()
23 sicte = shamrock.Constants(si)
24 kb = sicte.kb()
25 print(kb)
26 mu = 2.375
27 mh = 1.00784 * sicte.dalton()
28 print(mu * mh * kb)
29
30 rho_plot = np.logspace(-15, 5, 1000)
31 P_plot = []
32 cs_plot = []
33 T_plot = []
34 for rho in rho_plot:
35 P, _cs, T = shamrock.phys.eos.eos_Machida06(
36 cs=cs, rho=rho, rho_c1=rho_c1, rho_c2=rho_c2, rho_c3=rho_c3, mu=mu, mh=mh, kb=kb
37 )
38 P_plot.append(P)
39 cs_plot.append(_cs)
40 T_plot.append(T)
41
42 plt.figure()
43 plt.plot(rho_plot, P_plot, label="P")
44 plt.yscale("log")
45 plt.xscale("log")
46 plt.xlabel("$\\rho$ [kg.m^-3]")
47 plt.ylabel("$P$ [Pa]")
48 plt.axvspan(rho_c1, rho_c2, color="grey", alpha=0.5)
49 plt.axvspan(rho_c3, rho_plot[-1], color="grey", alpha=0.5)
50
51 plt.figure()
52 plt.plot(rho_plot, cs_plot, label="cs")
53 plt.yscale("log")
54 plt.xlabel("$\\rho$ [kg.m^-3]")
55 plt.xscale("log")
56 plt.ylabel("$c_s$ [m/s]")
57 plt.axvspan(rho_c1, rho_c2, color="grey", alpha=0.5)
58 plt.axvspan(rho_c3, rho_plot[-1], color="grey", alpha=0.5)
59
60 plt.figure()
61 plt.plot(rho_plot, T_plot, label="T")
62 plt.yscale("log")
63 plt.xscale("log")
64 plt.xlabel("$\\rho$ [kg.m^-3]")
65 plt.ylabel("$T$ [K]")
66 plt.axvspan(rho_c1, rho_c2, color="grey", alpha=0.5)
67 plt.axvspan(rho_c3, rho_plot[-1], color="grey", alpha=0.5)
68
69
70 plt.tight_layout()
1.380649e-23
5.487664926077112e-50
Fermi gas EoS#
76 rho_plot = np.logspace(1, 20, 1000)
77 P_plot = []
78 cs_plot = []
79
80 for rho in rho_plot:
81 P, _cs = shamrock.phys.eos.eos_Fermi(mu_e=2, rho=rho)
82 P_plot.append(P)
83 cs_plot.append(_cs)
84
85 plt.figure()
86 plt.suptitle("Fermi Gas EoS")
87 plt.plot(rho_plot, P_plot, label="P", color="blue")
88 plt.yscale("log")
89 plt.xscale("log")
90 plt.xlabel("$\\rho$ [kg.m^-3]")
91 plt.ylabel("$P$ [Pa]", color="blue")
92 plt.legend()
93
94
95 ax = plt.twinx()
96 ax.plot(rho_plot, cs_plot, label="cs", color="orange")
97 ax.set_yscale("log")
98 ax.set_ylabel("$c_s$ [m/s]", color="orange")
99 ax.legend(loc="lower right")
<matplotlib.legend.Legend object at 0x7fb9bc47c560>
Total running time of the script: (0 minutes 1.071 seconds)
Estimated memory usage: 159 MB