Add files via upload

HR_Maker_v2-2.py

@@ -0,0 +1,444 @@
+"""
+* HR Maker: Building HR diagram from BVR.csv files and fit the main sequence
+* Version 2.2 - 2024-04-16
+@ Yaël Moussouni
+@ Unistra, P&E, MSc1-MoFP
+@ Observatory of Strasbourg (Intern)
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+import os
+from scipy.interpolate import LinearNDInterpolator
+from amuse.community.seba.interface import SeBa
+from amuse.datamodel import Particles
+from amuse.units import units as u
+from astropy import units as u2
+import scipy.constants as const
+from astropy.io import fits
+from astropy.coordinates import SkyCoord
+from dustmaps.bayestar import BayestarQuery
+
+# * Parameters
+
+if "YII" in plt.style.available: plt.style.use("YII")
+
+# * Variables
+in_directory = "./Cluster_files/"
+fig_directory = "../Figs/HR_Maker_figs/"
+
+Cluster_prefix = "Cluster_"
+MS_prefix = "MS_"
+BVR_suffix = "_BVR.csv"
+
+metallicity = 0
+ZSun = 0.0134
+bessel_directory = "YBC_tables/ubvrijhklm/regrid/" # Bessel
+atlas9_prefix = "Avodonnell94Rv3.1f" # atlas9
+atlas9_suffix = "k2odfnew.BC.fits" # atlas9
+
+phoenix_filename = "Avodonnell94Rv3.1BT-Settl_M-0.0_a+0.0.BC.fits" # phoenix
+atlas12_filename = "Avodonnell94Rv3.147TucP1.BC.fits" # atlas12
+
+if metallicity >= 0: atlas9_filename = atlas9_prefix + "p{:02d}".format(metallicity) + atlas9_suffix
+else: atlas9_filename = atlas9_prefix + "m{:02d}".format(metallicity) + atlas9_suffix
+
+M_bol_sun = 4.74 # Solar bolometric magnitude
+sigma = const.Stefan_Boltzmann | u.W * u.m**-2 * u.K**-4 # Stefan-Boltzmann constant
+
+YES = ["Y", "YES", "1", "OUI"]
+
+stellar = SeBa()
+t_min = 0 # Myr
+t_max = 2e5 # Myr
+N_t = int(t_max//10+1)
+N_color = 1000
+N_stars = 500
+
+print("\033[94m"+"Info: loading extinction map..."+"\033[0m")
+ext_query = BayestarQuery(version='bayestar2019')
+print("\033[94m"+"Info: done!"+"\033[0m")
+
+coeff_ext_BV = 3.626-2.742
+coeff_ext_V = 2.742
+R_ext = 3.1
+
+# * Observational data import
+Cluster_filelist = np.sort(os.listdir(in_directory))
+Cluster_filelist = [f for f in Cluster_filelist if Cluster_prefix in f]
+MS_filelist = np.sort(os.listdir(in_directory))
+MS_filelist = [f for f in MS_filelist if MS_prefix in f]
+Cluster_sources = []
+MS_sources = []
+sources = []
+for i in range(len(Cluster_filelist)):
+   k = Cluster_filelist[i][len(Cluster_prefix):].find("_BVR")
+   source = Cluster_filelist[i][len(Cluster_prefix):k+len(Cluster_prefix)]
+   if source not in Cluster_sources: Cluster_sources.append(source)
+for i in range(len(MS_filelist)):
+   k = MS_filelist[i][len(MS_prefix):].find("_BVR")
+   source = MS_filelist[i][len(MS_prefix):k+len(MS_prefix)]
+   if source not in MS_sources: MS_sources.append(source)
+
+sources = np.intersect1d(Cluster_sources, MS_sources)
+
+print("\033[36m"+"Available source:"+"\033[0m")
+for i in range(len(sources)):
+   print("\033[36m"+"\t{}: ".format(i+1) + "\033[0m" + "{}".format(sources[i]))
+k = int(input("\033[32m"+"Select source number: "+"\033[0m"))-1
+print("")
+source = sources[k]
+
+BVR_Cluster = np.loadtxt(in_directory+Cluster_prefix+source+BVR_suffix, skiprows=1, usecols=(0,1,2,3,4,5,6,7), comments="#", delimiter=",")
+BVR_MS = np.loadtxt(in_directory+MS_prefix+source+BVR_suffix, skiprows=1, usecols=(0,1,2,3,4,5,6,7), comments="#", delimiter=",")
+
+RA_Cluster = BVR_Cluster[:,0]
+DE_Cluster = BVR_Cluster[:,1]
+B_Cluster = BVR_Cluster[:,2]
+V_Cluster = BVR_Cluster[:,3]
+R_Cluster = BVR_Cluster[:,4]
+dB_Cluster = BVR_Cluster[:,5]
+dV_Cluster = BVR_Cluster[:,6]
+dR_Cluster = BVR_Cluster[:,7]
+
+B_MS = BVR_MS[:,2]
+V_MS = BVR_MS[:,3]
+R_MS = BVR_MS[:,4]
+dB_MS = BVR_MS[:,5]
+dV_MS = BVR_MS[:,6]
+dR_MS = BVR_MS[:,7]
+
+N_MS = len(B_MS)
+N_Cluster = len(B_Cluster)
+
+ra = np.mean(RA_Cluster)
+de = np.mean(DE_Cluster)
+
+Coords_Cluster = SkyCoord(ra*u2.deg, de*u2.deg, distance=0*u2.pc, frame='icrs')
+print("\033[36m"+"Cluster coordinates: "+"\033[0m" + Coords_Cluster.to_string('hmsdms'))
+# * Interpolation data import
+
+hdul = fits.open(bessel_directory+atlas9_filename)
+# hdul = fits.open(bessel_directory+phoenix_filename)
+# hdul = fits.open(bessel_directory+atlas12_filename)
+interp_data = hdul[1].data
+
+# * Main sequence generation
+print("\033[94m"+"Info: generating and evolving stars..."+"\033[0m")
+
+
+def stars_MS(M_min = .5, M_max = 5, z=ZSun, N = 100):
+   M_stars = np.linspace(M_min, M_max, N) | u.MSun
+   stars = Particles(len(M_stars), mass=M_stars)
+   stellar = SeBa()
+   stellar.particles.add_particles(stars)
+   stellar.set_metallicity(z)
+   stellar.commit_particles()
+   stellar.update_particles(stars)
+   return stars
+
+# * Interpolation functions
+
+interp_logTeff = interp_data.field(0)
+interp_logg = interp_data.field(1)
+BC_B_bessel = interp_data.field(3)
+BC_V_bessel = interp_data.field(4)
+BC_R_bessel = interp_data.field(5)
+
+"""
+Interpolation using YBC tables (Yang, 2019)
+@ https://ui.adsabs.harvard.edu/abs/2019A%26A...632A.105C/abstract
+@ https://sec.center/YBC/
+"""
+interp_BC_B = LinearNDInterpolator(list(zip(interp_logTeff, interp_logg)), BC_B_bessel)
+interp_BC_V = LinearNDInterpolator(list(zip(interp_logTeff, interp_logg)), BC_V_bessel) 
+interp_BC_R = LinearNDInterpolator(list(zip(interp_logTeff, interp_logg)), BC_R_bessel)
+
+# * Interpolation quantities computation
+
+def extract_BVR_stars(stars):
+   L_stars = stars.luminosity.value_in(u.LSun)
+   Teff_stars = ((stars.luminosity / (4*np.pi * stars.radius**2 * sigma))**(1/4)).value_in(u.K)
+   G_stars = (u.constants.G * stars.mass / stars.radius**2).value_in(u.m * u.s**(-2))
+   # type_stars = stars.stellar_type.value_in(u.stellar_type)
+   BC_B_stars = interp_BC_B(np.log10(Teff_stars),np.log10(G_stars))
+   BC_V_stars = interp_BC_V(np.log10(Teff_stars),np.log10(G_stars))
+   BC_R_stars = interp_BC_R(np.log10(Teff_stars),np.log10(G_stars))
+
+   M_bol_stars = M_bol_sun - 2.5*np.log10(L_stars)
+
+   B_abs_stars = M_bol_stars-BC_B_stars
+   V_abs_stars = M_bol_stars-BC_V_stars
+   R_abs_stars = M_bol_stars-BC_R_stars
+
+   # w1_not = np.where(type_stars != 1)
+
+   # B_abs_stars[w1_not] = np.nan
+   # V_abs_stars[w1_not] = np.nan
+   # R_abs_stars[w1_not] = np.nan
+
+   return B_abs_stars, V_abs_stars, R_abs_stars
+
+
+stars = stars_MS(N = N_stars, z=0.2)
+N_stars = len(stars)
+stellar.particles.add_particles(stars)
+stellar.commit_particles()
+
+Time_stars = np.linspace(t_min, t_max, N_t)
+
+def BVR_abs_stars(stellar, stars):
+   B_abs_stars = np.empty((N_t,N_stars))
+   V_abs_stars = np.empty((N_t,N_stars))
+   R_abs_stars = np.empty((N_t,N_stars))
+   for i in range(N_t):
+      t = Time_stars[i]
+      stellar.evolve_model(t|u.Myr)
+      stellar.update_particles(stars)
+      B_abs_stars[i], V_abs_stars[i], R_abs_stars[i] = extract_BVR_stars(stars)
+      
+   return B_abs_stars, V_abs_stars, R_abs_stars
+ 
+
+B_abs_stars, V_abs_stars, R_abs_stars = BVR_abs_stars(stellar, stars)
+
+def V(t, mu):
+   i = np.argmin(np.abs(Time_stars-t))
+   return V_abs_stars[i]+mu
+
+# * Theoretical values
+BV_MS = B_MS-V_MS
+
+BV_Cluster = B_Cluster-V_Cluster
+
+print("\033[94m"+"Info: done!"+"\033[0m")
+
+# * Figure and fitting
+
+def print_results(t_fit, dt_fit, d_fit, dd_fit, ext, ext_BV, ext_V):
+   print("\033[36m"+"t = "+"\033[0m" + "{:.2f}±{:.2f}".format(t_fit, dt_fit)+"\033[36m"+" Myr"+"\033[0m")
+   print("\033[36m"+"d = "+"\033[0m" + "{:.2f}±{:.2f}".format(d_fit/1000, dd_fit/1000)+"\033[36m"+" kpc"+"\033[0m")
+   print("\033[36m"+"ext = "+"\033[0m" + "{:.2f}".format(ext)+"\033[36m"+" mag"+"\033[0m")
+   print("\033[36m"+"E(B-V) = "+"\033[0m" + "{:.2f}".format(ext_BV)+"\033[36m"+" mag"+"\033[0m")
+   print("\033[36m"+"A_V = "+"\033[0m" + "{:.2f}".format(ext_V)+"\033[36m"+" mag"+"\033[0m")
+   return None
+
+def mu(d): return 5*np.log10(d)-5
+fig1, ax1 = plt.subplots(1)
+fig4, ax4 = plt.subplots(1)
+
+w = np.argwhere(ext_query.distances.to(u2.pc).value < 5000)
+ext_tot = ext_query(SkyCoord(ra*u2.deg, de*u2.deg, frame='icrs'), mode='mean')
+ax4.plot(ext_query.distances.to(u2.pc).value[w], ext_tot[w], color="C0")
+
+sct = ax4.plot(0,0, color="C4", marker="+", linestyle="", alpha=1, markersize=5)
+sctp = ax4.plot(0,0, color="C1", marker="+", linestyle="", alpha=0.65, markersize=3)
+sctm = ax4.plot(0,0, color="C1", marker="+", linestyle="", alpha=0.65, markersize=3)
+
+
+ax4.set_xlabel("Distance [pc]")
+ax4.set_ylabel("Reddening")
+
+
+t_fit = 0
+dt_fit = 0
+d_fit = 10
+dd_fit = 0
+ext = 0
+ext_BV = 0
+ext_V = 0
+
+continue_fit = True
+auto_ext = True
+
+ax1.scatter(BV_Cluster, V_Cluster, s=1, marker='.', color="C3")
+sct_BV = ax1.scatter(BV_MS, V_MS, s=1, marker='.', color="C0")
+
+
+X_BV = B_abs_stars[0] - V_abs_stars[0]
+
+Y = V(t_fit, mu(d_fit))
+Yp = V(t_fit+dt_fit, mu(d_fit-dd_fit))
+Ym = V(t_fit-dt_fit, mu(d_fit+dd_fit))
+
+fit_BV = ax1.plot(X_BV, Y, color="C4", marker=".", linestyle="", markersize=2)
+err_BVp = ax1.plot(X_BV, Yp, color="C1", marker=".", linestyle="", alpha=0.65, markersize=1.5)
+err_BVm = ax1.plot(X_BV, Ym, color="C1", marker=".", linestyle="", alpha=0.65, markersize=1.5)
+
+ax1.set_xlabel("Color $B-V$")
+ax1.set_ylabel("Magnitude $V$")
+ax1.invert_yaxis()
+ax1.legend([sct_BV, fit_BV[0]], ["Observational data (main sequence)", "$t =$ {:.2f} Myr, $d =$ {:.2f} kpc".format(t_fit, d_fit/1000)], loc="upper right")
+fig1.tight_layout()
+
+plt.show(block=False)
+
+print("\033[36m"+"Actions:"+"\033[0m")
+print("\033[36m"+"\t- t+/t- [int/float, Myr] (time increase/decrease)"+"\033[0m")
+print("\033[36m"+"\t- d+/d- [int/float, pc] (distance increase/decrease)"+"\033[0m")
+print("\033[36m"+"\t- vt+/vt- [int/float, Myr] (time uncertainty increase/decrease)"+"\033[0m")
+print("\033[36m"+"\t- vd+/vd- [int/float, pc] (distance uncertainty increase/decrease)"+"\033[0m")
+print("\033[36m"+"\t- ext+/ext- [int/float] (extinction value)"+"\033[0m")
+print("\033[36m"+"\t- ext auto [int/float] (auto extinction value)"+"\033[0m")
+print("\033[36m"+"\t- adjust (adjust plots)"+"\033[0m")
+print("\033[36m"+"\t- print (show results)"+"\033[0m")
+print("\033[36m"+"\t- ok (end)"+"\033[0m")
+print("\033[36m"+"\t- cancel (end)"+"\033[0m")
+
+
+while continue_fit:
+   var = input(("\033[32m"+"Enter action: "+"\033[0m"))
+   if var == "":
+      continue
+   elif var == "ok":
+      continue_fit = False
+   elif var == "cancel":
+      continue_fit = False
+      plt.close()
+      exit()
+   elif var == "print":
+      print_results(t_fit, dt_fit, d_fit, dd_fit, ext, ext_BV, ext_V)
+   elif var == "adjust":
+      ax1.set_ylim(np.max(V_Cluster), np.min(V_Cluster))
+      ax1.set_xlim(np.min(B_Cluster-V_Cluster), np.max(B_Cluster-V_Cluster))
+   elif var[0:3] == "ext":
+      if "auto" in var:
+         auto_ext = True
+      else:
+         auto_ext = False
+         if var[3] == "=":
+            val = float(var[4:])
+            ext = val
+         else:
+            val = float(var[3:])
+            ext += val
+         ext = np.clip(ext, 0, 10)
+         extp = ext
+         extm = ext
+   elif var[0:2] == "dd": 
+      print("\033[33m"+"Warning: please use vd instead!"+ "\033[0m")
+      if var[2] == "=":
+         val = float(var[3:])
+         dd_fit = val
+      else:
+         val = float(var[2:])
+         dd_fit += val
+      dd_fit = np.clip(dd_fit, 0, d_fit)
+   elif var[0:2] == "dt": 
+      print("\033[33m"+"Warning: please use vt instead!"+ "\033[0m")
+      if var[2] == "=":
+         val = float(var[3:])
+         dt_fit = val
+      else:
+         val = float(var[2:])
+         dt_fit += val
+      dt_fit = np.clip(dt_fit, 0, t_fit)
+   elif var[0] == "d": 
+      if var[1] == "=":
+         val = float(var[2:])
+         d_fit = val
+      else:
+         val = float(var[1:])
+         d_fit += val
+      d_fit = np.clip(d_fit, 1, 1e12)
+   elif var[0] == "t":
+      if var[1] == "=":
+         val = float(var[2:])
+         t_fit = val
+      else:
+         val = float(var[1:])
+         t_fit += val
+      t_fit = np.clip(t_fit, t_min, t_max)
+   elif var[0] == "v":
+      if var[1] == "d":
+         if var[2] == "=":
+            val = float(var[3:])
+            dd_fit = val
+         else:
+            val = float(var[2:])
+            dd_fit += val
+         dd_fit = np.clip(dd_fit, 0, d_fit)
+      if var[1] == "t":
+         if var[2] == "=":
+            val = float(var[3:])
+            dt_fit = val
+         else:
+            val = float(var[2:])
+            dt_fit += val
+         dt_fit = np.clip(dt_fit, 0, t_fit)
+   else:
+      print("\033[33m"+"Warning: not found!"+ "\033[0m")
+
+
+   # * Extinction computation
+   i = np.argmin(np.abs(Time_stars-t_fit))
+   ip = np.argmin(np.abs(Time_stars-(t_fit+dt_fit)))
+   im = np.argmin(np.abs(Time_stars-(t_fit-dt_fit)))
+   X_BV = B_abs_stars[i] - V_abs_stars[i]
+   Xp_BV = B_abs_stars[ip] - V_abs_stars[ip]
+   Xm_BV = B_abs_stars[im] - V_abs_stars[im]
+   
+   if auto_ext:
+      Coords_Cluster = SkyCoord(ra*u2.deg, de*u2.deg, distance=d_fit*u2.pc, frame='icrs')
+      ext = ext_query(Coords_Cluster, mode='mean') 
+      Coords_Clusterp = SkyCoord(ra*u2.deg, de*u2.deg, distance=(d_fit-dd_fit)*u2.pc, frame='icrs')
+      extp = ext_query(Coords_Clusterp, mode='mean') 
+      Coords_Clusterm = SkyCoord(ra*u2.deg, de*u2.deg, distance=(d_fit+dd_fit)*u2.pc, frame='icrs')
+      extm = ext_query(Coords_Clusterm, mode='mean') 
+
+   ext_BV = coeff_ext_BV * ext
+   ext_V = coeff_ext_V * ext
+   X_BV_corr = X_BV + ext_BV
+
+   extp_BV = coeff_ext_BV * extp
+   extp_V = coeff_ext_V * extp
+   Xp_BV_corr = Xp_BV + extp_BV
+   
+   extm_BV = coeff_ext_BV * extm
+   extm_V = coeff_ext_V * extm
+   Xm_BV_corr = Xm_BV + extm_BV
+
+   Y = V(t_fit, mu(d_fit)) + ext_V
+   Yp = V(t_fit+dt_fit, mu(d_fit-dd_fit)) + extp_V
+   Ym = V(t_fit-dt_fit, mu(d_fit+dd_fit)) + extm_V
+
+   fit_BV[0].set_xdata(X_BV_corr)
+
+   fit_BV[0].set_ydata(Y)
+
+   err_BVp[0].set_xdata(Xp_BV_corr)
+   err_BVm[0].set_xdata(Xm_BV_corr)
+
+   err_BVp[0].set_ydata(Yp)
+   err_BVm[0].set_ydata(Ym)
+
+   sct[0].set_xdata([d_fit])
+   sct[0].set_ydata([ext])
+   sctp[0].set_xdata([d_fit-dd_fit])
+   sctp[0].set_ydata([extp])
+   sctm[0].set_xdata([d_fit+dd_fit])
+   sctm[0].set_ydata([extm])
+
+
+   ax1.legend([sct_BV, fit_BV[0]], ["Observational data", "$t =$ {:.2f} $\\pm$ {:.2f} Myr, $d =$ {:.2f} $\\pm$ {:.2f} kpc, $\\mathrm{{ext}} =$ {:.2f}".format(t_fit, dt_fit, d_fit/1000, dd_fit/1000, ext)], loc="upper right")
+
+   fig1.canvas.draw_idle()
+   fig1.canvas.flush_events()
+
+   fig4.canvas.draw_idle()
+   fig4.canvas.flush_events()
+
+
+print_results(t_fit, dt_fit, d_fit, dd_fit, ext, ext_BV, ext_V)
+
+
+save_fig = input("\033[32m"+"Save figure (y/n): "+"\033[0m").upper() in YES
+
+if save_fig:
+   fig1.savefig(fig_directory+source+"_HR_B-V_fit")
+   fig4.savefig(fig_directory+source+"_HR_ext_fit")
+
+plt.show(block=False)
+input(("\033[33m"+"(Press enter to quit)"+"\033[0m"))
+plt.close()
+

JointSEx_v2.py

@@ -0,0 +1,152 @@
+"""
+* JointSEx : Joint SExtractor files
+* Version 2 - 2024-03-21
+@ Yaël Moussouni
+@ Unistra, P&E, MSc1-MoFP
+@ Observatory of Strasbourg (Intern)
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+import os
+from matplotlib.ticker import FuncFormatter, MaxNLocator
+
+# * Parameters
+
+if "YII_light" in plt.style.available: plt.style.use("YII_light")
+
+# * Variables
+in_directory = "./SEx_files/"
+out_directory = "./Joint_files/"
+fig_directory = "../Figs/JointSEx_figs/"
+SEx_prefix = "S-ex "
+B_suffix = "-B_cal.cat"
+V_suffix = "-G_cal.cat"
+R_suffix = "-R_cal.cat"
+BVR_suffix = "_BVR.csv"
+out_prefix = "Tab_"
+precision = 1/3600 # °
+
+YES = ["Y", "YES", "1", "OUI"]
+
+i_ra = 13
+i_de = 14
+i_mag = 1
+i_dmag = 2
+# MAG_ISO:  1, 2
+# MAG_BEST: 9, 10
+# TODO APERTURE
+
+# * Functions
+def calibration_function(I,c): return I + c 
+
+# deg = FuncFormatter(lambda x, pos=None: f"{x:.0f}\N{DEGREE SIGN}{(x-np.floor(x))*60:02.0f}\'00\"")
+# hrs = FuncFormatter(lambda x, pos=None: f"{x*24/360:.00f}:{(x*24/360-np.floor(x*24/360))*60:02.0f}:00")
+
+def DMS_from_deg_format(x, pos=None):
+   D = int(np.floor(x))
+   M = int(np.floor((x-D)*60))
+   S = int(np.floor(((x-D)*60-M)*60))
+   return "{D:02d}°{M:02d}\'{S:02d}\"".format(D=D,M=M,S=S)
+def HMS_from_deg_format(x, pos=None):
+   y = x/15
+   H = int(np.floor(y))
+   M = int(np.floor((y-H)*60))
+   S = int(np.floor(((y-H)*60-M)*60))
+   return "{H:02d}:{M:02d}:{S:02d}".format(H=H,M=M,S=S)
+
+DMS_from_deg = FuncFormatter(DMS_from_deg_format)
+HMS_from_deg = FuncFormatter(HMS_from_deg_format)
+
+def openfile(in_dir, SEx_prefix, source, suffix):
+   tab = []
+   with open(in_dir+SEx_prefix+source+suffix) as file:
+      for i in range(20):
+         file.readline()
+      for line in file:
+         columns = line.strip().split()
+         tab.append([float(columns[i_ra]), float(columns[i_de]), float(columns[i_mag]), float(columns[i_dmag])])
+      file.close()
+   return np.array(tab)
+
+# * Import 
+filelist = np.sort(os.listdir(in_directory))
+filelist = [f for f in filelist if SEx_prefix in f]
+sources = []
+for i in range(len(filelist)):
+   k = filelist[i][5:].find("-")
+   source = filelist[i][5:k+5]
+   if source not in sources: sources.append(source)
+
+print("\033[36m"+"Available source:"+"\033[0m")
+for i in range(len(sources)):
+   print("\033[36m"+"\t{}: ".format(i+1) + "\033[0m" + "{}".format(sources[i]))
+k = int(input("\033[32m"+"Select source number: "+"\033[0m"))-1
+print("")
+source = sources[k]
+
+tab_B = openfile(in_directory, SEx_prefix, source, B_suffix)
+tab_V = openfile(in_directory, SEx_prefix, source, V_suffix)
+tab_R = openfile(in_directory, SEx_prefix, source, R_suffix)
+
+w1 = np.where(tab_B[:,0] -  tab_V[:,0] < precision)
+w2 = np.where(tab_V[:,0] -  tab_R[:,0] < precision)
+w3 = np.where(tab_B[:,0] -  tab_R[:,0] < precision)
+
+
+if not (len(w1[0]) == len(w2[0]) == len(w3[0]) == len(tab_B) == len(tab_V) == len(tab_R)):
+   print("\033[33m"+"Warning: not all tables have the same size, or a star may have moved !"+ "\033[0m")
+   print("\033[33m"+"\tB table size:"+ "\033[0m"+"{:d}".format(len(tab_B)))
+   print("\033[33m"+"\tV table size:"+ "\033[0m"+"{:d}".format(len(tab_V)))
+   print("\033[33m"+"\tR table size:"+ "\033[0m"+"{:d}".format(len(tab_R)))
+   print("\033[33m"+"\tSame position in B-V:"+ "\033[0m"+"{:d}".format(len(w1[0])))
+   print("\033[33m"+"\tSame position in V-R:"+ "\033[0m"+"{:d}".format(len(w2[0])))
+   print("\033[33m"+"\tSame position in B-R:"+ "\033[0m"+"{:d}".format(len(w3[0])))
+else:
+   N = len(tab_R)
+   tab_BVR = [] # ra, dec, B, V, R, dB, dV, dR
+   for i in range(N):
+      ra = np.mean([tab_B[i,0], tab_V[i,0], tab_R[i,0]])
+      de = np.mean([tab_B[i,1], tab_V[i,1], tab_R[i,1]])
+      B = tab_B[i,2]
+      V = tab_V[i,2]
+      R = tab_R[i,2]
+      dB = tab_B[i,3]
+      dV = tab_V[i,3]
+      dR = tab_R[i,3]
+      tab_BVR.append([ra, de, B, V, R, dB, dV, dR])
+   tab_BVR = np.array(tab_BVR)
+
+# * Output choices
+
+plot_map = input("\033[32m"+"Draw map (y/n): "+"\033[0m").upper() in YES
+plot_save = input("\033[32m"+"Save figures (y/n): "+"\033[0m").upper() in YES
+print("")
+save_output = input("\033[32m"+"Save output file (y/n): "+"\033[0m").upper() in YES
+print("")
+
+# * Map
+if plot_map:
+   fig, ax = plt.subplots(1)
+   ax.scatter(tab_BVR[:,0], tab_BVR[:,1], s=2, color='C3', alpha=1)
+   ax.xaxis.set_major_formatter(HMS_from_deg)
+   ax.xaxis.set_major_locator(MaxNLocator(5))
+   ax.yaxis.set_major_formatter(DMS_from_deg)
+   ax.set_aspect("equal")
+   ax.set_xlabel("Right ascension")
+   ax.set_ylabel("Declination")
+   if plot_save: 
+      fig.savefig(fig_directory+source+"_map")
+      print("\033[94m"+"Info: Figure saved as " + "\033[0m" + source+"_map")
+
+# * File output
+if save_output:
+   file = open(out_directory+out_prefix+source+BVR_suffix, "w")
+   file.write("ra [deg],dec [deg],B_inst [mag],V_inst [mag],R_inst [mag],dB_inst [mag],dV_inst [mag],dR_inst [mag]\n")
+   for i in range(len(tab_BVR)):
+      ra, dec, B, V, R, dB, dV, dR = tab_BVR[i,:]
+      file.write("{ra:e},{dec:e},{B:e},{V:e},{R:e},{dB:e},{dV:e},{dR:e}\n".format(ra=ra, dec=dec, B=B, V=V, R=R, dB=dB, dV=dV, dR=dR))
+   print("\033[94m"+"Info: File saved as " + "\033[0m" + out_prefix+source+BVR_suffix)
+   file.close()
+plt.show()
+
+

