Add files via upload

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

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)

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", "")

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
+

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

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))
 	 

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

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):