commit e5bb7b3ebd4c2840a4987e7f3928832dcfd16b06
Author: Yaël M <154762804+Yael-II@users.noreply.github.com>
Date:   Wed Jun 19 15:45:57 2024 +0200

    Add files via upload

diff --git a/COSMIC/cosmic.sh b/COSMIC/cosmic.sh
new file mode 100644
index 0000000..4678bfb
--- /dev/null
+++ b/COSMIC/cosmic.sh
@@ -0,0 +1,6 @@
+if [[ "$VIRTUAL_ENV" == "/Users/yael/.venv/amuse" ]]
+then
+     mpirun -n 1 python source/COSMIC_v3.py
+else
+     echo "please activate amuse"
+fi
diff --git a/COSMIC/input/default.cfg b/COSMIC/input/default.cfg
new file mode 100644
index 0000000..6749645
--- /dev/null
+++ b/COSMIC/input/default.cfg
@@ -0,0 +1,20 @@
+N_stars 10
+timestep 1
+output_times 1 2 5 10 20 50 100 200 500 1000 2000 5000
+R0_king 1
+W0_king 9
+metallicity_stellar 0.02
+M_min_salpeter 0.1
+M_max_salpeter 125
+a_salpeter -2.35
+binary_fraction 0.13
+mean_period 4.54
+std_period 2.4
+X_coefficient_OBA 1.4e-7
+X_coefficient_GKM 1.4e27
+out_directory ./output/
+in_directory ./input/
+config default.cfg
+workers_stellar 1 
+workers_gravity 1
+format_type csv
diff --git a/COSMIC/old/COSMIC-VIS_v2.py b/COSMIC/old/COSMIC-VIS_v2.py
new file mode 100644
index 0000000..8f551b2
--- /dev/null
+++ b/COSMIC/old/COSMIC-VIS_v2.py
@@ -0,0 +1,385 @@
+"""
+* COSMIC-VIS : Cluster Orbital SysteM Integration Code - Visualisation
+* Version 2 - 2024-02-27
+@ Yaël Moussouni
+@ Unistra, P&E, MSc1-MoFP
+@ Observatory of Strasbourg (Intern)
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.animation as animation
+import os
+
+# * Parameters
+if "YII" in plt.style.available: plt.style.use("YII")
+
+if "YII_light" in plt.style.available: light = "YII_light"
+else: light = "default"
+
+# * Variables
+
+in_directory = "./COSMIC_output/"
+out_directory = "../Figs/COSMIC_figs/"
+
+cluster_suffix = "_cluster.csv"
+stars_suffix = "_stars.csv"
+
+YES = ["Y", "YES", "1", "OUI"]
+MARKER = ["o", "s", "D", "v", "^", "<", ">"]
+
+Tmin = 2000  # K - Temperature truncation
+Tmax = 10000 # K - Temperature truncation
+
+# * Functions
+
+def density_hist(X,Y,M,T,bins,i):
+   x_cm = np.average(X[i], weights=M[i])
+   y_cm = np.average(Y[i], weights=M[i])
+   Radius = np.sqrt((X[i]-x_cm)**2 + (Y[i]-y_cm)**2)
+   Rho = np.zeros(len(bins)-1)
+   for k in range(len(bins)-1):
+      r1 = bins[k]
+      r2 = bins[k+1]
+      w = np.where((Radius > r1) & (Radius < r2))
+      Rho[k] = np.sum(M[i,w])/(4*np.pi*r2**3/3 - 4*np.pi*r1**3/3)
+   return Rho
+
+# * Listing files for selection 
+filelist = np.sort(os.listdir(in_directory))
+filelist = [f for f in filelist if f[0] != "."]
+dates = []
+for f in filelist:
+   if len(f) == len("YYYY-MM-DD_HH:MM_cluster.csv") or len(f) == len("YYYY-MM-DD_HH:MM_stars.csv"):
+      date = f[0:16]
+      if not date in dates:
+         dates.append(date)
+   
+# * File selection
+print("\033[36m"+"Available dates:"+"\033[0m")
+for i in range(len(dates)):
+   print("\033[36m"+"\t{}: ".format(i+1) + "\033[0m" + "{}".format(dates[i].replace("h", ":").replace("_", " ")))
+j = int(input("\033[32m"+"Select a date number: "+"\033[0m"))-1
+print("")
+date_prefix = dates[j]
+
+# * Data collection
+cluster = np.loadtxt(in_directory+date_prefix+cluster_suffix, skiprows=1, comments="#", delimiter=",")
+stars = np.loadtxt(in_directory+date_prefix+stars_suffix, skiprows=1, comments="#", delimiter=",")
+
+T = cluster[:,1]
+E = cluster[:,2]
+DE = cluster[:,3]
+L = cluster[:,4]
+DL = cluster[:,5]
+
+M = np.zeros(np.int32(stars[-1,0:2])+1)
+X = np.zeros(np.int32(stars[-1,0:2])+1)
+Y = np.zeros(np.int32(stars[-1,0:2])+1)
+Z = np.zeros(np.int32(stars[-1,0:2])+1)
+VX = np.zeros(np.int32(stars[-1,0:2])+1)
+VY = np.zeros(np.int32(stars[-1,0:2])+1)
+VZ = np.zeros(np.int32(stars[-1,0:2])+1)
+Lumi = np.zeros(np.int32(stars[-1,0:2])+1)
+Radi = np.zeros(np.int32(stars[-1,0:2])+1)
+Temp = np.zeros(np.int32(stars[-1,0:2])+1)
+
+for star in stars:
+   i = np.int32(star[0])
+   j = np.int32(star[1])
+   M[i,j] = star[2]
+   X[i,j] = star[3]
+   Y[i,j] = star[4]
+   Z[i,j] = star[5]
+   VX[i,j] = star[6]
+   VY[i,j] = star[7]
+   VZ[i,j] = star[8]
+   Lumi[i,j] = star[9]
+   Radi[i,j] = star[10]
+   Temp[i,j] = star[11]
+N_iter = i
+N_part = j
+
+x_cm = np.array([np.average(X[i], weights=M[i]) for i in range(N_iter)])
+y_cm = np.array([np.average(Y[i], weights=M[i]) for i in range(N_iter)])
+z_cm = np.array([np.average(Z[i], weights=M[i]) for i in range(N_iter)])
+
+lim_min = np.min((X,Y))
+lim_max = np.max((X,Y))
+Alpha = np.clip((np.log(Lumi) - np.min(np.log(Lumi)))/(np.max(np.log(Lumi)) - np.min(np.log(Lumi))), 0.1, 1)
+
+# * Output choice
+print("\033[36m"+"Choose figures:"+"\033[0m")
+figure_1 = input("\033[32m"+"\tFigure 1 (energy, angular momentum, time): "+"\033[0m").upper() in YES
+figure_2 = input("\033[32m"+"\tFigure 2a and 2b (positions, speed): "+"\033[0m").upper() in YES
+figure_3 = input("\033[32m"+"\tFigure 3 (mass distribution): "+"\033[0m").upper() in YES
+figure_4 = input("\033[32m"+"\tFigure 4 (density distribution): "+"\033[0m").upper() in YES
+figure_5 = input("\033[32m"+"\tFigure 5 (HR diagram): "+"\033[0m").upper() in YES
+print("")
+
+print("\033[36m"+"General settings:"+"\033[0m")
+time =  np.int32(input("\033[32m"+"\tStatic images drawing instant ({} < int < {}): ".format(0,N_iter)+"\033[0m"))
+if input("\033[32m"+"\tWindow size restriction (y/n): "+"\033[0m").upper() in YES:
+      lim = np.float32(input("\033[32m"+"\t\tLimit [pc]: "+"\033[0m"))
+      lim_min = -lim/2
+      lim_max = +lim/2
+print("")
+
+draw_time = scatter_time = distrib_time = diagram_time = time
+
+if figure_1:
+   print("\033[36m"+"Figure 1 (energy, angular momentum, time):"+"\033[0m")
+   plot_E = input("\033[32m"+"\tPlot energy vs. time (y/n): "+"\033[0m").upper() in YES
+   plot_L = input("\033[32m"+"\tPlot angular momentum vs. time (y/n): "+"\033[0m").upper() in YES
+   plot_save = input("\033[32m"+"\tSave figure (y/n): "+"\033[0m").upper() in YES
+   print("")
+else: plot_E = plot_L = plot_save = False
+
+if figure_2:
+   print("\033[36m"+"Figure 2a and 2b (positions, speed):"+"\033[0m")
+   scatter_2D = input("\033[32m"+"\tScatter positions in 2D (y/n): "+"\033[0m").upper() in YES
+   scatter_3D = input("\033[32m"+"\tScatter positions in 3D (y/n): "+"\033[0m").upper() in YES
+   scatter_speed = input("\033[32m"+"\tShow speed vectors (y/n): "+"\033[0m").upper() in YES
+   scatter_save = input("\033[32m"+"\tSave figure (y/n): "+"\033[0m").upper() in YES
+   scatter_animate = input("\033[32m"+"\tAnimate over time (y/n): "+"\033[0m").upper() in YES
+   print("")
+else: scatter_2D = scatter_3D = scatter_speed = scatter_save = scatter_animate = False
+
+if figure_3:
+   print("\033[36m"+"Figure 3 (mass distribution):"+"\033[0m")
+   distrib_mass = input("\033[32m"+"\tShow mass distribution (y/n): "+"\033[0m").upper() in YES
+   distrib_bins = np.int32(input("\033[32m"+"\tNumber of bins (int): "+"\033[0m"))
+   distrib_save = input("\033[32m"+"\tSave figure (y/n): "+"\033[0m").upper() in YES
+   print("")
+else: distrib_mass = distrib_bins = distrib_save = False
+
+if figure_4:
+   print("\033[36m"+"Figure 4 (density distribution):"+"\033[0m")
+   draw_density = input("\033[32m"+"\tShow density distribution (y/n): "+"\033[0m").upper() in YES
+   draw_bins = np.float32(input("\033[32m"+"\tWidth of bins [pc]: "+"\033[0m"))
+   draw_max = np.float32(input("\033[32m"+"\tMaximum distance [pc]: "+"\033[0m"))
+   draw_save = input("\033[32m"+"\tSave figure (y/n): "+"\033[0m").upper() in YES
+   draw_animate = input("\033[32m"+"\tAnimate over time (y/n): "+"\033[0m").upper() in YES
+   print("")
+else: draw_density = draw_bins = draw_max = draw_save = draw_animate = False
+
+if figure_5:
+   print("\033[36m"+"Figure 5 (HR diagram):"+"\033[0m")
+   diagram_HT = input("\033[32m"+"\tDraw HR diagram (y/n): "+"\033[0m").upper() in YES
+   diagram_save = input("\033[32m"+"\tSave figure (y/n): "+"\033[0m").upper() in YES
+   diagram_animate = input("\033[32m"+"\tAnimate over time (y/n): "+"\033[0m").upper() in YES
+else : diagram_HT = diagram_save = diagram_animate = False
+
+# * Plotting
+if plot_E or plot_L:
+   fig1, axs1 = plt.subplots(2, sharex=True)
+   if plot_E and not plot_L: # Just the energy
+      axs1[1].set_xlabel("$t\\ [\\mathrm{{Myr}}]$")
+      axs1[0].set_ylabel("$E\\ [\\mathrm{{pc^2\\ M_\\odot\\ Myr^{{-2}}}}]$")
+      axs1[1].set_ylabel("$\\Delta E/E_0$")
+      axs1[0].plot(T, E, color="C3")
+      axs1[1].plot(T, DE, color="C3")
+   elif plot_L and not plot_E: # Just the angular momentum
+      axs1[1].set_xlabel("$t\\ [\\mathrm{{Myr}}]$")
+      axs1[0].set_ylabel("$L\\ [\\mathrm{{pc^2\\ M_\\odot\\ Myr^{{-1}}}}]$")
+      axs1[1].set_ylabel("$\\Delta L/L_0$")
+      axs1[0].plot(T, L, color="C4")
+      axs1[1].plot(T, DL, color="C4")
+   elif plot_L and plot_E: # ...both
+      xas1 = [axs1[i].twinx() for i in range(len(axs1))]
+      axs1[1].set_xlabel("$t\\ [\\mathrm{{Myr}}]$")
+      axs1[0].set_ylabel("$E\\ [\\mathrm{{pc^2\\ M_\\odot\\ Myr^{{-2}}}}]$")
+      axs1[1].set_ylabel("$\\Delta E/E_0$")
+      xas1[0].set_ylabel("$L\\ [\\mathrm{{pc^2\\ M_\\odot\\ Myr^{{-1}}}}]$")
+      xas1[1].set_ylabel("$\\Delta L/L_0$")
+      axs1[0].plot(T, E, color="C3")
+      axs1[1].plot(T, DE, color="C3")
+      xas1[0].plot(T, L, color="C4")
+      xas1[1].plot(T, DL, color="C4")
+      axs1[0].tick_params(axis='y', which="both", colors='C3')
+      axs1[0].yaxis.label.set_color('C3')
+      xas1[0].tick_params(axis='y', which="both", colors='C4')
+      xas1[0].yaxis.label.set_color('C4')
+      axs1[1].tick_params(axis='y', which="both", colors='C3')
+      axs1[1].yaxis.label.set_color('C3')
+      xas1[1].tick_params(axis='y', which="both", colors='C4')
+      xas1[1].yaxis.label.set_color('C4')
+   if plot_save:
+      fig1.tight_layout()
+      fig1.savefig(out_directory+date_prefix+"_figure_1")
+      print("\033[94m"+"Info: Figure saved as " + "\033[0m" + date_prefix+"_figure_1")
+
+if scatter_2D or scatter_3D:
+   if scatter_2D:
+      fig2a, ax2a = plt.subplots(1)
+      # ax2a.set_facecolor('k')
+      i = scatter_time
+      ax2a.set_title("$t =$ {:.2f} Myr".format(T[i]))
+      if scatter_speed: 
+         qvr2a = ax2a.quiver(X[i]-x_cm[i],Y[i]-y_cm[i],VX[i],VY[i], alpha=0.1, headwidth=2, headlength=2, headaxislength=2, color="k")
+         ax2a.quiverkey(qvr2a, 0.95, 0.05, 1, "Velocity scale: $1\mathrm{{~km\\cdot s^{{-1}}}}$", labelpos='W', coordinates='figure', color="k", alpha=1)
+      sct2a = ax2a.scatter(X[i]-x_cm[i],Y[i]-y_cm[i],s=np.clip(Radi[i]**2,0.5,30), c=Temp[i], alpha=Alpha[i], cmap="RdYlBu", vmin=Tmin, vmax=Tmax)
+      ax2a.set_xlabel("$x$ [pc]")
+      ax2a.set_ylabel("$y$ [pc]")
+      ax2a.set_xlim(lim_min, lim_max)
+      ax2a.set_ylim(lim_min, lim_max)
+      ax2a.set_aspect("equal")
+      fig2a.colorbar(sct2a, label="Stellar effective temperature $T$ [K]", extend='both')
+      fig2a.tight_layout()
+      if scatter_save:
+         fig2a.savefig(out_directory+date_prefix+"_figure_2a")
+         print("\033[94m"+"Info: Figure saved as " + "\033[0m" + date_prefix+"_figure_2a")
+   if scatter_3D:
+      fig2b, ax2b = plt.subplots(1, subplot_kw={"projection": "3d"})
+      # ax2b.set_facecolor('k')
+      i = scatter_time
+      ax2b.set_title("$t =$ {:.2f} Myr".format(T[i]))
+      if scatter_speed: 
+         qvr2b = ax2b.quiver(X[i]-x_cm[i],Y[i]-y_cm[i],Z[i]-z_cm[i],VX[i],VY[i],VZ[i], alpha=0.1, color="k")
+         # ax2b.quiverkey(qvr2b, 0.95, 0.05, 2, "Velocity scale: $2\mathrm{{~km\\cdot s^{{-1}}}}$", labelpos='W', coordinates='figure')
+      sct2b = ax2b.scatter(X[i]-x_cm[i],Y[i]-y_cm[i],Z[i]-z_cm[i],s=np.clip(Radi[i]**2,0.5,30), c=Temp[i], alpha=Alpha[i], cmap="RdYlBu", vmin=Tmin, vmax=Tmax)
+      ax2b.set_xlabel("$x$ [pc]")
+      ax2b.set_ylabel("$y$ [pc]")
+      ax2b.set_zlabel("$z$ [pc]")
+      ax2b.set_xlim(lim_min, lim_max)
+      ax2b.set_ylim(lim_min, lim_max)
+      ax2b.set_zlim(lim_min, lim_max)
+      ax2b.set_aspect("equal")
+      fig2b.colorbar(sct2b, label="Stellar effective temperature $T$ [K]", location="left", extend='both')
+      fig2b.tight_layout()
+      if scatter_save:
+         fig2b.savefig(out_directory+date_prefix+"_figure_2b")
+         print("\033[94m"+"Info: Figure saved as " + "\033[0m" + date_prefix+"_figure_2b")
+
+if distrib_mass:
+   fig3, ax3 = plt.subplots(1)
+   i = distrib_time
+   ax3.set_title("$t =$ {:.2f} Myr".format(T[i]))
+   nbins = distrib_bins
+   hist, lbin = np.histogram(M[i], nbins)
+   cbins = 0.5*(lbin[1:]+ lbin[: -1])
+   ax3.scatter(cbins, hist, s=5)
+   ax3.set_xscale ("log")
+   ax3.set_yscale ("log")
+   ax3.set_xlabel("Masses $M\\ [\\mathrm{{M_\\odot}}]$")
+   ax3.set_ylabel("Star distribution")
+   if distrib_save:
+      fig3.savefig(out_directory+date_prefix+"_figure_3")
+      print("\033[94m"+"Info: Figure saved as " + "\033[0m" + date_prefix+"_figure_3")
+
+if draw_density:
+   fig4, ax4 = plt.subplots(1)
+   i = draw_time
+   bin_width = draw_bins
+   bin_max = draw_max
+   bins = np.arange(0,bin_max,bin_width)
+   Rho = density_hist(X,Y,M,T,bins,i)
+   cbins = 0.5*(bins[1:]+ bins[: -1])
+   ax4.set_title("$t =$ {:.2f} Myr".format(T[i]))
+   ax4.set_xscale ("log")
+   ax4.set_yscale ("log")
+   ax4.set_xlabel("Distance from center $r$ [pc]")
+   ax4.set_ylabel("Stellar masses density $\\rho\\ [\\mathrm{{M_\\odot\\ pc^{{-3}}}}]$")
+   ax4.scatter(cbins, Rho, s=5)
+   fig4.tight_layout()
+   if draw_save:
+      fig4.savefig(out_directory+date_prefix+"_figure_4")
+      print("\033[94m"+"Info: Figure saved as " + "\033[0m" + date_prefix+"_figure_4")
+   
+if diagram_HT:
+   fig5, ax5 = plt.subplots(1)
+   ax5.invert_xaxis()
+   i = diagram_time
+   ax5.set_title("$t =$ {:.2f} Myr".format(T[i]))
+   ax5.set_xlabel("Temperature $T$ [K]")
+   ax5.set_ylabel("Luminosity $L\\ [\\mathrm{{L_\\odot}}]$")
+   ax5.set_xscale ("log")
+   ax5.set_yscale ("log")
+   ax5.scatter(Temp[i], Lumi[i], s=5, alpha=Alpha[i], c=Temp[i], cmap="RdYlBu", vmin=Tmin, vmax=Tmax)
+   fig5.tight_layout()
+   if diagram_save:
+      fig5.savefig(out_directory+date_prefix+"_figure_5")
+      print("\033[94m"+"Info: Figure saved as " + "\033[0m" + date_prefix+"_figure_5")
+
+if scatter_animate:
+   figA2a,axA2a = plt.subplots(1)
+   # ax2a.set_facecolor('k')
+   axA2a.set_title("$t =$ {:.2f} Myr".format(T[0]))
+   if scatter_speed:
+      qvr = axA2a.quiver(X[0]-x_cm[0],Y[0]-y_cm[0],VX[0],VY[0], alpha=0.5, headwidth=2, headlength=2, headaxislength=2, color="k")
+      axA2a.quiverkey(qvr, 0.95, 0.05, 1, "Velocity scale: $1\mathrm{{~km\\cdot s^{{-1}}}}$", labelpos='W', coordinates='figure', color="k", alpha=1)
+   sctA2a = axA2a.scatter(X[0]-x_cm[0],Y[0]-y_cm[0],s=np.clip(Radi[0]**2,0.5,30), c=Temp[0], alpha=Alpha[0], cmap="RdYlBu", vmin=Tmin, vmax=Tmax)
+   figA2a.colorbar(sctA2a, label="Stellar effective temperature $T$ [K]", extend='both')
+   axA2a.set_xlabel("$x$ [pc]")
+   axA2a.set_ylabel("$y$ [pc]")
+   axA2a.set_aspect("equal")
+   axA2a.set_xlim(lim_min, lim_max)
+   axA2a.set_ylim(lim_min, lim_max)
+   figA2a.tight_layout()
+   def anim2a(i):
+      if scatter_speed:
+         qvr.set_offsets(np.array([X[i]-x_cm[i],Y[i]-y_cm[i]]).T)
+         qvr.set_UVC(VX[i],VY[i])
+      axA2a.set_title("$t =$ {:.2f} Myr".format(T[i]))
+      sctA2a.set_offsets(np.array([X[i]-x_cm[i],Y[i]-y_cm[i]]).T)
+      sctA2a.set_array(Temp[i])
+      sctA2a.set_sizes(np.clip(Radi[i]**2,0.5,30))
+      sctA2a.set_alpha(Alpha[i])
+      return sctA2a,
+      
+   ani2a = animation.FuncAnimation(fig=figA2a, func=anim2a, frames=N_iter, interval=100)
+   if scatter_save:
+      ani2a.save(out_directory+date_prefix+"_animation_2a.png", dpi=300, fps=10)
+      ani2a.save(out_directory+date_prefix+"_animation_2a.pdf", dpi=300, fps=10)
+      print("\033[94m"+"Info: Animation saved as " + "\033[0m" + date_prefix+"_animation_2a")
+
+
+if draw_animate:
+   figA4, axA4 = plt.subplots(1)
+   bin_width = draw_bins
+   bin_max = draw_max
+   bins = np.arange(0,bin_max,bin_width)
+   cbins = 0.5*(bins[1:]+ bins[: -1])
+   Rho = np.array([density_hist(X,Y,M,T,bins,k) for k in range(N_iter)])
+   axA4.set_title("$t =$ {:.2f} Myr".format(T[0]))
+   sctA4 = axA4.scatter(cbins, Rho[0], s=5)
+   axA4.set_xscale ("log")
+   axA4.set_yscale ("log")
+   axA4.set_ylim(top=np.max(Rho))
+   axA4.set_xlabel("Distance from center $r$ [pc]")
+   axA4.set_ylabel("Stellar masses density $\\rho\\ [\\mathrm{{M_\\odot\\ pc^{{-3}}}}]$")
+   figA4.tight_layout()
+   def anim4(i):
+      axA4.set_title("$t =$ {:.2f} Myr".format(T[i]))
+      sctA4.set_offsets(np.array([cbins, Rho[i]]).T)
+      return sctA4
+   ani4 = animation.FuncAnimation(fig=figA4, func=anim4, frames=N_iter, interval=100)
+   if scatter_save:
+      ani4.save(out_directory+date_prefix+"_animation_4.png", dpi=300, fps=10)
+      ani4.save(out_directory+date_prefix+"_animation_4.pdf", dpi=300, fps=10)
+      print("\033[94m"+"Info: Animation saved as " + "\033[0m" + date_prefix+"_animation_4")
+
+if diagram_animate:
+   figA5, axA5 = plt.subplots(1)
+   axA5.invert_xaxis()
+   axA5.set_title("$t =$ {:.2f} Myr".format(T[0]))
+   axA5.set_xlabel("Temperature $T$ [K]")
+   axA5.set_ylabel("Luminosity $L\\ [\\mathrm{{L_\\odot}}]$")
+   axA5.set_xscale ("log")
+   axA5.set_yscale ("log")
+   axA5.set_xlim(np.max(Temp), np.min(Temp))
+   axA5.set_ylim(np.min(Lumi), np.max(Lumi))
+   sctA5 = axA5.scatter(Temp[0], Lumi[0], s=5, alpha=Alpha[0], c=Temp[0], cmap="RdYlBu", vmin=Tmin, vmax=Tmax)
+   fig5.tight_layout()
+   def anim5(i):
+      axA5.set_title("$t =$ {:.2f} Myr".format(T[i]))
+      sctA5.set_offsets(np.array([Temp[i], Lumi[i]]).T)
+      sctA5.set_array(Temp[i])
+   ani5 = animation.FuncAnimation(fig=figA5, func=anim5, frames=N_iter, interval=100)
+   if diagram_save:
+      ani5.save(out_directory+date_prefix+"_animation_5.png", dpi=300, fps=10)
+      ani5.save(out_directory+date_prefix+"_animation_5.pdf", dpi=300, fps=10)
+      print("\033[94m"+"Info: Animation saved as " + "\033[0m" + date_prefix+"_animation_5")
+
+# * SHOW !
+plt.show(block=False)
+input("(press enter to quit)")
diff --git a/COSMIC/old/COSMIC_v2.py b/COSMIC/old/COSMIC_v2.py
new file mode 100644
index 0000000..e1a83d3
--- /dev/null
+++ b/COSMIC/old/COSMIC_v2.py
@@ -0,0 +1,164 @@
+"""
+* COSMIC : Cluster Orbital SysteM Integration Code
+* Version 2 - 2024-02-27
+@ Yaël Moussouni
+@ Unistra, P&E, MSc1-MoFP
+@ Observatory of Strasbourg (Intern)
+"""
+import numpy as np
+import amuse.units.units as u
+import datetime
+
+
+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.datamodel.particles import Channels
+from amuse.community.ph4.interface import Ph4
+from amuse.community.seba.interface import SeBa
+from amuse.couple.bridge import Bridge
+
+
+# * Parameters
+N_part = 10  # Number of particle
+N_iter = 30  # Number of iteration
+dt = 1|u.Myr   # Time increment
+
+R0 = 1|u.parsec       # Radius of the system
+W0 = 9                # King parameter
+z = 0.02                 # Metallicity
+
+# * Variables
+date = datetime.datetime.now()
+date = "{:04d}-{:02d}-{:02d}_{:02d}h{:02d}".format(date.year, date.month, date.day, date.hour, date.minute)
+
+directory = "./output/"
+filename_cluster = "{}_cluster.csv".format(date)
+filename_stars = "{}_stars.csv".format(date)
+
+# * Function
+def get_total_momentum(M, X, V): 
+   """
+   Compute the total momentum of particles with mass M, positions X and velocities V.
+   """
+   l = 0
+   for j in range(len(X)):
+      m = M[j].value_in(u.MSun)
+      x = X[j].value_in(u.pc)
+      v = V[j].value_in(u.pc/u.Myr)
+      l += m*np.cross(x,v)
+   return np.sqrt(np.sum(l**2))|(u.pc**2 * u.MSun / u.Myr)
+
+# * Initialization
+def create_cluster(N:int, R=1|u.pc, W=9):
+   """
+   Create cluster of N stars of radius R (default: 1 pc) following King model of parameter W (default: 9).
+   """
+   print("\033[94m"+"Info: Simulation with "+"\033[0m"+"{}".format(N)+"\033[94m"+" stars"+"\033[0m")
+   M_stars = new_salpeter_mass_distribution(N)
+   converter = nbody_system.nbody_to_si(M_stars.sum(), R)
+   stars = new_king_model(N, W, converter)
+   stars.mass = M_stars
+   stars.zams_mass = M_stars
+   return stars, converter
+
+# * Evolution
+def stellar_gravity_evolution(stars, converter, t_max=16|u.Myr, dt=1|u.Myr, out_dt=1|u.Myr):
+   """
+   Evolve a model with both gravity and stellar evolution. dt is the simulation time step, while out_dt is the delay between printing in file.
+   """
+   print("\033[94m"+"Info: Simulation duration of "+"\033[0m"+"{:.2f} Myr".format(t_max.value_in(u.Myr)))
+   init_files()
+
+   gravity = Ph4(converter) # gravity model
+   stellar = SeBa()         # stellar model
+   bridge = Bridge()        # bridge between models
+
+   # TODO gravity.initialize_code()
+
+   gravity.particles.add_particles(stars)
+   stellar.particles.add_particles(stars)
+
+   stellar.set_metallicity(z)
+   stellar.commit_particles()
+   stellar.update_particles(stars) # TODO pourquoi ?
+   
+   bridge.add_system(gravity)
+   bridge.add_system(stellar)
+
+   bridge.synchronize_model()
+
+   bridge.channels.add_channel(stellar.particles.new_channel_to(gravity.particles, attributes=["mass", "radius"])) 
+   # Sync the mass and radius between the two models
+
+   out_channel = Channels()
+   out_channel.add_channel(stellar.particles.new_channel_to(stars))
+   out_channel.add_channel(gravity.particles.new_channel_to(stars))
+
+   bridge.timestep = dt # TODO timestep dynamique
+   print(gravity.get_initial_timestep_median()) # TODO ????
+   i = 0
+   t = 0 | u.Myr
+   E0 = gravity.get_total_energy()
+   L0 = get_total_momentum(stars.mass, stars.position, stars.velocity)
+   while t <= t_max:
+      print("\033[36m"+"t: "+"\033[0m"+"{:02.2f} Myr/{:02.2f} Myr".format(t.value_in(u.Myr),t_max.value_in(u.Myr)))
+      bridge.evolve_model(t)
+      out_channel.copy()
+      E = gravity.get_total_energy()
+      L = get_total_momentum(stars.mass, stars.position, stars.velocity)
+      write_files(stars, i, t, E, E0, L, L0)
+      i += 1
+      t += out_dt
+   
+   print("\033[36m"+"t: "+"\033[0m"+"{:02.2f} Myr/{:02.2f} Myr".format(t_max.value_in(u.Myr),t_max.value_in(u.Myr)))
+   print("\033[36m"+"Simulation completed"+"\033[0m")
+   return 0
+
+# * Files initialization and write 
+def init_files():
+   file = open(directory+filename_cluster, "w")
+   file.write("#i (time), t [Myr], E [pc2 MSun Myr-2], (E-E0)/E0, L [pc2 MSun Myr-1], (L-L0)/L0\n")
+   file.close()
+
+   file = open(directory+filename_stars, "w")
+   file.write("#i (time), j (star), m [MSun], x [pc], y [pc], z [pc], vx [km s-1], vy [km s-1], vz [km s-1], L [LSun], R [RSun], T [K]\n")
+   file.close()
+
+   print("\033[94m"+"Info: Created two output files with date " + "\033[0m" + date.replace("h", ":").replace("_", " "))
+   return 0
+
+def write_files(stars, i, t, E, E0, L, L0):
+   file = open(directory+filename_cluster, "a")
+   file.write("{:d}, ".format(i))
+   file.write("{:+e}, ".format(t.value_in(u.Myr)))
+   file.write("{:+e}, ".format(E.value_in(u.pc**2 * u.MSun * u.Myr**-2)))
+   file.write("{:+e}, ".format((E-E0)/E0))
+   file.write("{:+e}, ".format(L.value_in(u.pc**2 * u.MSun * u.Myr**-1)))
+   file.write("{:+e}\n".format((L-L0)/L0))
+   file.close()
+
+   file = open(directory+filename_stars, "a")
+   for j in range(N_part):
+      file.write("{:d}, ".format(i))
+      file.write("{:d}, ".format(j))
+      file.write("{:+e}, ".format(stars.mass[j].value_in(u.MSun)))
+      file.write("{:+e}, ".format(stars.position[j,0].value_in(u.pc)))
+      file.write("{:+e}, ".format(stars.position[j,1].value_in(u.pc)))
+      file.write("{:+e}, ".format(stars.position[j,2].value_in(u.pc)))
+      file.write("{:+e}, ".format(stars.velocity[j,0].value_in(u.km * u.s**-1)))
+      file.write("{:+e}, ".format(stars.velocity[j,1].value_in(u.km * u.s**-1)))
+      file.write("{:+e}, ".format(stars.velocity[j,2].value_in(u.km * u.s**-1)))
+      file.write("{:+e}, ".format(stars.luminosity[j].value_in(u.LSun)))
+      file.write("{:+e}, ".format(stars.radius[j].value_in(u.RSun)))
+      file.write("{:+e}\n".format(stars.temperature[j].value_in(u.K)))
+   file.close()
+   return 0
+
+# * Main
+cluster, nbody_converter = create_cluster(N_part, R0, W0)
+stellar_gravity_evolution(cluster, nbody_converter, N_iter*dt, dt, dt)
+
+
+
+
diff --git a/COSMIC/source/COSMIC_v3.py b/COSMIC/source/COSMIC_v3.py
new file mode 100644
index 0000000..56b1dd9
--- /dev/null
+++ b/COSMIC/source/COSMIC_v3.py
@@ -0,0 +1,28 @@
+"""
+* 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)
+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, 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
new file mode 100644
index 0000000..9519bc8
--- /dev/null
+++ b/COSMIC/source/COSMIC_v3_config.py
@@ -0,0 +1,96 @@
+"""
+* 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", "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["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 == "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/source/COSMIC_v3_emission.py b/COSMIC/source/COSMIC_v3_emission.py
new file mode 100644
index 0000000..76039db
--- /dev/null
+++ b/COSMIC/source/COSMIC_v3_emission.py
@@ -0,0 +1,27 @@
+"""
+* 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
+    return stars
diff --git a/COSMIC/source/COSMIC_v3_evolve.py b/COSMIC/source/COSMIC_v3_evolve.py
new file mode 100644
index 0000000..a38e148
--- /dev/null
+++ b/COSMIC/source/COSMIC_v3_evolve.py
@@ -0,0 +1,44 @@
+"""
+* 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
+
+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()
+        if compute_X_emission: 
+            stars = emission.rotation(stars, params)
+            stars = emission.X_emission(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/source/COSMIC_v3_galaxy.py b/COSMIC/source/COSMIC_v3_galaxy.py
new file mode 100644
index 0000000..e69de29
diff --git a/COSMIC/source/COSMIC_v3_init.py b/COSMIC/source/COSMIC_v3_init.py
new file mode 100644
index 0000000..ae7f33c
--- /dev/null
+++ b/COSMIC/source/COSMIC_v3_init.py
@@ -0,0 +1,193 @@
+"""
+* 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
+
+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.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
+    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()
+    vc = galaxy_model.vel_circ(stars.center_of_mass().length())
+    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)
+    stars.vx += - vc * np.sin(phi)
+    stars.vy += vc * np.cos(phi)
+    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/source/COSMIC_v3_output.py b/COSMIC/source/COSMIC_v3_output.py
new file mode 100644
index 0000000..0841c5f
--- /dev/null
+++ b/COSMIC/source/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