PhotoCal_v2.py

@@ -0,0 +1,180 @@
+"""
+* PhotoCal: Extract the instrumental magnitude coefficients and perform photometry calibration
+* Version 2 - 2024-03-21
+@ Yaël Moussouni
+@ Unistra, P&E, MSc1-MoFP
+@ Observatory of Strasbourg (Intern)
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+import os
+from scipy.optimize import curve_fit
+
+if "YII_light" in plt.style.available: plt.style.use("YII_light")
+
+Simbad_directory = "./Simbad_files/"
+Joint_directory = "./Joint_files/"
+PhotoCal_directory = "./PhotoCal_files/"
+fig_directory = "../Figs/PhotoCal_figs/"
+Joint_prefix = "Joint_"
+Simbad_prefix = "Simbad_"
+PhotoCal_prefix = "PhotoCal_"
+BVR_suffix = "_BVR"
+filtered_suffix = "_filtered"
+file_ext = ".csv"
+
+precision = 1/3600 # °
+
+YES = ["Y", "YES", "1", "OUI"]
+
+# * Functions
+
+def calibration_function(I,c): return I + c 
+
+# * Import
+
+Joint_filelist = np.sort(os.listdir(Joint_directory))
+Joint_filelist = [f for f in Joint_filelist if Joint_prefix in f]
+filtered_filelist = np.sort(os.listdir(Joint_directory))
+filtered_filelist = [f for f in filtered_filelist if (Joint_prefix in f and filtered_suffix in f)]
+Simbad_filelist = np.sort(os.listdir(Simbad_directory))
+Simbad_filelist = [f for f in Simbad_filelist if Simbad_prefix in f]
+Joint_sources = []
+filtered_sources = []
+Simbad_sources = []
+
+for i in range(len(Joint_filelist)):
+   k = Joint_filelist[i][len(Joint_prefix):].find("_BVR")
+   source = Joint_filelist[i][len(Joint_prefix):k+len(Joint_prefix)]
+   if source not in Joint_sources: Joint_sources.append(source)
+for i in range(len(filtered_filelist)):
+   k = filtered_filelist[i][len(Joint_prefix):].find("_BVR")
+   source = filtered_filelist[i][len(Joint_prefix):k+len(Joint_prefix)]
+   if source not in filtered_sources: filtered_sources.append(source)
+for i in range(len(Simbad_filelist)):
+   k = Simbad_filelist[i][len(Simbad_prefix):].find("_BVR")
+   source = Simbad_filelist[i][len(Simbad_prefix):k+len(Simbad_prefix)]
+   if source not in Simbad_sources: Simbad_sources.append(source)
+sources = np.intersect1d(Joint_sources, Simbad_sources)
+
+print("\033[36m"+"Available source:"+"\033[0m")
+for i in range(len(sources)):
+   if sources[i] in filtered_sources:  print("\033[36m"+"\t{}: ".format(i+1) + "\033[0m" + "{} (filtered)".format(sources[i]))
+   else: print("\033[36m"+"\t{}: ".format(i+1) + "\033[0m" + "{}".format(sources[i]))
+k = int(input("\033[32m"+"Select source number: "+"\033[0m"))-1
+print("")
+source = sources[k]
+
+i_ra = 0
+i_de = 1
+i_B_inst = 2
+i_V_inst = 3
+i_R_inst = 4
+i_dB_inst = 5
+i_dV_inst = 6
+i_dR_inst = 7
+i_B_ref = 37
+i_V_ref = 38
+i_R_ref = 39
+
+cols = [i_ra, i_de, i_B_inst, i_V_inst, i_R_inst, i_dB_inst, i_dV_inst, i_dR_inst, i_B_ref, i_V_ref, i_R_ref]
+
+
+Simbad_file = np.genfromtxt(Simbad_directory+Simbad_prefix+source+BVR_suffix+file_ext, skip_header=1, comments="#", delimiter=",", usecols=cols)
+if source in filtered_sources: Joint_file = np.loadtxt(Joint_directory+Joint_prefix+source+BVR_suffix+filtered_suffix+file_ext, skiprows=1, comments="#", delimiter=",")
+else: Joint_file = np.loadtxt(Joint_directory+Joint_prefix+source+BVR_suffix+file_ext, skiprows=1, comments="#", delimiter=",")
+
+# * Calibration
+w_B = np.intersect1d(np.argwhere(np.isfinite(Simbad_file[:,8])), np.argwhere(Simbad_file[:,2]!=99))
+w_V = np.intersect1d(np.argwhere(np.isfinite(Simbad_file[:,9])), np.argwhere(Simbad_file[:,3]!=99))
+w_R = np.intersect1d(np.argwhere(np.isfinite(Simbad_file[:,10])), np.argwhere(Simbad_file[:,4]!=99))
+
+c_B, dc_B = curve_fit(calibration_function, Simbad_file[w_B,2].flatten(), Simbad_file[w_B,8].flatten())
+c_V, dc_V = curve_fit(calibration_function, Simbad_file[w_V,3].flatten(), Simbad_file[w_V,9].flatten())
+c_R, dc_R = curve_fit(calibration_function, Simbad_file[w_R,4].flatten(), Simbad_file[w_R,10].flatten())
+
+dc_B = np.sqrt(np.diag(dc_B))
+dc_V = np.sqrt(np.diag(dc_V))
+dc_R = np.sqrt(np.diag(dc_R))
+
+Simbad_file[:,2] += c_B[0]
+Simbad_file[:,3] += c_V[0]
+Simbad_file[:,4] += c_R[0]
+
+Simbad_file[:,5] += dc_B[0]
+Simbad_file[:,6] += dc_V[0]
+Simbad_file[:,7] += dc_R[0]
+
+Joint_file[:,2] += c_B[0]
+Joint_file[:,3] += c_V[0]
+Joint_file[:,4] += c_R[0]
+
+Joint_file[:,5] += dc_B[0]
+Joint_file[:,6] += dc_V[0]
+Joint_file[:,7] += dc_R[0]
+
+# * Output choices
+
+plot_calibration = input("\033[32m"+"Calibration plot (y/n): "+"\033[0m").upper() in YES
+plot_magnitude =  input("\033[32m"+"Show magnitude distribution (y/n): "+"\033[0m").upper() in YES
+plot_save = input("\033[32m"+"Save figures (y/n): "+"\033[0m").upper() in YES
+print("")
+save_output = input("\033[32m"+"Save output file (y/n): "+"\033[0m").upper() in YES
+print("")
+
+# * Calibration plot
+tab_mag = np.concatenate([Simbad_file[w_B,2].flatten(),Simbad_file[w_B,8].flatten(),Simbad_file[w_V,3].flatten(),Simbad_file[w_V,9].flatten(),Simbad_file[w_R,3].flatten(),Simbad_file[w_R,10].flatten()])
+
+if plot_calibration:
+   fig, ax = plt.subplots(1)
+   ax.errorbar(Simbad_file[w_B,8].flatten(), Simbad_file[w_B,2].flatten(), Simbad_file[w_B,5].flatten(), markersize=2, color="#0000FF", linestyle="", marker="o", label='B filter', alpha=0.5)
+   ax.errorbar(Simbad_file[w_V,9].flatten(), Simbad_file[w_V,3].flatten(), Simbad_file[w_V,6].flatten(), markersize=2, color="#00FF00", linestyle="", marker="s", label="V filter", alpha=0.5)
+   ax.errorbar(Simbad_file[w_R,10].flatten(), Simbad_file[w_R,4].flatten(), Simbad_file[w_R,7].flatten(), markersize=2, color="#FF0000", linestyle="", marker="^", label="R filter", alpha=0.5)
+   
+   min_CAL = np.min(tab_mag)
+   max_CAL = np.max(tab_mag)
+   ax.plot([min_CAL-0.5,max_CAL+0.5], [min_CAL-0.5,max_CAL+0.5], "--", alpha=0.5)
+   ax.set_xlabel("Reference magnitude")
+   ax.set_ylabel("Calibrated magnitude")
+   ax.legend()
+   fig.tight_layout()
+   if plot_save: 
+      fig.savefig(fig_directory+source+"_calibration_plot")
+      print("\033[94m"+"Info: Figure saved as " + "\033[0m" + source+"_calibration_plot")
+
+# * Magnitude histograms
+
+if plot_magnitude:
+   fig, ax = plt.subplots(1)
+   bins = np.linspace(np.quantile(tab_mag, 0.1), np.quantile(tab_mag,0.9), 20)
+   hist_B, bins_B = np.histogram(Joint_file[:,2],bins)
+   hist_V, bins_V = np.histogram(Joint_file[:,3],bins)
+   hist_R, bins_R = np.histogram(Joint_file[:,4],bins)
+
+   cbins_B = 0.5*(bins_B[1:] + bins_B[:-1])
+   cbins_V = 0.5*(bins_V[1:] + bins_V[:-1])
+   cbins_R = 0.5*(bins_R[1:] + bins_R[:-1])
+
+   ax.scatter(cbins_B, hist_B, color="#0000FF", label="B filter", marker="o", alpha=0.5)
+   ax.scatter(cbins_V, hist_V, color="#00FF00", label="V filter", marker="s", alpha=0.5)
+   ax.scatter(cbins_R, hist_R, color="#FF0000", label="R filter", marker="^", alpha=0.5)
+
+   ax.set_xlabel("Relative magnitude")
+   ax.set_ylabel("Star distribution")
+   ax.legend()
+   fig.tight_layout()
+   if plot_save: 
+      fig.savefig(fig_directory+source+"_magnitude_hist")
+      print("\033[94m"+"Info: Figure saved as " + "\033[0m" + source+"_magnitude_hist")
+
+if save_output:
+   file = open(PhotoCal_directory+PhotoCal_prefix+source+BVR_suffix+file_ext, "w")
+   file.write("ra [deg],dec [deg],B_obs [mag],V_obs [mag],R_obs [mag],dB_obs [mag],dV_obs [mag],dR_obs [mag]\n")
+   for i in range(len(Joint_file)):
+      ra, dec, B, V, R, dB, dV, dR = Joint_file[i,:]
+      file.write("{ra:e},{dec:e},{B:e},{V:e},{R:e},{dB:e},{dV:e},{dR:e}\n".format(ra=ra, dec=dec, B=B, V=V, R=R, dB=dB, dV=dV, dR=dR))
+   print("\033[94m"+"Info: File saved as " + "\033[0m" + PhotoCal_prefix+source+BVR_suffix+file_ext)
+   file.close()
+plt.show()
+

SCOPE_v1.py

@@ -0,0 +1,73 @@
+"""
+* SCOPE : Sky Coordinates for Observations Python Estimator
+* Version 1 - 2024-02-09
+@ Yaël Moussouni
+@ Unistra, P&E, MSc1-MoFP
+@ Observatory of Strasbourg (Intern)
+"""
+import astropy.time as time
+import astropy.coordinates as coord
+import astropy.units as u
+
+# * Params
+eps = coord.Angle("0d0m1s")
+# * Location
+obs_name = "Strasbourg Observatory"
+obs_lat = coord.Latitude("""48d34m59s""")  # @ 48°34'59" N
+obs_lon = coord.Longitude("""07d46m07s""") # @ 07°46'07" E
+
+# * User input
+print("")
+print("\033[94m"+"Info: Location is manually set."+"\033[0m")
+print("")
+min_alt = coord.Angle(input("\033[32m"+"Minimum altitude above horizon (deg): "+"\033[0m")+"d")
+min_dec = coord.Angle(input("\033[32m"+"Minimum declination around the north (deg): "+"\033[0m")+"d")
+window_east = coord.Angle(input("\033[32m"+"Observation window before crossing meridian (h): "+"\033[0m")+" hours")
+window_west = coord.Angle(input("\033[32m"+"Observation window after crossing meridian (h): "+"\033[0m")+" hours")
+obs_date = input("\033[32m"+"Date of observation (YYYY-MM-DD format): "+"\033[0m")
+obs_utc = input("\033[32m"+"Time of observation (hh:mm or hh:mm:ss formats - UTC): "+"\033[0m")
+obs_duration = coord.Angle(input("\033[32m"+"Observation duration (h): "+"\033[0m")+" hours")
+obs_sky_rotation = obs_duration * (1*u.sday/u.day).decompose()
+
+# * Manual input
+# min_alt = coord.Angle("45"+"d")
+# min_dec = coord.Angle("30"+"d")
+# window_east = coord.Angle("3"+" hours")
+# window_west = coord.Angle("1"+" hours")
+# obs_date ="2024-02-08"
+# obs_utc = "20:00"
+# obs_duration = coord.Angle("1"+" hours")
+   
+# * Date/time/location objects
+obs_loc = coord.EarthLocation(lon=obs_lon, lat=obs_lat)
+obs_utc = time.Time(obs_date + " " + obs_utc, scale='utc')
+obs = time.Time(obs_utc, location=obs_loc)
+
+# * Local sidereal time
+st = obs.sidereal_time("apparent") # ? apparent or absolute
+
+# * Borders and corners coordinates
+dec_upper = min_dec
+dec_lower = + min_alt + obs_loc.lat - coord.Angle("90d") 
+a_west = st + obs_sky_rotation + window_west
+a_east = st - window_east
+
+# * Output
+print("")
+print("\033[36m"+"Location: "+"\033[0m"+ obs_name)
+print("\033[36m"+"\tlat: "+"\033[0m" + obs_lat.to_string())
+print("\033[36m"+"\tlon: "+"\033[0m" + obs_lon.to_string())
+print("\033[36m"+"Observation "+"\033[0m")
+print("\033[36m"+"\tdate and time: "+"\033[0m"+ obs_utc.to_string())
+print("\033[36m"+"\tsidereal time: "+"\033[0m"+ st.to_string(unit=u.hour))
+print("\033[36m"+"\tduration:      "+"\033[0m"+ obs_duration.to_string(unit=u.hour))
+print("\033[36m"+"\tsky rotation:  "+"\033[0m"+ obs_sky_rotation.to_string(unit=u.hour))
+print("\033[36m"+"Sky coordinates window:")
+print("\033[36m"+"\teast ra:   "+"\033[0m" + a_east.to_string(unit=u.hour))
+print("\033[36m"+"\twest ra:   "+"\033[0m" + a_west.to_string(unit=u.hour))
+print("\033[36m"+"\tupper dec: "+"\033[0m" + dec_upper.to_string(unit=u.degree))
+print("\033[36m"+"\tlower dec: "+"\033[0m" + dec_lower.to_string(unit=u.degree))
+
+print("\033[36m"+"Simbad constraints:"+"\033[0m")
+print("\033[36m"+"\tRA: "+"\033[0m" + "<" + a_west.to_string(unit=u.hour, sep=":") + " && " + ">" + a_east.to_string(unit=u.hour, sep=":"))
+print("\033[36m"+"\tDE: "+"\033[0m" + "<" + dec_upper.to_string(unit=u.degree, sep=":") + " && " + ">" + dec_lower.to_string(unit=u.degree, sep=":"))

YBC_tables/ubvrijhklm/extinctionodonnell94Rv3.1.al

@@ -0,0 +1 @@
+1.55814  1.32616  1.00096  0.80815  0.59893  0.28688  0.18103  0.11265  0.05415  0.03306  #extinctions with Av=1.

YBC_tables/ubvrijhklm/filter.info

@@ -0,0 +1,11 @@
+#  I   Filter  File                Lamb_C       D_Lamb     Flam_0     Sys   Dev   Mag_zp      NP comments
+   1 U          u_maiz.dat          3598.5445     584.0610   4.2687e-09 Vega  CCD    0.0300      0 #
+   2 B          b_maiz.dat          4385.9244     912.7500   6.3606e-09 Vega  CCD    0.0300      0 #
+   3 V          v_maiz.dat          5490.5551     857.4373   3.6168e-09 Vega  CCD    0.0300      0 #
+   4 R          r_bessell.dat       6594.7209    1591.1579   2.0973e-09 Vega  CCD    0.0300      0 #
+   5 I          i_bessell.dat       8059.8822    1495.0994   1.1123e-09 Vega  CCD    0.0300      0 #
+   6 J          j_bessell.dat      12369.2625    2046.5809   3.0916e-10 Vega  CCD    0.0300      0 #
+   7 H          h_bessell.dat      16464.4483    2854.1417   1.1256e-10 Vega  CCD    0.0300      0 #
+   8 K          k_bessell.dat      22105.4456    3646.1890   3.8958e-11 Vega  CCD    0.0300      0 #
+   9 L          l_bessell.dat      34840.7510    4595.4536   6.9824e-12 Vega  CCD    0.0300      0 #
+  10 M          m_bessell.dat      47339.3251    3559.9992   2.1226e-12 Vega  CCD    0.0300      0 #