diff --git a/COSMIC/input/default.cfg b/COSMIC/input/default.cfg index 6749645..d5e06de 100644 --- a/COSMIC/input/default.cfg +++ b/COSMIC/input/default.cfg @@ -1,4 +1,4 @@ -N_stars 10 +N_stars 255 timestep 1 output_times 1 2 5 10 20 50 100 200 500 1000 2000 5000 R0_king 1 diff --git a/COSMIC/source/COSMIC_v3.py b/COSMIC/source/COSMIC_v3.py index 56b1dd9..e6733d4 100644 --- a/COSMIC/source/COSMIC_v3.py +++ b/COSMIC/source/COSMIC_v3.py @@ -19,10 +19,12 @@ params, params_list = config.default() params = config.from_cfg_file(params) stars, binaries, converter = init.king_salpeter_binaries(params) # stars, converter = init.king_salpeter(params) +stars, galaxy_model = init.generate_galaxy(stars, params) channels = init.create_channels() codes = init.gravity_stellar_bridge(converter, params) codes = init.add_stars(codes, stars) codes, channels = init.add_binaries(codes, stars, binaries, channels) -codes, channels = init.commit(codes, stars, channels, params) +# codes, channels = init.commit(codes, stars, channels, params) +codes, channels = init.commit_with_potential(codes, stars, channels, galaxy_model, params) codes, stars, channels = evolve.stellar_gravity_binaries(codes, stars, binaries, channels, params, save_stars = True, save_binaries = True, compute_X_emission = True) # evolve.stellar_gravity(codes, stars, channels, params, save_stars = True) diff --git a/COSMIC/source/COSMIC_v3_config.py b/COSMIC/source/COSMIC_v3_config.py index 9519bc8..b0c9e16 100644 --- a/COSMIC/source/COSMIC_v3_config.py +++ b/COSMIC/source/COSMIC_v3_config.py @@ -11,7 +11,7 @@ import amuse.units.units as u import datetime def params_dic(): - params_list = ["N_stars", "timestep", "output_times", "R0_king", "W0_king", "metallicity_stellar", "M_min_salpeter", "M_max_salpeter", "a_salpeter", "binary_fraction", "mean_period", "std_period", "X_coefficient_OBA", "X_coefficient_GKM", "out_directory", "in_directory", "config", "filename", "workers_stellar", "workers_gravity", "format_type"] + params_list = ["N_stars", "timestep", "output_times", "R0_king", "W0_king", "metallicity_stellar", "M_min_salpeter", "M_max_salpeter", "a_salpeter", "binary_fraction", "mean_period", "std_period", "X_coefficient_OBA", "X_coefficient_GKM", "cluster_position_x", "cluster_position_y", "cluster_position_z", "out_directory", "in_directory", "config", "filename", "workers_stellar", "workers_gravity", "format_type"] params = {i:0 for i in params_list} return params, params_list @@ -42,6 +42,9 @@ def default(): params["std_period"] = np.float32(2.4) # Orbital period deviation of binary systems params["X_coefficient_OBA"] = np.float32(1.4e-7) # X-ray luminosity coefficient for O, B and A type stars params["X_coefficient_GKM"] = np.float32(1.4e27) # X-ray luminosity coefficient for G, K and M type stars + params["cluster_position_x"] = np.float32(7000)|u.pc # Position of the cluster in the disk along the sun-center axis + params["cluster_position_y"] = np.float32(0)|u.pc # Position of the cluster in the disk in the perpendicular direction + params["cluster_position_z"] = np.float32(0)|u.pc # Position of the cluster above or below the disk params["out_directory"] = "./output/" params["in_directory"] = "./input/" params["config"] = "default.cfg" @@ -68,6 +71,9 @@ def change_params(params, cfg_param, cfg_values): elif cfg_param == "std_period": params["std_period"] = np.float32(cfg_values[0]) elif cfg_param == "X_coefficient_OBA": params["X_coefficient_OBA"] = np.float32(cfg_values[0]) elif cfg_param == "X_coefficient_GKM": params["X_coefficient_GKM"] = np.float32(cfg_values[0]) + elif cfg_param == "cluster_position_x": params["cluster_position_x"] = np.float32(cfg_values[0])|u.pc + elif cfg_param == "cluster_position_y": params["cluster_position_y"] = np.float32(cfg_values[0])|u.pc + elif cfg_param == "cluster_position_z": params["cluster_position_z"] = np.float32(cfg_values[0])|u.pc elif cfg_param == "out_directory": params["out_directory"] = cfg_values[0].replace("\n", "") elif cfg_param == "in_directory": params["in_directory"] = cfg_values[0].replace("\n", "") elif cfg_param == "config": params["config"] = cfg_values[0].replace("\n", "") diff --git a/COSMIC/source/COSMIC_v3_coordinates.py b/COSMIC/source/COSMIC_v3_coordinates.py new file mode 100644 index 0000000..9996883 --- /dev/null +++ b/COSMIC/source/COSMIC_v3_coordinates.py @@ -0,0 +1,28 @@ +""" +* COSMIC - COORDINATES +* Version 3 +@ Yaël Moussouni +@ Unistra, P&E, MSc1-MoFP +@ Observatory of Strasbourg (intern) +""" + +import numpy as np +import amuse.units.units as u_amuse +import astropy.units as u_astropy + +import astropy.coordinates as c + +def xyz2radecdist(stars): + X = stars.x.value_in(u_amuse.pc) + Y = stars.y.value_in(u_amuse.pc) + Z = stars.z.value_in(u_amuse.pc) + Galactocentric = c.SkyCoord(x = X * u_astropy.pc, y = Y * u_astropy.pc, z = Z * u_astropy.pc, frame=c.Galactocentric) + ICRS = Galactocentric.transform_to(c.ICRS) + RA = ICRS.ra.deg + DEC = ICRS.dec.deg + DIST = ICRS.distance.pc + stars.ra = RA | u_amuse.deg + stars.dec = DEC | u_amuse.deg + stars.dist = DIST | u_amuse.pc + return stars + diff --git a/COSMIC/source/COSMIC_v3_emission.py b/COSMIC/source/COSMIC_v3_emission.py index 76039db..6a59ebe 100644 --- a/COSMIC/source/COSMIC_v3_emission.py +++ b/COSMIC/source/COSMIC_v3_emission.py @@ -23,5 +23,31 @@ def X_emission(stars, params): elif stars[i].temperature.value_in(u.K) < 6000: # G K M stars[i].X_luminosity = params["X_coefficient_GKM"] * (stars[i].X_vsini.value_in(u.km * u.s**(-1)))**(1.9) | u.erg * u.s**(-1) else: # F - stars[i].X_luminosity = 0 + stars[i].X_luminosity = 0 | u.erg * u.s**(-1) + return stars + +def plasma_temperature(stars, params): + for i in range(len(stars)): + if stars[i].temperature.value_in(u.K) > 7300 : # O B A + if stars[i].age.value_in(u.Myr) < 150: # young + stars[i].X_temperature_0 = 0.6 | u.kilo(u.eV) + stars[i].X_temperature_1 = 0.0 | u.kilo(u.eV) + else: # old + stars[i].X_temperature_0 = 0.5 | u.kilo(u.eV) + stars[i].X_temperature_1 = 0.0 | u.kilo(u.eV) + elif stars[i].temperature.value_in(u.K) < 6000: # G K M + if stars[i].age.value_in(u.Myr) < 150: # young + stars[i].X_temperature_0 = 0.4 | u.kilo(u.eV) + stars[i].X_temperature_1 = 1.0 | u.kilo(u.eV) + else: # old + stars[i].X_temperature_0 = 0.2 | u.kilo(u.eV) + stars[i].X_temperature_1 = 0.8 | u.kilo(u.eV) + else: # F + if stars[i].age.value_in(u.Myr) < 150: # young + stars[i].X_temperature_0 = 0.6 | u.kilo(u.eV) + stars[i].X_temperature_1 = 0.0 | u.kilo(u.eV) + else: # old + stars[i].X_temperature_0 = 0.5 | u.kilo(u.eV) + stars[i].X_temperature_1 = 0.0 | u.kilo(u.eV) + return stars diff --git a/COSMIC/source/COSMIC_v3_evolve.py b/COSMIC/source/COSMIC_v3_evolve.py index a38e148..dec35fa 100644 --- a/COSMIC/source/COSMIC_v3_evolve.py +++ b/COSMIC/source/COSMIC_v3_evolve.py @@ -10,6 +10,7 @@ import numpy as np import amuse.units.units as u import COSMIC_v3_output as output import COSMIC_v3_emission as emission +import COSMIC_v3_coordinates as coords def stellar_gravity(codes, stars, channels, params, save_stars = False): """ @@ -32,9 +33,11 @@ def stellar_gravity_binaries(codes, stars, binaries, channels, params, save_star print("\033[36m"+"t: "+"\033[0m"+"{:02.2f} Myr/{:02.2f} Myr".format(t.value_in(u.Myr),np.max(params["output_times"].value_in(u.Myr)))) codes["b"].evolve_model(t) channels.copy() + stars = coords.xyz2radecdist(stars) if compute_X_emission: stars = emission.rotation(stars, params) stars = emission.X_emission(stars, params) + stars = emission.plasma_temperature(stars, params) if save_stars: output.save_stars(stars, params, t.value_in(u.Myr)) if save_binaries: output.save_binaries(binaries, params, t.value_in(u.Myr)) diff --git a/COSMIC/source/COSMIC_v3_galaxy.py b/COSMIC/source/COSMIC_v3_galaxy.py index e69de29..b646e44 100644 --- a/COSMIC/source/COSMIC_v3_galaxy.py +++ b/COSMIC/source/COSMIC_v3_galaxy.py @@ -0,0 +1,87 @@ +""" +* COSMIC - GALAXY +* Version 3 +@ Yaël Moussouni +@ Unistra, P&E, MSc1-MoFP +@ Observatory of Strasbourg (Intern) +""" + +import numpy as np +import matplotlib.pyplot as plt +import amuse.units.units as u + +class MilkyWay_AMUSE(object): + # THIS IS A COPY OF https://github.com/amusecode/amuse/blob/cdd21cc5cb06e40ccf5ecb86d513d211634e2689/examples/textbook/solar_cluster_in_galaxy_potential.py#L37 + + def __init__(self, Mb=1.40592e10| u.MSun, + Md=8.5608e10| u.MSun, + Mh=1.07068e11 | u.MSun ): + self.Mb= Mb + self.Md= Md + self.Mh= Mh + + def get_potential_at_point(self,eps,x,y,z): + r=(x**2+y**2+z**2)**0.5 + R= (x**2+y**2)**0.5 + # buldge + b1= 0.3873 |u.kpc + pot_bulge= -u.constants.G*self.Mb/(r**2+b1**2)**0.5 + # disk + a2= 5.31 |u.kpc + b2= 0.25 |u.kpc + pot_disk = \ + -u.constants.G*self.Md/(R**2 + (a2+ (z**2+ b2**2)**0.5 )**2 )**0.5 + #halo + a3= 12.0 |u.kpc + cut_off=100 |u.kpc + d1= r/a3 + c=1+ (cut_off/a3)**1.02 + pot_halo= -u.constants.G*(self.Mh/a3)*d1**1.02/(1+ d1**1.02) \ + - (u.constants.G*self.Mh/(1.02*a3))\ + * (-1.02/c +numpy.log(c) + 1.02/(1+d1**1.02) \ + - numpy.log(1.0 +d1**1.02) ) + return 2*(pot_bulge+pot_disk+ pot_halo) # multiply by 2 because it + # is a rigid potential + + + def get_gravity_at_point(self, eps, x,y,z): + r= (x**2+y**2+z**2)**0.5 + R= (x**2+y**2)**0.5 + #bulge + b1= 0.3873 |u.kpc + force_bulge= -u.constants.G*self.Mb/(r**2+b1**2)**1.5 + #disk + a2= 5.31 |u.kpc + b2= 0.25 |u.kpc + d= a2+ (z**2+ b2**2)**0.5 + force_disk=-u.constants.G*self.Md/(R**2+ d**2 )**1.5 + #halo + a3= 12.0 |u.kpc + d1= r/a3 + force_halo= -u.constants.G*self.Mh*d1**0.02/(a3**2*(1+d1**1.02)) + + ax= force_bulge*x + force_disk*x + force_halo*x/r + ay= force_bulge*y + force_disk*y + force_halo*y/r + az= force_bulge*z + force_disk*d*z/(z**2 + b2**2)**0.5 + force_halo*z/r + + return ax,ay,az + + def vel_circ(self, r ): + z=0 | u.kpc + b1= 0.3873 |u.kpc + a2= 5.31 |u.kpc + b2= 0.25 |u.kpc + a3= 12.0 |u.kpc + + rdphi_b = u.constants.G*self.Mb*r**2/(r**2+b1**2)**1.5 + rdphi_d =u.constants.G*self.Md*r**2/(r**2+(a2+(z**2+b2**2)**0.5)**2 )**1.5 + rdphi_h = u.constants.G*self.Mh*(r/a3)**0.02*r/(a3**2*(1+(r/a3)**1.02)) + + vel_circb = rdphi_b + vel_circd = rdphi_d + vel_circh = rdphi_h + + return (vel_circb+ vel_circd+ vel_circh)**0.5 + + def stop(self): + return diff --git a/COSMIC/source/COSMIC_v3_init.py b/COSMIC/source/COSMIC_v3_init.py index ae7f33c..248b7bd 100644 --- a/COSMIC/source/COSMIC_v3_init.py +++ b/COSMIC/source/COSMIC_v3_init.py @@ -10,6 +10,7 @@ import numpy as np import amuse.units.units as u import COSMIC_v3_galaxy as galaxy +import COSMIC_v3_coordinates as coords from amuse.units import nbody_system from amuse.ic.kingmodel import new_king_model @@ -97,6 +98,7 @@ def king_salpeter_binaries(params): stars.ra = 0 | u.deg stars.dec = 0 | u.deg + stars.dist = 0 | u.pc stars.X_v = 0 | u.km * u.s**(-1) stars.X_vsini = 0 | u.km * u.s**(-1) stars.X_temperature_0 = 0 | u.kilo(u.eV) @@ -120,6 +122,8 @@ def king_salpeter_binaries(params): semi_major_axis, eccentricity = orbital_parameters(params, binaries[i].mass_tot) binaries[i].semi_major_axis = semi_major_axis binaries[i].eccentricity = eccentricity + + stars.move_to_center() return stars, binaries, converter @@ -140,9 +144,10 @@ def generate_galaxy(stars, params): stars.x += params["cluster_position_x"] stars.y += params["cluster_position_y"] stars.z += params["cluster_position_z"] - phi = np.arctan2(stars.center_of_mass().y, stars.center_of_mass().x) + phi = np.arctan2(stars.center_of_mass().y.value_in(u.pc), stars.center_of_mass().x.value_in(u.pc)) stars.vx += - vc * np.sin(phi) stars.vy += vc * np.cos(phi) + stars = coords.xyz2radecdist(stars) return stars, galaxy_model def add_stars(codes, stars):