Skip to content

Commit 59b8bd1

Browse files
riederrieder
committed
reworked example
1 parent 95ca82e commit 59b8bd1

1 file changed

Lines changed: 107 additions & 79 deletions

File tree

‎src/amuse/examples/plot_solar_comparison.py‎

Lines changed: 107 additions & 79 deletions
Original file line numberDiff line numberDiff line change
@@ -1,112 +1,140 @@
1-
import numpy
2-
from amuse.lab import *
3-
from amuse.units.optparse import OptionParser
4-
from matplotlib import pyplot
5-
from amuse.plot import scatter
6-
7-
from prepare_figure import single_frame
8-
from distinct_colours import get_distinct
9-
10-
from amuse.community.seba.interface import SeBa
11-
12-
def main(t_end, mass, z, Tstar, Lstar):
13-
14-
stellar_evolution_codes = [SeBa(), SSE(), MESA(), EVtwin()]
15-
label = ["SeBa", "SSE", "MESA", "EVtwin"]
1+
import argparse
2+
import numpy as np
3+
import matplotlib.pyplot as plt
4+
from amuse.units import units
5+
from amuse.datamodel import Particles
6+
from amuse.community.seba import Seba
7+
from amuse.community.sse import Sse
8+
from amuse.community.mesa import Mesa
9+
from amuse.community.evtwin import Evtwin
10+
11+
12+
def plot_solar_comparison(time_end, mass, z, Tstar, Lstar):
13+
stellar_evolution_codes = ["SeBa", "SSE", "MESA", "EVtwin"]
1614
marker = ["o", "v", "<", ">"]
1715

18-
pyplot.figure()
19-
pyplot.xlabel("$(T-T_\odot)/T_\odot)$")
20-
pyplot.ylabel("$(L-L_\odot)/L_\odot)$")
21-
pyplot.xlim(-0.006, 0.004)
22-
pyplot.ylim(-0.1, 0.1)
23-
pyplot.scatter([0], [0], marker="o", label="Sun", s=200, lw=0)
24-
25-
for si in range(len(stellar_evolution_codes)):
26-
stellar = stellar_evolution_codes[si]
16+
figure = plt.figure()
17+
ax = figure.add_subplot(1, 1, 1)
18+
plt.xlabel("$(T-T_\odot)/T_\odot)$")
19+
plt.ylabel("$(L-L_\odot)/L_\odot)$")
20+
plt.xlim(-0.006, 0.004)
21+
plt.ylim(-0.1, 0.1)
22+
plt.scatter([0], [0], marker="o", label="Sun", s=200, lw=0)
23+
24+
for si, code in enumerate(stellar_evolution_codes):
25+
if code == "SeBa":
26+
stellar = Seba()
27+
elif code == "SSE":
28+
stellar = Sse()
29+
elif code == "MESA":
30+
stellar = Mesa()
31+
elif code == "EVtwin":
32+
stellar = Evtwin()
33+
else:
34+
raise ValueError(f"Unknown stellar evolution code {code}")
2735
stellar.parameters.metallicity = z
2836

2937
star = Particles(1)
3038
star.mass = mass
31-
t_end = 6000.0 | units.Myr
39+
# time_end = 6000.0 | units.Myr
3240
stellar.particles.add_particles(star)
33-
attributes = ["temperature", "luminosity","age"]
34-
to_framework = stellar.particles.new_channel_to(star,
35-
attributes=attributes,
36-
target_names=attributes)
41+
attributes = ["temperature", "luminosity", "age"]
42+
to_framework = stellar.particles.new_channel_to(
43+
star, attributes=attributes, target_names=attributes
44+
)
3745
t = [] | units.Myr
3846
T = []
3947
L = []
4048
min_dist_sun = 10000.0
4149
current_time = 3000.0 | units.Myr
4250

43-
print(label[si])
44-
#dt = 50 | units.Myr
45-
#time = 4000 | units.Myr
51+
print(code)
52+
# dt = 50 | units.Myr
53+
# time = 4000 | units.Myr
4654
while stellar:
47-
print(stellar.model_time.value_in(units.Myr))
55+
print(stellar.model_time.value_in(units.Myr))
4856
current_time = current_time + stellar.particles[0].time_step
4957
stellar.evolve_model(current_time)
5058

5159
to_framework.copy()
5260

53-
if star[0].age >= t_end:
61+
if star[0].age >= time_end:
5462
stellar.stop()
5563
stellar = False
5664
else:
57-
L.append((star[0].luminosity - Lstar)/Lstar)
58-
T.append((star[0].temperature - Tstar)/Tstar)
65+
L.append((star[0].luminosity - Lstar) / Lstar)
66+
T.append((star[0].temperature - Tstar) / Tstar)
5967
t.append(star[0].age)
6068

