From 80c6d54f008d97f513646526f50916674aa00a76 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ya=C3=ABl=20M?= <154762804+Yael-II@users.noreply.github.com>
Date: Fri, 21 Jun 2024 16:19:00 +0200
Subject: [PATCH] Add files via upload

---
 COSMIC/COSMIC_v3.py             |  30 +++++
 COSMIC/COSMIC_v3_config.py      | 102 ++++++++++++++++
 COSMIC/COSMIC_v3_coordinates.py |  28 +++++
 COSMIC/COSMIC_v3_emission.py    |  53 +++++++++
 COSMIC/COSMIC_v3_evolve.py      |  47 ++++++++
 COSMIC/COSMIC_v3_galaxy.py      |  87 ++++++++++++++
 COSMIC/COSMIC_v3_init.py        | 198 ++++++++++++++++++++++++++++++++
 COSMIC/COSMIC_v3_output.py      |  20 ++++
 8 files changed, 565 insertions(+)
 create mode 100644 COSMIC/COSMIC_v3.py
 create mode 100644 COSMIC/COSMIC_v3_config.py
 create mode 100644 COSMIC/COSMIC_v3_coordinates.py
 create mode 100644 COSMIC/COSMIC_v3_emission.py
 create mode 100644 COSMIC/COSMIC_v3_evolve.py
 create mode 100644 COSMIC/COSMIC_v3_galaxy.py
 create mode 100644 COSMIC/COSMIC_v3_init.py
 create mode 100644 COSMIC/COSMIC_v3_output.py

