diff --git a/COSMIC/cexe.sh b/COSMIC/cexe.sh new file mode 100644 index 0000000..9551502 --- /dev/null +++ b/COSMIC/cexe.sh @@ -0,0 +1,6 @@ +if [[ "$VIRTUAL_ENV" == "/Users/yael/.venv/amuse" ]] +then + python source/CEXE_v3.py +else + echo "please activate amuse" +fi diff --git a/COSMIC/input/default.cfg b/COSMIC/input/default.cfg index a38ba77..15fedd6 100644 --- a/COSMIC/input/default.cfg +++ b/COSMIC/input/default.cfg @@ -1,4 +1,4 @@ -N_stars 10 +N_stars 128 timestep 1 output_times 1 2 5 10 20 50 100 200 500 1000 2000 5000 R0_king 1 @@ -18,6 +18,6 @@ cluster_position_z 0 out_directory ./output/ in_directory ./input/ config default.cfg -workers_stellar 1 -workers_gravity 1 +workers_stellar 1 +workers_gravity 3 format_type csv diff --git a/COSMIC/source/CEXE_v3.py b/COSMIC/source/CEXE_v3.py new file mode 100644 index 0000000..e883b94 --- /dev/null +++ b/COSMIC/source/CEXE_v3.py @@ -0,0 +1,143 @@ +""" +* CEXE: COSMIC Extension for X-ray Emission +* Version 3 +@ Yaël Moussouni +@ Unistra, P&E, MSc1-MoFP +@ Observatory of Strasbourg (Intern) +""" + +import numpy as np +import os +import subprocess as sp +from astropy import units as u + +directory = "./output/" + +# SIMPUT (sources) +val_Elow = 0.1 +val_Eup = 15 +val_NBins = 1000 +val_logEgrid = "yes" +val_Emin = 2 +val_Emax = 10 +val_nH = 0.2 +val_flux_min = 2 +val_flux_max = 10 +# SIXTE (instrument) +val_ra = 0 +val_dec = 0 +val_Exposure = 1000 +XML_file = "TODO" # TODO !! + +# * COSMIC +def open_file_terminal(): + filelist = np.sort([f for f in os.listdir(directory) if "stars" in f]) + print("\033[36m"+"Available files:"+"\033[0m") + for i in range(len(filelist)): + print("\033[36m"+"\t{}: ".format(i+1) + "\033[0m" + filelist[i]) + j = int(input("\033[32m"+"Select a file number: "+"\033[0m"))-1 + filename = filelist[j] + date = filename[filename.find("COSMIC_")+len("COSMIC_") : filename.find("_stars")] + time = filename[filename.find("stars_")+len("stars_") : filename.find("_Myr")] + print("") + file = np.loadtxt(directory+filename, delimiter=",", comments="#") + with open(directory+filename) as f: head = f.readline().strip('\n')[1:] + return file, head, date, time + +# * XSpec +def xspec_build(date, time): + with open(directory+"CEXE_XSPEC_{}_{}_Myr_0.6_0.0.sh".format(date, time), "w+") as xspec_file_0600: + xspec_file_0600.write("model phabs*TODO\n") + xspec_file_0600.write("phabs:nH>{}\n".format(val_nH)) + xspec_file_0600.write("TODO:...\n") + xspec_file_0600.write("flux {} {}\n".format(val_flux_min, val_flux_max)) + xspec_file_0600.write("save model {}XSPEC_{}_{}_Myr_0.6_0.0.xcm\n".format(directory, date, time)) + xspec_file_0600.write("quit") + with open(directory+"CEXE_XSPEC_{}_{}_Myr_0.5_0.0.sh".format(date, time), "w+") as xspec_file_0500: + xspec_file_0500.write("model phabs*TODO\n") + xspec_file_0500.write("phabs:nH>{}\n".format(val_nH)) + xspec_file_0500.write("TODO:...\n") + xspec_file_0500.write("flux {} {}\n".format(val_flux_min, val_flux_max)) + xspec_file_0500.write("save model {}XSPEC_{}_{}_Myr_0.5_0.0.xcm\n".format(directory, date, time)) + xspec_file_0500.write("quit") + with open(directory+"CEXE_XSPEC_{}_{}_Myr_0.4_1.0.sh".format(date, time), "w+") as xspec_file_0410: + xspec_file_0410.write("model phabs*TODO\n") + xspec_file_0410.write("phabs:nH>{}\n".format(val_nH)) + xspec_file_0410.write("TODO:...\n") + xspec_file_0410.write("flux {} {}\n".format(val_flux_min, val_flux_max)) + xspec_file_0410.write("save model {}XSPEC_{}_{}_Myr_0.4_1.0.xcm\n".format(directory, date, time)) + xspec_file_0410.write("quit") + with open(directory+"CEXE_XSPEC_{}_{}_Myr_0.2_0.8.sh".format(date, time), "w+") as xspec_file_0208: + xspec_file_0208.write("model phabs*TODO\n") + xspec_file_0208.write("phabs:nH>{}\n".format(val_nH)) + xspec_file_0208.write("TODO:...\n") + xspec_file_0208.write("flux {} {}\n".format(val_flux_min, val_flux_max)) + xspec_file_0208.write("save model {}XSPEC_{}_{}_Myr_0.2_0.8.xcm\n".format(directory, date, time)) + xspec_file_0208.write("quit") + return None + +def xspec_exec(date, time): + sp.run("echo \"xspec < {}CEXE_XSPEC_{}_{}_Myr_0.6_0.0.sh - TO EXECUTE\"".format(directory, date, time), shell=True) + sp.run("echo \"xspec < {}CEXE_XSPEC_{}_{}_Myr_0.5_0.0.sh - TO EXECUTE\"".format(directory, date, time), shell=True) + sp.run("echo \"xspec < {}CEXE_XSPEC_{}_{}_Myr_0.4_1.0.sh - TO EXECUTE\"".format(directory, date, time), shell=True) + sp.run("echo \"xspec < {}CEXE_XSPEC_{}_{}_Myr_0.2_0.8.sh - TO EXECUTE\"".format(directory, date, time), shell=True) + return None + +# * SIMPUT +def simput_build(file, col, date, time): + with open(directory+"CEXE_SIMPUT_{}_{}_Myr.sh".format(date, time), "w+") as simput_file: + Infiles = "" + for i in range(len(file)): + star_flux = file[i,col["X_luminosity"]] / (4*np.pi*(file[i,col["dist"]] * u.pc).to(u.cm).value**2) + star_ra = file[i,col["ra"]] + star_dec = file[i,col["dec"]] + star_T_0 = file[i,col["X_temperature_0"]] + star_T_1 = file[i,col["X_temperature_1"]] + + Infiles += "{}SIMPUT_{}_{}_Myr_{}.fits,".format(directory, date, time, i) + + simput_file.write("$SIXTE/bin/simputfile Simput={}SIMPUT_{}_{}_Myr_{}.fits ".format(directory, date, time, i)) + simput_file.write("Src_Name = star_{} ".format(i)) + simput_file.write("RA = {} Dec = {} ".format(star_ra, star_dec)) + simput_file.write("srcFlux = {} ".format(star_flux)) + simput_file.write("Elow={} Eup={} NBins={} logEgrid={} Emin={} Emax={} ".format(val_Elow, val_Eup, val_NBins, val_logEgrid, val_Emin, val_Emax)) + + simput_file.write("XSPECFile={}XSPEC_{}_{}_Myr_{:2.01f}_{:2.01f}.xcm".format(directory, date, time, star_T_0, star_T_1)) + simput_file.write("\n") + simput_file.write("$SIXTE/bin/simputmerge ") + simput_file.write("Infiles={} ".format(Infiles[:-1])) + simput_file.write("Outfile={}SIMPUT_{}_{}_Myr.fits".format(directory, date, time)) + simput_file.write("\n") + return 0 + +def simput_exec(date, time): + sp.run(["chmod", "a+x", "{}CEXE_SIMPUT_{}_{}_Myr.sh".format(directory, date, time)]) + sp.run(["echo" ,"{}CEXE_SIMPUT_{}_{}_Myr.sh - TO EXECUTE".format(directory, date, time)]) + return 0 + +def sixte_build(date, time): + with open(directory+"CEXE_SIXTE_{}_{}_Myr.sh".format(date, time), "w+") as sixte_file: + sixte_file.write("$SIXTE/bin/runsixt ") + sixte_file.write("XMLFile = {} ".format(XML_file)) + sixte_file.write("RA = {} DEC = {} ".format(val_ra, val_dec)) + # sixte_file.write("Prefix = SIXTE_ ") + sixte_file.write("Simput = {}SIMPUT_{}_{}_Myr.fits ".format(directory, date, time)) + sixte_file.write("EvtFile = {}SIXTE_Event_{}_{}_Myr.fits ".format(directory, date, time)) + sixte_file.write("Exposure = val_Exposure") + sixte_file.write("\n") + return 0 + +def sixte_exec(date, time): + sp.run(["chmod", "a+x", "{}CEXE_SIXTE_{}_{}_Myr.sh".format(directory, date, time)]) + sp.run(["echo", "{}CEXE_SIXTE_{}_{}_Myr.sh - TO EXECUTE".format(directory, date, time)]) + return 0 + +file, head, date, time = open_file_terminal() +col = {head.split(",")[i]:i for i in range(len(head.split(",")))} + +xspec_build(date, time) +xspec_exec(date, time) +simput_build(file, col, date, time) +simput_exec(date, time) +sixte_build(date, time) +sixte_exec(date, time) diff --git a/COSMIC/source/COSMIC_v3.py b/COSMIC/source/COSMIC_v3.py index e6733d4..659dfb1 100644 --- a/COSMIC/source/COSMIC_v3.py +++ b/COSMIC/source/COSMIC_v3.py @@ -1,5 +1,5 @@ """ -* COSMIC : Cluster Orbital SysteM Integration Code +* COSMIC: Cluster Orbital SysteM Integration Code * Version 3 @ Yaël Moussouni @ Unistra, P&E, MSc1-MoFP diff --git a/COSMIC/source/COSMIC_v3_emission.py b/COSMIC/source/COSMIC_v3_emission.py index 6a59ebe..ce7ca95 100644 --- a/COSMIC/source/COSMIC_v3_emission.py +++ b/COSMIC/source/COSMIC_v3_emission.py @@ -10,8 +10,9 @@ 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) + 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) # TODO ROTATION TABLES p = np.random.uniform(0, 1, N) stars.X_vsini = stars.X_v * np.sqrt(2*p-p**2) return stars @@ -19,9 +20,12 @@ def rotation(stars, params): 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)) + 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) + 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 diff --git a/COSMIC/source/COSMIC_v3_galaxy.py b/COSMIC/source/COSMIC_v3_galaxy.py index b646e44..fb190a3 100644 --- a/COSMIC/source/COSMIC_v3_galaxy.py +++ b/COSMIC/source/COSMIC_v3_galaxy.py @@ -9,7 +9,9 @@ import numpy as np import matplotlib.pyplot as plt import amuse.units.units as u +from amuse.ext.galactic_potentials import MWpotentialBovy2015 +MilkyWay_Bovy2015 = MWpotentialBovy2015 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 diff --git a/COSMIC/source/COSMIC_v3_init.py b/COSMIC/source/COSMIC_v3_init.py index 768b9c3..f699b42 100644 --- a/COSMIC/source/COSMIC_v3_init.py +++ b/COSMIC/source/COSMIC_v3_init.py @@ -139,11 +139,11 @@ def create_channels(): return channels def generate_galaxy(stars, params): - galaxy_model = galaxy.MilkyWay_AMUSE() + galaxy_model = galaxy.MilkyWay_Bovy2015() 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)) + vc = galaxy_model.circular_velocity(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) diff --git a/COSMIC/source/__pycache__/COSMIC_v3.cpython-39.pyc b/COSMIC/source/__pycache__/COSMIC_v3.cpython-39.pyc new file mode 100644 index 0000000..3795d63 Binary files /dev/null and b/COSMIC/source/__pycache__/COSMIC_v3.cpython-39.pyc differ diff --git a/COSMIC/source/__pycache__/COSMIC_v3_config.cpython-39.pyc b/COSMIC/source/__pycache__/COSMIC_v3_config.cpython-39.pyc new file mode 100644 index 0000000..b5fd384 Binary files /dev/null and b/COSMIC/source/__pycache__/COSMIC_v3_config.cpython-39.pyc differ diff --git a/COSMIC/source/__pycache__/COSMIC_v3_coordinates.cpython-39.pyc b/COSMIC/source/__pycache__/COSMIC_v3_coordinates.cpython-39.pyc new file mode 100644 index 0000000..f6b0ec3 Binary files /dev/null and b/COSMIC/source/__pycache__/COSMIC_v3_coordinates.cpython-39.pyc differ diff --git a/COSMIC/source/__pycache__/COSMIC_v3_emission.cpython-39.pyc b/COSMIC/source/__pycache__/COSMIC_v3_emission.cpython-39.pyc new file mode 100644 index 0000000..ca93ec9 Binary files /dev/null and b/COSMIC/source/__pycache__/COSMIC_v3_emission.cpython-39.pyc differ diff --git a/COSMIC/source/__pycache__/COSMIC_v3_evolve.cpython-39.pyc b/COSMIC/source/__pycache__/COSMIC_v3_evolve.cpython-39.pyc new file mode 100644 index 0000000..da24154 Binary files /dev/null and b/COSMIC/source/__pycache__/COSMIC_v3_evolve.cpython-39.pyc differ diff --git a/COSMIC/source/__pycache__/COSMIC_v3_galaxy.cpython-39.pyc b/COSMIC/source/__pycache__/COSMIC_v3_galaxy.cpython-39.pyc new file mode 100644 index 0000000..7bfe617 Binary files /dev/null and b/COSMIC/source/__pycache__/COSMIC_v3_galaxy.cpython-39.pyc differ diff --git a/COSMIC/source/__pycache__/COSMIC_v3_init.cpython-39.pyc b/COSMIC/source/__pycache__/COSMIC_v3_init.cpython-39.pyc new file mode 100644 index 0000000..d1ce3ec Binary files /dev/null and b/COSMIC/source/__pycache__/COSMIC_v3_init.cpython-39.pyc differ diff --git a/COSMIC/source/__pycache__/COSMIC_v3_output.cpython-39.pyc b/COSMIC/source/__pycache__/COSMIC_v3_output.cpython-39.pyc new file mode 100644 index 0000000..ccbf44c Binary files /dev/null and b/COSMIC/source/__pycache__/COSMIC_v3_output.cpython-39.pyc differ