added limit t12 (4c1cdc4c) · Commits · GUIBRUSHR / DataReductionGIANOB

DataReductionGIANOB/DataInterface.py

+202 −66

Original line number	Diff line number	Diff line
		import os

		import numpy as np
		import yaml
		from numpy import array
		from os import environ
		@@ -66,9 +67,10 @@ def get_target_info(path_targets_folder, target):
		A tuple containing target parameters in the following order:
		(name, radius_jup, mass_jup, gravity_cm_s2, t_eq_K, stellar_radius_sun,
		stellar_mass_sun, p0_log10_bar, hjd0_days, period_days, v_system_km_s,
		Limit_Phase, kp_km_s, ecc, periastron_argument, stellar_teff_K, ra, dec,
		a_Rs_ratio, projected_obliquity_deg, inclination_deg, v_sini_km_s).
		Returns a list of 22 None values if target is "None".
		Limit_Phase_T14, kp_km_s, ecc, periastron_argument, stellar_teff_K, ra, dec,
		a_Rs_ratio, projected_obliquity_deg, inclination_deg, v_sini_km_s,
		Limit_Phase_T12, t14_hours, t12_hours).
		Returns a list of 25 None values if target is "None".
		"""
		if target != "None":
		with open(str(Path(path_targets_folder, target, "info.yaml")), 'r') as f:
		@@ -85,7 +87,7 @@ def get_target_info(path_targets_folder, target):
		yaml_file["hjd0_days"],
		yaml_file["period_days"],
		yaml_file["v_system_km_s"],
		yaml_file["Limit_Phase"],
		yaml_file.get("Limit_Phase_T14", yaml_file.get("Limit_Phase", 0)),
		yaml_file["kp_km_s"],
		yaml_file["ecc"],
		yaml_file["periastron_argument"],
		@@ -96,9 +98,12 @@ def get_target_info(path_targets_folder, target):
		yaml_file.get("projected_obliquity_deg", ""),
		yaml_file.get("inclination_deg", ""),
		yaml_file.get("v_sini_km_s", ""),
		yaml_file.get("Limit_Phase_T12", 0),
		yaml_file.get("t14_hours", 0),
		yaml_file.get("t12_hours", 0),
		)
		else:
		return [None for _ in range(22)]
		return [None for _ in range(25)]


		def get_target_nights(path_targets, target, rad_mode, instrument="GIANO_B", simulated=0):
		@@ -134,6 +139,110 @@ def get_target_nights(path_targets, target, rad_mode, instrument="GIANO_B", simu
		return nights


		def compute_transit_durations(period, a_Rs, k, inc_deg, ecc=0.0, omega=np.pi / 2, t14_def=None):
		"""
		Compute transit durations T14, T23, T12 (=T34) for circular and eccentric orbits.

		Parameters
		----------
		period : float
		Orbital period in days.
		a_Rs : float
		Semi-major axis in units of stellar radii (a / R_star).
		k : float
		Planet-to-star radius ratio (Rp / Rs).
		inc_deg : float
		Orbital inclination in degrees.
		ecc : float, optional
		Orbital eccentricity (default 0 for circular).
		omega : float, optional
		Argument of periastron of the STAR in radians (default pi/2).
		t14_def : float, optional
		Default transit duration T14 in hours (default None).

		Returns
		-------
		dict with keys:
		'T14', 'T23', 'T12' : durations in days
		'T14_hours', 'T23_hours', 'T12_hours' : durations in hours
		'b' : impact parameter
		'ecc_factor' : eccentricity correction factor (1.0 if circular)
		'grazing' : True if grazing transit (T23 = 0)
		'limph_T14', 'limph_T12' : phase limits = T / (2*P)

		References
		----------
		Winn (2010), "Transits and Occultations", arXiv:1001.2010, Eq. 14-16.
		Kipping (2010), MNRAS 407, 301, arXiv:1004.3819.
		"""
		inc = np.radians(inc_deg)
		sin_i = np.sin(inc)

		# Impact parameter
		# Eccentric: b = (a/Rs) * cos(i) * (1 - e^2) / (1 + e*sin(omega))
		# Circular: b = (a/Rs) * cos(i)
		if ecc > 0:
		ecc_factor_b = (1.0 - ecc ** 2) / (1.0 + ecc * np.sin(omega))
		else:
		ecc_factor_b = 1.0

		b = a_Rs * np.cos(inc) * ecc_factor_b

		if b > (1.0 + k):
		raise ValueError(
		f"No transit: b = {b:.4f} > 1 + k = {1 + k:.4f}."
		)

		# Winn (2010) Eq. 14-15
		arg_T14_sq = (1.0 + k) 2 - b 2
		arg_T23_sq = (1.0 - k) 2 - b 2

		if arg_T14_sq < 0:
		raise ValueError(
		f"No transit: (1+k)^2 - b^2 = {arg_T14_sq:.6f} < 0"
		)

		grazing = False
		if arg_T23_sq < 0:
		grazing = True
		arg_T23_sq = 0.0

		sin_T14 = min((1.0 / a_Rs) * np.sqrt(arg_T14_sq) / sin_i, 1.0)
		sin_T23 = min((1.0 / a_Rs) * np.sqrt(arg_T23_sq) / sin_i, 1.0)

		T14_circ = (period / np.pi) * np.arcsin(sin_T14)
		T23_circ = (period / np.pi) * np.arcsin(sin_T23)

		# Eccentricity correction (Winn 2010, Eq. 16)
		# Psi = sqrt(1 - e^2) / (1 + e*sin(omega))
		if ecc > 0:
		ecc_factor = np.sqrt(1.0 - ecc ** 2) / (1.0 + ecc * np.sin(omega))
		else:
		ecc_factor = 1.0

		T14 = T14_circ * ecc_factor
		T23 = T23_circ * ecc_factor
		T12 = (T14 - T23) / 2.0

		if t14_def is not None:
		print(f"T14 set to {t14_def} hours (provided) instead of {T14 * 24} hours (calculated)")
		T14 = t14_def / 24.0

		return {
		'T14': T14,
		'T23': T23,
		'T12': T12,
		'T14_hours': T14 * 24.0,
		'T23_hours': T23 * 24.0,
		'T12_hours': T12 * 24.0,
		'b': b,
		'ecc_factor': ecc_factor,
		'grazing': grazing,
		'limph_T14': T14 / (2.0 * period),
		'limph_T12': T12 / (2.0 * period),
		}


		def insert_target_catalog(
		path_target,
		target_id,
		@@ -147,10 +256,13 @@ def insert_target_catalog(
		hjd0,
		period,
		vsys,
		limph,
		limph_T14,
		limph_T12,
		kp,
		ks,
		ecc,
		t14,
		t12,
		opi=1.57,
		stellar_teff=0,
		rad_mode="Transmission",
		@@ -193,14 +305,20 @@ def insert_target_catalog(
		Orbital period in days.
		vsys : float
		System velocity in km/s.
		limph : float
		Phase limit for observations.
		limph_T14 : float
		Phase limit for T14 transit.
		limph_T12 : float
		Phase limit for T12 transit.
		kp : float
		Planet's semi-amplitude velocity in km/s.
		ks : float
		Star's semi-amplitude velocity in km/s.
		ecc : float
		Orbital eccentricity.
		t14 : float
		Total transit duration T14 in hours.
		t12 : float
		Ingress/egress duration T12 in hours.
		opi : float, optional
		Periastron argument in radians. Default is 1.57 (π/2).
		stellar_teff : float, optional
		@@ -225,8 +343,7 @@ def insert_target_catalog(
		None
		Creates directories and writes info.yaml file.
		"""
		yaml_dict = (
		{
		yaml_dict = {
		'name': target_id,
		'radius_jup': float(radius),
		'mass_jup': float(mass),
		@@ -238,10 +355,13 @@ def insert_target_catalog(
		'p0_log10_bar': float(p0),
		'hjd0_days': float(hjd0),
		'period_days': float(period),
		't14_hours': float(t14),
		't12_hours': float(t12),
		'Limit_Phase_T14': float(limph_T14),
		'Limit_Phase_T12': float(limph_T12),
		'kp_km_s': float(kp),
		'ks_km_s': float(ks),
		'v_system_km_s': float(vsys),
		'Limit_Phase': float(limph),
		'ecc': float(ecc),
		'periastron_argument': float(opi),
		'ra': float(ra),
		@@ -249,10 +369,8 @@ def insert_target_catalog(
		'a_Rs_ratio': float(a_Rs_ratio),
		'projected_obliquity_deg': float(projected_obliquity),
		'inclination_deg': float(inclination),
		'v_sini_km_s': float(v_sini)

		'v_sini_km_s': float(v_sini),
		}
		)
		os.system("mkdir -p " + str(Path(path_target, target_id, "Models")))
		os.system("mkdir -p " + str(Path(path_target, target_id, "Star_Spectrum")))
		os.system("mkdir -p " + str(Path(path_target, target_id, "Retrievals")))
		@@ -274,7 +392,7 @@ def update_target_catalog(
		hjd0,
		period,
		vsys,
		limph,
		limph_T14,
		kp,
		ecc,
		opi=1.57,
		@@ -284,7 +402,10 @@ def update_target_catalog(
		a_Rs_ratio=0,
		projected_obliquity=0,
		inclination=0,
		v_sini=0
		v_sini=0,
		limph_T12=0,
		t14=0,
		t12=0
		):
		"""
		Update an existing target's info.yaml file with new parameters.
		@@ -319,8 +440,8 @@ def update_target_catalog(
		Orbital period in days.
		vsys : float
		System velocity in km/s.
		limph : float
		Phase limit for observations.
		limph_T14 : float
		Phase limit for T14 transit.
		kp : float
		Planet's semi-amplitude velocity in km/s.
		ecc : float
		@@ -333,14 +454,27 @@ def update_target_catalog(
		Right ascension in degrees. Default is 0.
		dec : float, optional
		Declination in degrees. Default is 0.
		a_Rs_ratio : float, optional
		Semi-major axis to stellar radius ratio. Default is 0.
		projected_obliquity : float, optional
		Projected obliquity in degrees. Default is 0.
		inclination : float, optional
		Inclination in degrees. Default is 0.
		v_sini : float, optional
		Stellar velocity in km/s. Default is 0.
		limph_T12 : float, optional
		Phase limit for T12 transit. Default is 0.
		t14 : float, optional
		Total transit duration T14 in hours. Default is 0.
		t12 : float, optional
		Ingress/egress duration T12 in hours. Default is 0.

		Returns
		-------
		None
		Updates the info.yaml file for the target.
		"""
		yaml_dict = (
		{
		yaml_dict = {
		'name': target_id,
		'radius_jup': float(radius),
		'mass_jup': float(mass),
		@@ -352,9 +486,12 @@ def update_target_catalog(
		'p0_log10_bar': float(p0),
		'hjd0_days': float(hjd0),
		'period_days': float(period),
		't14_hours': float(t14),
		't12_hours': float(t12),
		'Limit_Phase_T14': float(limph_T14),
		'Limit_Phase_T12': float(limph_T12),
		'kp_km_s': float(kp),
		'v_system_km_s': float(vsys),
		'Limit_Phase': float(limph),
		'ecc': float(ecc),
		'periastron_argument': float(opi),
		'ra': float(ra),
		@@ -362,9 +499,8 @@ def update_target_catalog(
		'a_Rs_ratio': float(a_Rs_ratio),
		'projected_obliquity_deg': float(projected_obliquity),
		'inclination_deg': float(inclination),
		'v_sini_km_s': float(v_sini)
		'v_sini_km_s': float(v_sini),
		}
		)
		with open(str(Path(path_target, target_id, "info.yaml")), 'w') as f:
		yaml.dump(yaml_dict, f, sort_keys=False)

DataReductionGIANOB/Frame_Gofio.py

+18 −8

Original line number	Diff line number	Diff line
		@@ -82,7 +82,7 @@ class Frame_Gofio:
		)
		except Exception as e:
		nights = ""
		self.target_name, self.target_radius, self.target_mass, self.target_gravity, self.target_t_eq, self.target_stellar_radius, self.target_stellar_mass, self.target_p0, self.target_hjd0, self.target_period, self.target_vsys, self.target_limph, self.target_kp, self.target_ecc, self.target_opi, self.target_stellar_teff, self.target_ra, self.target_dec, self.target_a_Rs_ratio, self.target_projected_obliquity, self.target_inclination, self.target_v_sini = DataInterface.get_target_info(
		self.target_name, self.target_radius, self.target_mass, self.target_gravity, self.target_t_eq, self.target_stellar_radius, self.target_stellar_mass, self.target_p0, self.target_hjd0, self.target_period, self.target_vsys, self.target_limph, self.target_kp, self.target_ecc, self.target_opi, self.target_stellar_teff, self.target_ra, self.target_dec, self.target_a_Rs_ratio, self.target_projected_obliquity, self.target_inclination, self.target_v_sini, self.target_limph_T12, self.target_t14_hours, self.target_t12_hours = DataInterface.get_target_info(
		self.path_target, target_list[0])
		else:
		nights = ""
		@@ -1170,7 +1170,7 @@ class Frame_Gofio:
		)
		except:
		nights = ""
		self.target_name, self.target_radius, self.target_mass, self.target_gravity, self.target_t_eq, self.target_stellar_radius, self.target_stellar_mass, self.target_p0, self.target_hjd0, self.target_period, self.target_vsys, self.target_limph, self.target_kp, self.target_ecc, self.target_opi, self.target_stellar_teff, self.target_ra, self.target_dec, self.target_a_Rs_ratio, self.target_projected_obliquity, self.target_inclination, self.target_v_sini = DataInterface.get_target_info(
		self.target_name, self.target_radius, self.target_mass, self.target_gravity, self.target_t_eq, self.target_stellar_radius, self.target_stellar_mass, self.target_p0, self.target_hjd0, self.target_period, self.target_vsys, self.target_limph, self.target_kp, self.target_ecc, self.target_opi, self.target_stellar_teff, self.target_ra, self.target_dec, self.target_a_Rs_ratio, self.target_projected_obliquity, self.target_inclination, self.target_v_sini, self.target_limph_T12, self.target_t14_hours, self.target_t12_hours = DataInterface.get_target_info(
		self.path_target, target_id)
		self.night.delete("1.0", tk.END)
		self.night.insert(tk.END, nights)
		@@ -1763,7 +1763,7 @@ class Frame_Gofio:

		def wcal_auto(self):
		target = self.chosen_target_var.get()
		_, _, _, _, _, _, _, _, _, _, vsys, _, _, _, _, stellar_teff, _, _, _, _, _, _ = DataInterface.get_target_info(
		_, _, _, _, _, _, _, _, _, _, vsys, _, _, _, _, stellar_teff, _, _, _, _, _, _, _, _, _ = DataInterface.get_target_info(
		self.path_target, target)
		instrument = self.chosen_instrument_var_post.get()
		with open(str(Path(self.path_default, "DataReductionGIANOB", "Instruments", instrument + ".yaml")), "r") as f:
		@@ -1849,7 +1849,7 @@ class Frame_Gofio:

		def test_wlen_calib(self):
		target = self.chosen_target_var.get()
		_, _, _, _, _, _, _, _, _, _, vsys, _, _, _, _, stellar_teff, _, _, _, _, _, _ = DataInterface.get_target_info(
		_, _, _, _, _, _, _, _, _, _, vsys, _, _, _, _, stellar_teff, _, _, _, _, _, _, _, _, _ = DataInterface.get_target_info(
		self.path_target, target)
		instrument = self.chosen_instrument_var_post.get()
		path_instrument = str(
		@@ -2007,13 +2007,23 @@ class Frame_Gofio:
		kp = ((mass_st_sun * Msun * ks) / (mass_pl_jup_min * MJ))
		gravity = phys_const.G * (mass_pl_jup * phys_const.m_jup) / np.power(radius_pl_jup * phys_const.r_jup_mean,
		2)
		limph = float(self.T14.get()) / (24 * 2 * period)
		t14_hours_input = float(self.T14.get())
		phase_dict = DataInterface.compute_transit_durations(
		period, float(a_Rs_ratio),
		(radius_pl_jup * phys_const.r_jup_mean) / (float(radius_st_sun) * phys_const.r_sun),
		float(inclination), ecc=float(ecc), omega=float(opi), t14_def=t14_hours_input
		)
		limph_T14 = phase_dict["limph_T14"]
		limph_T12 = phase_dict["limph_T12"]
		t14_hours = phase_dict["T14_hours"]
		t12_hours = phase_dict["T12_hours"]

		DataInterface.insert_target_catalog(
		self.path_target, target, str(radius_pl_jup), str(mass_pl_jup), str(gravity), str(T_eq),
		str(radius_st_sun), str(mass_st_sun),
		str(P0), str(hjd0), str(period), str(vsys), str(limph), str(kp), str(ks), str(ecc), str(opi),
		str(stellar_teff), rad_mode,
		str(P0), str(hjd0), str(period), str(vsys), str(limph_T14), str(limph_T12),
		str(kp), str(ks), str(ecc), str(t14_hours), str(t12_hours),
		str(opi), str(stellar_teff), rad_mode,
		str(ra), str(dec),
		a_Rs_ratio=str(a_Rs_ratio),
		projected_obliquity=str(projected_obliquity),
		@@ -2067,7 +2077,7 @@ class Frame_Gofio:
		target_period, target_vsys, target_limph, target_kp, target_ecc, target_opi, target_stellar_teff, target_ra,
		target_dec,
		target_a_Rs_ratio, target_projected_obliquity, target_inclination,
		target_v_sini) = DataInterface.get_target_info(self.path_target, target_id)
		target_v_sini, target_limph_T12, target_t14_hours, target_t12_hours) = DataInterface.get_target_info(self.path_target, target_id)
		self.radius_DB_up.delete(0, tk.END)
		self.radius_DB_up.insert(tk.END, target_radius)
		self.mass_DB_up.delete(0, tk.END)