diff --git a/COSMIC/COSMIC_v3.py b/COSMIC/COSMIC_v3.py
new file mode 100644
index 0000000..e6733d4
--- /dev/null
+++ b/COSMIC/COSMIC_v3.py
@@ -0,0 +1,30 @@
+"""
+* COSMIC : Cluster Orbital SysteM Integration Code
+* Version 3
+@ Yaël Moussouni
+@ Unistra, P&E, MSc1-MoFP
+@ Observatory of Strasbourg (Intern)
+"""
+
+import numpy as np
+import amuse.units.units as u
+# import matplotlib.pyplot as plt
+import COSMIC_v3_config as config
+import COSMIC_v3_init as init
+import COSMIC_v3_evolve as evolve
+# import COSMIC_v3_output as output
+# import COSMIC_v3_emission as emission
+
+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_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/COSMIC_v3_config.py b/COSMIC/COSMIC_v3_config.py
new file mode 100644
index 0000000..b0c9e16
--- /dev/null
+++ b/COSMIC/COSMIC_v3_config.py
@@ -0,0 +1,102 @@
+"""
+* COSMIC - CONFIG
+* Version 3
+@ Yaël Moussouni
+@ Unistra, P&E, MSc1-MoFP
+@ Observatory of Strasbourg (Intern)
+"""
+
+import numpy as np
+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", "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
+
+def verbatim(params):
+    print("\033[94m"+"Info: Simulation time of "+"\033[0m"+"{:.2f} Myr".format(np.max(params["output_times"].value_in(u.Myr))))
+    print("\033[94m"+"Info: Simulation of "+"\033[0m"+"{:d} stars".format(params["N_stars"]))
+    print("\033[94m"+"Info: Simulation will be saved with the prefix "+"\033[0m"+"{}".format(params["filename"]))
+    print("\033[94m"+"Info: Gravity: "+"\033[0m"+"{:d} workers".format(params["workers_gravity"])+"\033[94m"+", Stellar: "+"\033[0m"+"{:d} workers".format(params["workers_stellar"]))
+    return 0
+
+def default():
+    params, params_list = params_dic()
+
+    date = datetime.datetime.now()
+    date = "{:04d}-{:02d}-{:02d}_{:02d}h{:02d}".format(date.year, date.month, date.day, date.hour, date.minute)
+
+    params["N_stars"] = np.int32(10) # Number of particle
+    params["timestep"] = np.float32(1)|u.Myr # Time increment
+    params["output_times"] =  np.float32([1,2,5,10,20,50,100,200,500,1000,2000,5000])|u.Myr # Output intants
+    params["R0_king"] = np.float32(1)|u.pc # Cluster virial radius
+    params["W0_king"] = np.float32(9) # Cluster King's parameter
+    params["metallicity_stellar"] = np.float32(0.02) # Metallicity
+    params["M_min_salpeter"] = np.float32(0.1)|u.MSun # Minimum mass for Salpeter's law
+    params["M_max_salpeter"] = np.float32(125)|u.MSun # Maximum mass for Salpeter's law
+    params["a_salpeter"] = np.float32(-2.35) # Salpeter's law coefficient
+    params["binary_fraction"] = np.float32(0.13) # Fraction of binary stars
+    params["mean_period"] = np.float32(4.54) # Mean orbital period of binary systems (log10(T/days))
+    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"
+    params["filename"] = "COSMIC_{}".format(date)
+    params["workers_stellar"] = np.int32(1)
+    params["workers_gravity"] = np.int32(1)
+    params["format_type"] = "csv"
+    return params, params_list
+
+def change_params(params, cfg_param, cfg_values):
+    if cfg_param == "":
+        return params
+    elif cfg_param == "N_stars": params["N_stars"] = np.int32(cfg_values[0])
+    elif cfg_param == "timestep": params["timestep"] = np.float32(cfg_values[0])|u.Myr
+    elif cfg_param == "output_times": params["output_times"] = np.float32(cfg_values)|u.Myr
+    elif cfg_param == "R0_king": params["R0_king"] = np.float32(cfg_values[0])|u.pc
+    elif cfg_param == "W0_king": params["W0_king"] = np.float32(cfg_values[0])
+    elif cfg_param == "metallicity_stellar": params["metallicity_stellar"] = np.float32(cfg_values[0])
+    elif cfg_param == "M_min_salpeter": params["M_min_salpeter"] = np.float32(cfg_values[0])|u.MSun
+    elif cfg_param == "M_max_salpeter": params["M_max_salpeter"] = np.float32(cfg_values[0])|u.MSun
+    elif cfg_param == "a_salpeter": params["a_salpeter"] = np.float32(cfg_values[0])
+    elif cfg_param == "binary_fraction": params["binary_fraction"] = np.float32(cfg_values[0])
+    elif cfg_param == "mean_period": params["mean_period"] = np.float32(cfg_values[0])
+    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", "")
+    elif cfg_param == "filename": params["filename"] = cfg_values[0].replace("\n", "")
+    elif cfg_param == "workers_stellar": params["workers_stellar"] = np.int32(cfg_values[0])
+    elif cfg_param == "workers_gravity": params["workers_gravity"] = np.int32(cfg_values[0])
+    elif cfg_param == "format_type": params["format_type"] = cfg_values[0].replace("\n", "")
+    else:
+        print("\033[93m"+"Warning: cannot understand the parameter in config: "+"\033[0m" + cfg_param)
+    return params
+            
+
+def from_cfg_file(params):
+    cfg_path = params["in_directory"] + params["config"]
+    with open(cfg_path, "r") as cfg_file:
+        for cfg_line in cfg_file.readlines():
+            cfg_split = cfg_line.split(" ")
+            if len(cfg_split) == 1:
+                print("\033[93m"+"Warning: cannot understand the line in config: "+"\033[0m" + cfg_line)
+                continue
+            cfg_param = cfg_split[0]
+            cfg_values = cfg_split[1:]
+            params = change_params(params, cfg_param, cfg_values) 
+    verbatim(params)
+    return params
+
diff --git a/COSMIC/COSMIC_v3_coordinates.py b/COSMIC/COSMIC_v3_coordinates.py
new file mode 100644
index 0000000..9996883
--- /dev/null
+++ b/COSMIC/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/COSMIC_v3_emission.py b/COSMIC/COSMIC_v3_emission.py
new file mode 100644
index 0000000..6a59ebe
--- /dev/null
+++ b/COSMIC/COSMIC_v3_emission.py
@@ -0,0 +1,53 @@
+"""
+* COSMIC - EMISSION
+* Version 3
+@ Yaël Moussouni
+@ Unistra, P&E, MSc1-MoFP
+@ Observatory of Strasbourg (Intern)
+"""
+
+import numpy as np
+import amuse.units.units as u
+
+def rotation(stars, params):
+    N = params["N_stars"] + np.int32(np.round(params["N_stars"] * params["binary_fraction"]))
+    stars.X_v = np.random.uniform(0, 10, N) | u.km * u.s**(-1)
+    p = np.random.uniform(0, 1, N)
+    stars.X_vsini = stars.X_v * np.sqrt(2*p-p**2)
+    return stars
+
+def X_emission(stars, params):
+    for i in range(len(stars)):
+        if stars[i].temperature.value_in(u.K) > 7300: # O B A
+            stars[i].X_luminosity = params["X_coefficient_OBA"] * stars[i].luminosity.in_(u.erg * u.s**(-1))
+        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 | 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/COSMIC_v3_evolve.py b/COSMIC/COSMIC_v3_evolve.py
new file mode 100644
index 0000000..dec35fa
--- /dev/null
+++ b/COSMIC/COSMIC_v3_evolve.py
@@ -0,0 +1,47 @@
+"""
+* COSMIC - EVOLVE
+* Version 3
+@ Yaël Moussouni
+@ Unistra, P&E, MSc1-MoFP
+@ Observatory of Strasbourg (Intern)
+"""
+
+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):
+    """
+    Evolve a model with both gravity and stellar evolution.
+    """
+    for t in params["output_times"]:
+        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["g"].evolve_model(t)
+        channels.copy()
+        if save_stars: output.save_stars(stars, params, t.value_in(u.Myr))
+	 
+    print("\033[36m"+"Simulation completed"+"\033[0m")
+    return 0
+
+def stellar_gravity_binaries(codes, stars, binaries, channels, params, save_stars = False, save_binaries = False, compute_X_emission = False):
+    """
+    Evolve a model with both gravity and stellar evolution.
+    """
+    for t in params["output_times"]:
+        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))
+	 
+    print("\033[36m"+"Simulation completed"+"\033[0m")
+    return codes, stars, channels
+
+
diff --git a/COSMIC/COSMIC_v3_galaxy.py b/COSMIC/COSMIC_v3_galaxy.py
new file mode 100644
index 0000000..b646e44
--- /dev/null
+++ b/COSMIC/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/COSMIC_v3_init.py b/COSMIC/COSMIC_v3_init.py
new file mode 100644
index 0000000..768b9c3
--- /dev/null
+++ b/COSMIC/COSMIC_v3_init.py
@@ -0,0 +1,198 @@
+"""
+* COSMIC - INIT
+* Version 3 
+@ Yaël Moussouni
+@ Unistra, P&E, MSc1-MoFP
+@ Observatory of Strasbourg (Intern)
+"""
+
+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
+from amuse.ic.salpeter import new_salpeter_mass_distribution
+# from amuse.community.huayno.interface import Huayno # FIXME DOES NOT WORK
+from amuse.community.ph4.interface import Ph4
+from amuse.community.seba.interface import SeBa
+from amuse.couple.bridge import Bridge
+from amuse.datamodel.particles import Channels
+from amuse.datamodel import Particles
+
+# def IMF_Besancon(M, alpha_1=1.6, alpha_2=3.0):
+#     return (M.value_in(u.MSun))**(alpha_1-1) * ((M.value_in(u.MSun))<1) + (((M.value_in(u.MSun))-1)**(alpha_2-1) + 1) * ((M.value_in(u.MSun))>=1)
+
+# def Besancon_disk(params):
+#     """
+#     Generate two IMF of stars and binaries
+#     """
+#     M_1 = []
+#     M_2 = []
+#     N_min = 0
+#     N_max = IMF_Besancon(params["M_max_besancon"], params["a_besancon_1"], params["a_besancon_2"])
+#     while len(M_1) < params["N_stars"]:
+#         X = np.random.uniform(params["M_min_besancon"].value_in(u.MSun), params["M_max_besancon"].value_in(u.MSun))|u.MSun
+#         Y = np.random.uniform(N_min, N_max)
+#         if Y > IMF_Besancon(X, params["a_besancon_1"], params["a_besancon_2"]):
+#             M_1.append(X)
+#     while len(M_2) < params["N_stars"] * params["binary_fraction"]:
+#         X = np.random.uniform(params["M_min_besancon"].value_in(u.MSun), params["M_max_besancon"].value_in(u.MSun))|u.MSun
+#         Y = np.random.uniform(N_min, N_max)
+#         if Y > IMF_Besancon(X, params["a_besancon_1"], params["a_besancon_2"]):
+#             M_2.append(X)
+#     return M_1, M_2
+
+def orbital_parameters(params, m_tot):
+    x_0 = params["mean_period"] * np.log(2)
+    sigma_0 = params["std_period"] * np.log(2)
+    T = np.random.lognormal(mean = x_0, sigma = sigma_0)
+    a = (T**2 * u.constants.G.value_in(u.au**3*u.day**(-2)*u.MSun**(-1)) * m_tot.value_in(u.MSun) / (4*np.pi**2) )**(1/3) | u.au
+    e = 1
+    return a, e
+
+def king_salpeter(params):
+    """
+    Create a cluster of N stars, with a King's law of parameter W0, with a (virial) radius R, stellar masses following a Salpeter law ranging between M_min and M_max with exponent a_salpeter.
+    """
+    M_stars = new_salpeter_mass_distribution(params["N_stars"], mass_min=params["M_min_salpeter"], mass_max=params["M_max_salpeter"], alpha=params["a_salpeter"])
+    converter = nbody_system.nbody_to_si(M_stars.sum(), params["R0_king"])
+    stars = new_king_model(params["N_stars"], params["W0_king"], converter)
+    stars.mass = M_stars
+    stars.zams_mass = M_stars
+    return stars, converter
+
+def king_salpeter_binaries(params):
+    """
+    
+    """
+    M_primary = new_salpeter_mass_distribution(params["N_stars"], mass_min=params["M_min_salpeter"], mass_max=params["M_max_salpeter"], alpha=params["a_salpeter"])
+    M_secondary = new_salpeter_mass_distribution(np.int32(np.round(params["N_stars"]*params["binary_fraction"])), mass_min=params["M_min_salpeter"], mass_max=params["M_max_salpeter"], alpha=params["a_salpeter"])
+
+    M_stars = np.concatenate([M_primary, M_secondary])
+    N_binaries = len(M_secondary)
+    converter = nbody_system.nbody_to_si(M_stars.sum(), params["R0_king"])
+    stars_primary = new_king_model(params["N_stars"], params["W0_king"], converter)
+    stars_primary.mass = M_primary
+    stars_primary.zams_mass = M_primary
+     
+    stars_secondary = Particles(N_binaries)
+    for i in range(N_binaries):
+        stars_secondary[i].mass = M_secondary[i]
+        stars_secondary[i].radius = 0|u.m
+        stars_secondary[i].vx = stars_primary[i].vx
+        stars_secondary[i].vy = stars_primary[i].vy
+        stars_secondary[i].vz = stars_primary[i].vz
+        stars_secondary[i].x = stars_primary[i].x
+        stars_secondary[i].y = stars_primary[i].y
+        stars_secondary[i].z = stars_primary[i].z
+        stars_secondary[i].zams_mass = M_secondary[i]
+
+    stars = Particles()
+    binaries = Particles(N_binaries)
+
+    stars.add_particles(stars_primary)
+    stars.add_particles(stars_secondary)
+
+    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)
+    stars.X_temperature_1 = 0 | u.kilo(u.eV)
+    stars.X_luminosity = 0 | u.erg * u.s**(-1)
+
+    if N_binaries == 0:
+        binaries.child1 = [0]
+        binaries.child2 = [0]
+        binaries.mass1 = 0
+        binaries.mass2 = 0
+        binaries.mass_tot = 0
+    else:
+        binaries.child1 = list(stars_primary[:N_binaries])
+        binaries.child2 = list(stars_secondary)
+        binaries.mass1 = list(stars_primary[:N_binaries].mass)
+        binaries.mass2 = list(stars_secondary.mass.in_(u.kg))
+        binaries.mass_tot = binaries.mass1 + binaries.mass2
+
+    for i in range(N_binaries):
+        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
+
+
+def gravity_stellar_bridge(converter, params):
+    gravity = Ph4(converter, number_of_workers=params["workers_gravity"]) # TODO A CHANGER SI BESOIN
+    stellar = SeBa(number_of_worker=params["workers_stellar"])
+    bridge = Bridge()
+    codes = {"g":gravity, "s":stellar, "b":bridge}
+    return codes
+
+def create_channels():
+    channels = Channels()
+    return channels
+
+def generate_galaxy(stars, params):
+    galaxy_model = galaxy.MilkyWay_AMUSE()
+    stars.x += params["cluster_position_x"]
+    stars.y += params["cluster_position_y"]
+    stars.z += params["cluster_position_z"]
+    vc = galaxy_model.vel_circ(stars.center_of_mass().length()).in_(u.km * u.s**(-1))
+    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):
+    codes["g"].particles.add_particles(stars)
+    codes["s"].particles.add_particles(stars)
+    return codes
+
+def add_binaries(codes, stars, binaries, channels):
+    codes["s"].binaries.add_particles(binaries)
+    channels.add_channel(codes["s"].binaries.new_channel_to(stars))
+    return codes, channels
+
+def commit(codes, stars, channels, params):
+    codes["s"].set_metallicity(params["metallicity_stellar"])
+    codes["s"].commit_particles()
+    codes["g"].commit_particles()
+    # codes["g"].initialize_code() # NOT REQUIRED WITH Ph4
+
+    codes["b"].add_system(codes["g"])
+    codes["b"].add_system(codes["s"])
+
+    codes["b"].synchronize_model()
+    codes["b"].timestep = params["timestep"]
+
+    codes["b"].channels.add_channel(codes["s"].particles.new_channel_to(codes["g"].particles, attributes=["mass", "radius"])) 
+
+    channels.add_channel(codes["s"].particles.new_channel_to(stars))
+    channels.add_channel(codes["g"].particles.new_channel_to(stars))
+    return codes, channels
+
+def commit_with_potential(codes, stars, channels, galaxy_model, params):
+    codes["s"].set_metallicity(params["metallicity_stellar"])
+    codes["s"].commit_particles()
+    codes["g"].commit_particles()
+    # codes["g"].initialize_code() # NOT REQUIRED WITH Ph4
+
+    codes["b"].add_system(codes["g"], (galaxy_model,))
+    codes["b"].add_system(codes["s"])
+
+    codes["b"].synchronize_model()
+    codes["b"].timestep = params["timestep"]
+
+    codes["b"].channels.add_channel(codes["s"].particles.new_channel_to(codes["g"].particles, attributes=["mass", "radius"])) 
+
+    channels.add_channel(codes["s"].particles.new_channel_to(stars))
+    channels.add_channel(codes["g"].particles.new_channel_to(stars))
+    return codes, channels
+
diff --git a/COSMIC/COSMIC_v3_output.py b/COSMIC/COSMIC_v3_output.py
new file mode 100644
index 0000000..0841c5f
--- /dev/null
+++ b/COSMIC/COSMIC_v3_output.py
@@ -0,0 +1,20 @@
+"""
+* COSMIC - SAVE
+* Version 3
+@ Yaël Moussouni
+@ Unistra, P&E, MSc1-MoFP
+@ Observatory of Strasbourg (Intern)
+"""
+
+
+import numpy as np
+
+from amuse.io.base import write_set_to_file
+
+def save_stars(stars, params, t=0):
+    write_set_to_file(set=stars, filename=params["out_directory"]+params["filename"]+"_stars_{:04d}_Myr.{}".format(np.int32(t), params["format_type"]), format=params["format_type"])
+    return 0
+
+def save_binaries(binaries, params, t=0):
+    write_set_to_file(set=binaries, filename=params["out_directory"]+params["filename"]+"_binaries_{:04d}_Myr.{}".format(np.int32(t), params["format_type"]), format=params["format_type"])
+    return 0