61-
deltaL = numpy.abs((star[0].luminosity - Lstar)/Lstar)
62-
deltaT = numpy.abs((star[0].temperature - Tstar)/Tstar)
69+
deltaL = np.abs((star[0].luminosity - Lstar) / Lstar)
70+
deltaT = np.abs((star[0].temperature - Tstar) / Tstar)
71+
72+
dist = np.sqrt(deltaL * deltaL + deltaT * deltaT)
6373

64-
dist = numpy.sqrt(deltaL*deltaL + deltaT*deltaT)
65-
6674
if min_dist_sun > dist:
6775

6876
min_dist_sun = dist
6977

70-
L_sim_sun = (star[0].luminosity - Lstar)/Lstar
71-
T_sim_sun = (star[0].temperature - Tstar)/Tstar
78+
L_sim_sun = (star[0].luminosity - Lstar) / Lstar
79+
T_sim_sun = (star[0].temperature - Tstar) / Tstar
7280

7381
eta = star[0].age
7482
print(eta)
75-
if si==3:
76-
pyplot.plot(T, L,ls='-', marker=marker[si], markersize=10)
77-
pyplot.scatter(T_sim_sun, L_sim_sun, marker=marker[si],
78-
label=label[si], s=300, lw=1)
83+
if si == 3:
84+
plt.plot(T, L, ls="-", marker=marker[si], markersize=10)
85+
plt.scatter(
86+
T_sim_sun, L_sim_sun, marker=marker[si], label=code, s=300, lw=1
87+
)
7988
else:
80-
pyplot.plot(T, L,ls='-', marker=marker[si], markersize=10)
81-
pyplot.scatter(T_sim_sun, L_sim_sun, marker=marker[si],
82-
label=label[si], s=300, lw=1)
83-
84-
pyplot.legend(scatterpoints=1, loc='best')
85-
86-
save_file = 'fig_SunComparison.png'
87-
pyplot.savefig(save_file)
88-
print('\nSaved figure in file', save_file,'\n')
89-
pyplot.show()
90-
91-
def new_option_parser():
92-
result = OptionParser()
93-
result.add_option("-T", unit=units.K,
94-
dest="Tstar", type="float",default = 5778 |units.K,
95-
help="stellar temperature [%defailt]")
96-
result.add_option("-L", unit=units.LSun,
97-
dest="Lstar", type="float",default = 1 |units.LSun,
98-
help="stellar luminosity [%defailt]")
99-
result.add_option("-m", unit=units.MSun,
100-
dest="mass", type="float",default = 1.0 |units.MSun,
101-
help="stellar mass [%defailt]")
102-
result.add_option("-t", unit=units.Myr,
103-
dest="t_end", type="float", default = 4.6 |units.Gyr,
104-
help="end time of the simulation [%defailt]")
105-
result.add_option("-z", dest="z", type="float", default = 0.02,
106-
help="metalicity [%defailt]")
107-
return result
108-
109-
if __name__ in ('__main__', '__plot__'):
110-
o, arguments = new_option_parser().parse_args()
111-
main(**o.__dict__)
112-
89+
plt.plot(T, L, ls="-", marker=marker[si], markersize=10)
90+
plt.scatter(
91+
T_sim_sun, L_sim_sun, marker=marker[si], label=code, s=300, lw=1
92+
)
93+
94+
plt.legend(scatterpoints=1, loc="best")
95+
96+
save_file = "fig_SunComparison.png"
97+
plt.savefig(save_file)
98+
print("\nSaved figure in file", save_file, "\n")
99+
plt.show()
100+
101+
102+
def new_argument_parser():
103+
parser = argparse.ArgumentParser(
104+
formatter_class=argparse.ArgumentDefaultsHelpFormatter
105+
)
106+
parser.add_argument(
107+
"-T",
108+
"--Tstar",
109+
type=units.K,
110+
default=5778 | units.K,
111+
help="stellar temperature",
112+
)
113+
parser.add_argument(
114+
"-L",
115+
"--Lstar",
116+
type=units.LSun,
117+
default=1 | units.LSun,
118+
help="stellar luminosity",
119+
)
120+
parser.add_argument(
121+
"-m", "--mass", type=units.MSun, default=1.0 | units.MSun, help="stellar mass"
122+
)
123+
parser.add_argument(
124+
"-t",
125+
"--time_end",
126+
type=units.Myr,
127+
default=4.6 | units.Gyr,
128+
help="end time of the simulation",
129+
)
130+
parser.add_argument("-z", type=float, default=0.02, help="metallicity")
131+
return parser
132+
133+
134+
def main():
135+
args = new_argument_parser().parse_args()
136+
plot_solar_comparison(**args.__dict__)
137+
138+
139+
if __name__ == "__main__":
140+
main()

0 commit comments

Comments
 (0)