diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 405afb28..70da2154 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -12,7 +12,11 @@ jobs: python-version: ['3.10', '3.11', '3.12'] steps: - - uses: actions/checkout@v2 + - name: Checkout Project + uses: actions/checkout@v2 + + - name: Install CMake + uses: ssrobins/install-cmake@v1 - name: Setup Python ${{ matrix.python-version }} uses: actions/setup-python@v2 diff --git a/requirements.txt b/requirements.txt index eb86c4a0..a1b4653e 100644 --- a/requirements.txt +++ b/requirements.txt @@ -12,7 +12,6 @@ numba numpy < 2.0.0 pandas pooch -PyAstronomy PyMieScatt pymultinest requests diff --git a/species/data/companion_data/companion_data.json b/species/data/companion_data/companion_data.json index 2b09b11c..a0d97c29 100644 --- a/species/data/companion_data/companion_data.json +++ b/species/data/companion_data/companion_data.json @@ -1,6 +1,7 @@ { "AF Lep b": { "simbad": "V* AF Lep", + "gaia_dr3": 3009908378049913216, "parallax": [ 37.2539, 0.0195 @@ -46,6 +47,7 @@ }, "beta Pic b": { "simbad": "beta Pic", + "gaia_dr3": 4792774797545800832, "parallax": [ 50.9307, 0.1482 @@ -121,6 +123,7 @@ }, "beta Pic c": { "simbad": "beta Pic", + "gaia_dr3": 4792774797545800832, "parallax": [ 50.9307, 0.1482 @@ -159,6 +162,7 @@ }, "HIP 65426 b": { "simbad": "HIP 65426", + "gaia_dr3": 6070080754075553792, "parallax": [ 9.3031, 0.0346 @@ -223,6 +227,7 @@ }, "HIP 99770 b": { "simbad": "* b03 Cyg", + "gaia_dr3": 2060523518197844608, "parallax": [ 24.5456, 0.0911 @@ -254,6 +259,7 @@ }, "51 Eri b": { "simbad": "51 Eri", + "gaia_dr3": 3205095125321700480, "parallax": [ 33.439, 0.0777 @@ -323,6 +329,7 @@ }, "HR 8799 b": { "simbad": "HR 8799", + "gaia_dr3": 2832463659640297472, "parallax": [ 24.462, 0.0455 @@ -407,6 +414,7 @@ }, "HR 8799 c": { "simbad": "HR 8799", + "gaia_dr3": 2832463659640297472, "parallax": [ 24.462, 0.0455 @@ -486,6 +494,7 @@ }, "HR 8799 d": { "simbad": "HR 8799", + "gaia_dr3": 2832463659640297472, "parallax": [ 24.462, 0.0455 @@ -565,6 +574,7 @@ }, "HR 8799 e": { "simbad": "HR 8799", + "gaia_dr3": 2832463659640297472, "parallax": [ 24.462, 0.0455 @@ -634,6 +644,7 @@ }, "HD 95086 b": { "simbad": "HD 95086", + "gaia_dr3": 5231963962676292224, "parallax": [ 11.5659, 0.0187 @@ -681,6 +692,7 @@ }, "PDS 70 b": { "simbad": "PDS 70", + "gaia_dr3": 6110141563309613056, "parallax": [ 8.8975, 0.0191 @@ -768,6 +780,7 @@ }, "PDS 70 c": { "simbad": "PDS 70", + "gaia_dr3": 6110141563309613056, "parallax": [ 8.8975, 0.0191 @@ -812,6 +825,7 @@ }, "2M 1207 B": { "simbad": "2MASSWJ 1207334-393254", + "gaia_dr3": 3459372646830687104, "parallax": [ 15.4624, 0.1163 @@ -880,6 +894,7 @@ }, "AB Pic B": { "simbad": "AB Pic", + "gaia_dr3": 5495052596695570816, "parallax": [ 19.9452, 0.0124 @@ -920,6 +935,7 @@ }, "HD 206893 B": { "simbad": "HD 206893", + "gaia_dr3": 6843672087120107264, "parallax": [ 24.5275, 0.0354 @@ -982,6 +998,7 @@ }, "HD 206893 c": { "simbad": "HD 206893", + "gaia_dr3": 6843672087120107264, "parallax": [ 24.5275, 0.0354 @@ -1014,6 +1031,7 @@ }, "RZ Psc B": { "simbad": "RZ Psc", + "gaia_dr3": 307981178401165184, "parallax": [ 5.3661, 0.0274 @@ -1055,6 +1073,7 @@ }, "GQ Lup B": { "simbad": "GQ Lup", + "gaia_dr3": 6011522757643074304, "parallax": [ 6.4893, 0.0289 @@ -1180,6 +1199,7 @@ }, "PZ Tel B": { "simbad": "PZ Tel", + "gaia_dr3": 6655168686921108864, "parallax": [ 21.1621, 0.0223 @@ -1274,6 +1294,7 @@ }, "kappa And b": { "simbad": "kappa And", + "gaia_dr3": 1926476042681491712, "parallax": [ 19.4064, 0.2104 @@ -1327,6 +1348,7 @@ }, "HD 1160 B": { "simbad": "HD 1160", + "gaia_dr3": 2741090498161113344, "parallax": [ 8.2721, 0.0355 @@ -1375,7 +1397,8 @@ ] }, "ROXs 12 B": { - "simbad": "ROXs 12 B", + "simbad": "ROXs 12", + "gaia_dr3": 6048935358761628288, "parallax": [ 7.217, 0.0172 @@ -1416,6 +1439,7 @@ }, "ROXs 42 Bb": { "simbad": "ROXs 42 B", + "gaia_dr3": 6047587937321868288, "parallax": [ 6.8284, 0.0309 @@ -1471,6 +1495,7 @@ }, "GJ 504 b": { "simbad": "GJ 504", + "gaia_dr3": 3732539683617410816, "parallax": [ 56.8577, 0.1224 @@ -1544,41 +1569,42 @@ ] }, "GJ 758 B": { - "simbad": "GJ 758", - "distance": [15.61, 0.005], - "parallax": [64.04, 0.02], - "app_mag": { - "Paranal/SPHERE.IRDIS_D_Y23_2": [20.15, 0.20], - "Paranal/SPHERE.IRDIS_D_Y23_3": [19.39, 0.10], - "Paranal/SPHERE.IRDIS_D_J23_2": [20.02, 0.25], - "Paranal/SPHERE.IRDIS_D_J23_3": [17.79, 0.18], - "Paranal/SPHERE.IRDIS_D_H23_2": [17.55, 0.12], - "Paranal/SPHERE.IRDIS_D_H23_3": [19.84, 0.42], - "Paranal/SPHERE.IRDIS_D_K12_1": [17.99, 0.21], - "Paranal/SPHERE.IRDIS_D_K12_2": [18.74, 0.35], - "Subaru/CIAO.J": [17.58, 0.20], - "Gemini/NIRI.H-G0203w": [18.16, 0.2], - "Gemini/NIRI.Kcont209-G0217": [17.12, 0.2], - "Keck/NIRC2.Lp": [15.0, 0.1], - "Gemini/NIRI.CH4short-G0228": [17.74, 0.2] - }, - "mass_star": [0.96, 0.03], - "mass_companion": [38.0, 0.8], - "semi_major": [ - 29.7, - -5.3, - 5.3 - ], - "accretion": false, - "references": [ - "Gaia Early Data Release 3", - "Vigan et al. 2016", - "Janson et al. 2011", - "Brandt, G. M. et al. 2021" - ] + "simbad": "GJ 758", + "gaia_dr3": 2046207670631964928, + "parallax": [64.04, 0.02], + "app_mag": { + "Paranal/SPHERE.IRDIS_D_Y23_2": [20.15, 0.20], + "Paranal/SPHERE.IRDIS_D_Y23_3": [19.39, 0.10], + "Paranal/SPHERE.IRDIS_D_J23_2": [20.02, 0.25], + "Paranal/SPHERE.IRDIS_D_J23_3": [17.79, 0.18], + "Paranal/SPHERE.IRDIS_D_H23_2": [17.55, 0.12], + "Paranal/SPHERE.IRDIS_D_H23_3": [19.84, 0.42], + "Paranal/SPHERE.IRDIS_D_K12_1": [17.99, 0.21], + "Paranal/SPHERE.IRDIS_D_K12_2": [18.74, 0.35], + "Subaru/CIAO.J": [17.58, 0.20], + "Gemini/NIRI.H-G0203w": [18.16, 0.2], + "Gemini/NIRI.Kcont209-G0217": [17.12, 0.2], + "Keck/NIRC2.Lp": [15.0, 0.1], + "Gemini/NIRI.CH4short-G0228": [17.74, 0.2] + }, + "mass_star": [0.96, 0.03], + "mass_companion": [38.0, 0.8], + "semi_major": [ + 29.7, + -5.3, + 5.3 + ], + "accretion": false, + "references": [ + "Gaia Early Data Release 3", + "Vigan et al. 2016", + "Janson et al. 2011", + "Brandt, G. M. et al. 2021" + ] }, "GU Psc b": { "simbad": "GU Psc", + "gaia_dr3": 2784120641627589760, "parallax": [ 21.0051, 0.026 @@ -1634,6 +1660,7 @@ }, "2M0103 ABb": { "simbad": "2MASS J01033563-5515561", + "gaia_dr3": 4914632365579875200, "parallax": [ 21.18, 1.37 @@ -1688,6 +1715,7 @@ }, "1RXS 1609 B": { "simbad": "1RXS J160929.1-210524", + "gaia_dr3": 6243841249531772800, "parallax": [ 7.2444, 0.0176 @@ -1734,6 +1762,7 @@ }, "GSC 06214 B": { "simbad": "GSC 06214-00210", + "gaia_dr3": 6244440552088537600, "parallax": [ 9.1923, 0.0214 @@ -1794,6 +1823,7 @@ }, "HD 72946 B": { "simbad": "HD 72946", + "gaia_dr3": 594989417312890368, "parallax": [ 38.9809, 0.0412 @@ -1829,6 +1859,7 @@ }, "HIP 64892 B": { "simbad": "HIP 64892", + "gaia_dr3": 6088264477373916032, "parallax": [ 8.3595, 0.0483 @@ -1876,6 +1907,7 @@ }, "HD 13724 B": { "simbad": "HD 13724", + "gaia_dr3": 4942867480584948352, "parallax": [ 23.0157, 0.0178 @@ -1927,6 +1959,7 @@ }, "YSES 1 b": { "simbad": "TYC 8998-760-1", + "gaia_dr3": 5864061893213196032, "parallax": [ 10.6124, 0.0116 @@ -2007,6 +2040,7 @@ }, "YSES 1 c": { "simbad": "TYC 8998-760-1", + "gaia_dr3": 5864061893213196032, "parallax": [ 10.6124, 0.0116 @@ -2071,6 +2105,7 @@ }, "YSES 2 b": { "simbad": "TYC 8984-2245-1", + "gaia_dr3": 5236792880333011968, "parallax": [ 9.1115, 0.031 @@ -2106,6 +2141,7 @@ }, "HD 142527 B": { "simbad": "HD 142527", + "gaia_dr3": 5994826707951507200, "parallax": [ 6.2791, 0.0284 @@ -2225,6 +2261,7 @@ }, "CS Cha B": { "simbad": "CS Cha", + "gaia_dr3": 5201179389431267584, "parallax": [ 5.9251, 0.0675 @@ -2271,6 +2308,7 @@ }, "CT Cha B": { "simbad": "CT Cha", + "gaia_dr3": 5201360671411974912, "parallax": [ 5.2645, 0.0116 @@ -2315,6 +2353,7 @@ }, "SR 12 C": { "simbad": "2MASS J16271951-2441403", + "gaia_dr3": 6049072007438134272, "parallax": [ 8.9034, 0.4288 @@ -2371,6 +2410,7 @@ }, "DH Tau B": { "simbad": "DH Tau", + "gaia_dr3": 151374202498079872, "parallax": [ 7.4936, 0.0255 @@ -2412,6 +2452,7 @@ }, "HD 4747 B": { "simbad": "HD 4747", + "gaia_dr3": 2348830516542653824, "parallax": [ 53.0526, 0.0282 @@ -2448,6 +2489,7 @@ }, "HR 3549 B": { "simbad": "HR 3549", + "gaia_dr3": 5304593027881569024, "parallax": [ 10.5509, 0.0383 @@ -2502,6 +2544,7 @@ }, "CHXR 73 B": { "simbad": "CHXR 73", + "gaia_dr3": 5201175987817179136, "parallax": [ 5.2256, 0.0804 @@ -2537,6 +2580,7 @@ }, "HD 19467 B": { "simbad": "HD 19467", + "gaia_dr3": 5155981730587168384, "parallax": [ 31.2191, 0.024 @@ -2584,6 +2628,7 @@ }, "b Cen (AB)b": { "simbad": "b Cen", + "gaia_dr3": 6105713692544123520, "parallax": [ 10.0339, 0.3142 @@ -2627,6 +2672,7 @@ }, "VHS 1256 B": { "simbad": "SIPS J1256-1257", + "gaia_dr3": 3526198184723289472, "parallax": [ 47.2733, 0.4731 diff --git a/species/data/database.py b/species/data/database.py index 5721e7a9..f6646222 100644 --- a/species/data/database.py +++ b/species/data/database.py @@ -821,7 +821,7 @@ def add_custom_model( ``logkzz`` (for :math:`\\log\\,K_\\mathrm{zz}`), ``adindex`` (for :math:`\\gamma_\\mathrm{ad}`), and ``log_co_iso`` (for - :math:`\\log\,^{12}\\mathrm{CO}/^{13}\\mathrm{CO}`). + :math:`\\log\\,^{12}\\mathrm{CO}/^{13}\\mathrm{CO}`). Please contact the code maintainer if support for other parameters should be added. wavel_range : tuple(float, float), None @@ -2419,7 +2419,7 @@ def get_probable_sample( print("\nParameters:") for param_key, param_value in prob_sample.items(): if isinstance(param_value, float): - if param_key in ["luminosity", "flux_scaling", "flux_offset"]: + if param_value < 0.01: print(f" - {param_key} = {param_value:.2e}") else: print(f" - {param_key} = {param_value:.2f}") @@ -2514,7 +2514,7 @@ def get_median_sample( print("\nParameters:") for param_key, param_value in median_sample.items(): if isinstance(param_value, float): - if param_key in ["luminosity", "flux_scaling", "flux_offset"]: + if param_value < 0.01: print(f" - {param_key} = {param_value:.2e}") else: print(f" - {param_key} = {param_value:.2f}") @@ -2581,11 +2581,11 @@ def get_compare_sample(self, tag: str, verbose: bool = True) -> Dict[str, float] if verbose: print("\nParameters:") - for key, value in model_param.items(): - if key in ["luminosity", "flux_scaling", "flux_offset"]: - print(f" - {key} = {value:.2e}") + for param_key, param_value in model_param.items(): + if param_value < 0.01: + print(f" - {param_key} = {param_value:.2e}") else: - print(f" - {key} = {value:.2f}") + print(f" - {param_key} = {param_value:.2f}") return model_param diff --git a/species/data/model_data/model_data.json b/species/data/model_data/model_data.json index 73fbe8df..565110e8 100644 --- a/species/data/model_data/model_data.json +++ b/species/data/model_data/model_data.json @@ -176,29 +176,29 @@ "url": "https://ui.adsabs.harvard.edu/abs/2008ApJ...675L.105H", "checksum": "b59b7711a4ed9dcf9284ef0c3f2b614598eec89bf53f548ceb361b806db70757" }, - "exo-rem-cloud-r500": { + "exo-rem-lowres": { "parameters": ["teff", "logg", "feh", "c_o_ratio"], "name": "Exo-REM", - "file size": "703 MB", - "wavelength range": [0.35, 250.0], + "file size": "493 MB", + "wavelength range": [0.34, 250.0], "lambda/d_lambda": 500, "teff range": [400, 2000], "reference": "Charnay et al. (2018)", "url": "https://ui.adsabs.harvard.edu/abs/2018ApJ...854..172C", - "checksum": "d91ca3a8bcde40c744ce7d9712ae74be822c2d4b394ce3394e89b8ef6ad79a74" + "checksum": "58df9a1227286644a3c4c505d573228c1ec6f10dd652ecbbee4f875c3d8ca86a" }, - "exo-rem-nocloud-r500": { + "exo-rem-nocloud": { "parameters": ["teff", "logg", "feh", "c_o_ratio"], "name": "Exo-REM", - "file size": "508 MB", + "file size": "356 MB", "wavelength range": [0.35, 250.0], "lambda/d_lambda": 500, "teff range": [400, 2000], "reference": "Charnay et al. (2018)", "url": "https://ui.adsabs.harvard.edu/abs/2018ApJ...854..172C", - "checksum": "af38e6a7243f4e11712363f776efedffa57055c22e471509a6276020ec70fc78" + "checksum": "52addf881edbeba5f2ced0628f30325c3890e2b7dfcdc7ec000565c93fdd0e60" }, - "exo-rem-cloud-r20000": { + "exo-rem-highres": { "parameters": ["teff", "logg", "feh", "c_o_ratio"], "name": "Exo-REM", "file size": "21 GB", @@ -209,7 +209,7 @@ "url": "https://ui.adsabs.harvard.edu/abs/2018ApJ...854..172C", "checksum": "4dd73d2e5a12756c05b9439cf7f4c03d2cc084e33b97af0efc3da3409c83a092" }, - "exo-rem-cloud-r20000-isotopes": { + "exo-rem-highres-isotopes": { "parameters": ["teff", "logg", "feh", "c_o_ratio", "log_co_iso"], "name": "Exo-REM", "file size": "239 MB", diff --git a/species/fit/fit_evolution.py b/species/fit/fit_evolution.py index 530b267b..f56ea876 100644 --- a/species/fit/fit_evolution.py +++ b/species/fit/fit_evolution.py @@ -481,18 +481,18 @@ def ln_prior(cube, n_dim, n_param) -> None: a_trunc = ( self.bounds["age"][0] - self.age_prior[0] - ) / np.abs(self.age_prior[2]) + ) / self.age_prior[2] b_trunc = ( self.bounds["age"][1] - self.age_prior[0] - ) / np.abs(self.age_prior[2]) + ) / self.age_prior[2] cube[cube_index["age"]] = truncnorm.ppf( cube[cube_index["age"]], a_trunc, b_trunc, loc=self.age_prior[0], - scale=np.abs(self.age_prior[2]), + scale=self.age_prior[2], ) else: diff --git a/species/fit/fit_model.py b/species/fit/fit_model.py index 11416b8d..c7d0cd36 100644 --- a/species/fit/fit_model.py +++ b/species/fit/fit_model.py @@ -143,14 +143,14 @@ def __init__( - Radial velocity can be included with the ``rad_vel`` parameter (km/s). This parameter will only be relevant - if the radial velocity shift can be spectrally - resolved given the instrument resolution. When - including ``rad_vel``, a single RV will be fitted - for all spectra. Or, it is also possible by fitting - an RV for individual spectra by for example including - the parameter as ``rad_vel_SPHERE`` in case the - spectrum name is ``SPHERE``, that is, the name used - as tag when adding the spectrum to the database with + if the radial velocity shift can be detected given + the instrument resolution. When including ``rad_vel``, + a single RV will be fitted for all spectra. Or, it is + also possible by fitting an RV for individual spectra + by for example including the parameter as + ``rad_vel_CRIRES`` in case the spectrum name is + ``CRIRES``, that is, the name used as tag when adding + the spectrum to the database with :func:`~species.data.database.Database.add_object`. - Rotational broadening can be fitted by including the @@ -161,23 +161,15 @@ def __init__( resolution. The resolution is set when adding a spectrum to the database with :func:`~species.data.database.Database.add_object`. - Note that the broadening is applied with the - `fastRotBroad `_ function from ``PyAstronomy``. - The rotational broadening is only accurate if the - wavelength range of the data is somewhat narrow. - For example, when fitting a medium- or - high-resolution spectrum across multiple bands - (e.g. $JHK$ bands) then it is best to split up the - data into the separate bands when adding them with - :func:`~species.data.database.Database.add_object`. + The broadening is applied with the function from + `Carvalho & Johns-Krull (2023) `_. A single broadening parameter, ``vsini``, can be fitted, so it is applied for all spectra. Or, it is also possible to fit the broadening for individual spectra by for example including the parameter as - ``vsini_SPHERE`` in case the spectrum name is - ``SPHERE``, that is, the name used as tag when + ``vsini_CRIRES`` in case the spectrum name is + ``CRIRES``, that is, the name used as tag when adding the spectrum to the database with :func:`~species.data.database.Database.add_object`. @@ -601,7 +593,7 @@ def __init__( "The use of 'ism_ext' for fitting the extinction is " "deprecated. Please use the 'ext_model' parameter of " "'FitModel' in combination with setting 'ext_av' in " - "the 'model_param' dictionary", + "the 'bounds' dictionary.", DeprecationWarning, ) @@ -1202,6 +1194,7 @@ def __init__( if bounds[spec_item][0] is not None: # Add the flux scaling parameter self.modelpar.append(f"scaling_{spec_item}") + self.bounds[f"scaling_{spec_item}"] = ( bounds[spec_item][0][0], bounds[spec_item][0][1], @@ -1210,6 +1203,7 @@ def __init__( if len(bounds[spec_item]) > 1 and bounds[spec_item][1] is not None: # Add the error inflation parameters self.modelpar.append(f"error_{spec_item}") + self.bounds[f"error_{spec_item}"] = ( bounds[spec_item][1][0], bounds[spec_item][1][1], @@ -2207,7 +2201,7 @@ def _lnlike_func( # # model_flux = all_param["veil_a"] * model_flux + veil_flux - # Optionally scale the data to account for calibration + # Optional flux scaling to account for calibration inaccuracy if f"scaling_{spec_item}" in all_param: spec_scaling = all_param[f"scaling_{spec_item}"] @@ -2215,35 +2209,41 @@ def _lnlike_func( spec_scaling = 1.0 data_flux = spec_scaling * self.spectrum[spec_item][0][:, 1] + data_var = spec_scaling**2 * self.spectrum[spec_item][0][:, 2] ** 2 - # Optionally inflate the data uncertainties + # Optional variance inflation (see Piette & Madhusudhan 2020) if f"error_{spec_item}" in all_param: - # Variance with error inflation (see Piette & Madhusudhan 2020) - data_var = ( - self.spectrum[spec_item][0][:, 2] ** 2 - + (all_param[f"error_{spec_item}"] * model_flux) ** 2 - ) - else: - # Variance without error inflation - data_var = self.spectrum[spec_item][0][:, 2] ** 2 + data_var += (all_param[f"error_{spec_item}"] * model_flux) ** 2 # Select the inverted covariance matrix if self.spectrum[spec_item][2] is not None: if f"error_{spec_item}" in all_param: # Ratio of the inflated and original uncertainties - sigma_ratio = np.sqrt(data_var) / self.spectrum[spec_item][0][:, 2] + sigma_ratio = np.sqrt(data_var) / ( + spec_scaling * self.spectrum[spec_item][0][:, 2] + ) sigma_j, sigma_i = np.meshgrid(sigma_ratio, sigma_ratio) # Calculate the inverted matrix of the inflated covariances data_cov_inv = np.linalg.inv( - self.spectrum[spec_item][1] * sigma_i * sigma_j + spec_scaling**2 + * sigma_i + * sigma_j + * self.spectrum[spec_item][1] ) else: - # Use the inverted covariance matrix directly - data_cov_inv = self.spectrum[spec_item][2] + if spec_scaling == 1.0: + # Use the inverted covariance matrix directly + data_cov_inv = self.spectrum[spec_item][2] + else: + # Apply flux scaling to covariance matrix + # and then invert the covariance matrix + data_cov_inv = np.linalg.inv( + spec_scaling**2 * self.spectrum[spec_item][1] + ) # Calculate the log-likelihood diff --git a/species/fit/retrieval.py b/species/fit/retrieval.py index b28a0ca5..c65cd889 100644 --- a/species/fit/retrieval.py +++ b/species/fit/retrieval.py @@ -33,7 +33,6 @@ ) from molmass import Formula -from PyAstronomy.pyasl import fastRotBroad from schwimmbad import MPIPool from scipy.integrate import simpson from scipy.interpolate import interp1d @@ -44,9 +43,10 @@ from species.phot.syn_phot import SyntheticPhotometry from species.read.read_filter import ReadFilter from species.read.read_object import ReadObject +from species.util.convert_util import logg_to_mass from species.util.core_util import print_section from species.util.dust_util import apply_ism_ext -from species.util.convert_util import logg_to_mass +from species.util.model_util import rot_int_cmj from species.util.retrieval_util import ( calc_metal_ratio, calc_spectrum_clear, @@ -2860,26 +2860,13 @@ def _lnlike( # Apply optional rotational broadening if "vsini" in cube_index and spec_item in self.apply_vsini: - spec_interp = interp1d(model_wavel, model_flux) - - wavel_new = np.linspace( - model_wavel[0], - model_wavel[-1], - 2 * model_wavel.shape[0], - ) - - flux_broad = fastRotBroad( - wvl=wavel_new, - flux=spec_interp(wavel_new), - epsilon=0.0, + model_flux = rot_int_cmj( + wavel=model_wavel, + flux=model_flux, vsini=cube[cube_index["vsini"]], - effWvl=None, + eps=0.0, ) - spec_interp = interp1d(wavel_new, flux_broad) - - model_flux = spec_interp(model_wavel) - # Shift the wavelengths of the data with # the fitted calibration parameter data_wavel = self.spectrum[spec_item][0][:, 0] + wavel_cal[spec_item] @@ -3209,23 +3196,15 @@ def setup_retrieval( - Rotational broadening can be fitted by including the ``vsini`` parameter (km/s). This parameter will only - be relevant if the rotational broadening is stronger - than or comparable to the instrumental broadening, + be relevant if the rotational broadening is + significant compared to the instrumental broadening, so typically when the data has a high spectral resolution. The resolution is set when adding a spectrum to the database with :func:`~species.data.database.Database.add_object`. - Note that the broadening is applied with the - `fastRotBroad `_ function from ``PyAstronomy``. - The rotational broadening is only accurate if the - wavelength range of the data is somewhat narrow. - For example, when fitting a medium- or - high-resolution spectrum across multiple bands - (e.g. $JHK$ bands) then it is best to split up the - data into the separate bands when adding them with - :func:`~species.data.database.Database.add_object`. + The broadening is applied with the function from + `Carvalho & Johns-Krull (2023) `_. Calibration parameters (optional): @@ -3390,7 +3369,7 @@ def setup_retrieval( the spectra (i.e. excluding low-resolution spectra for which the broadening will not matter), the computation will be a bit faster. This parameter is only used when - the ``vsini`` model parameter has been include in + the ``vsini`` model parameter has been included in ``bounds``. The :math:`v \\sin(i)` is applied to all spectra by setting the argument of ``apply_vsini`` to ``None``. diff --git a/species/read/read_model.py b/species/read/read_model.py index a2078891..539912bc 100644 --- a/species/read/read_model.py +++ b/species/read/read_model.py @@ -13,7 +13,6 @@ import numpy as np from astropy import units as u -from PyAstronomy.pyasl import rotBroad, fastRotBroad from typeguard import typechecked from scipy.integrate import simpson from scipy.interpolate import interp1d, RegularGridInterpolator @@ -34,7 +33,7 @@ from species.read.read_planck import ReadPlanck from species.util.convert_util import logg_to_mass from species.util.dust_util import check_dust_database, ism_extinction, convert_to_av -from species.util.model_util import binary_to_single +from species.util.model_util import binary_to_single, rot_int_cmj from species.util.spec_util import smooth_spectrum @@ -599,7 +598,6 @@ def get_model( spec_res: Optional[float] = None, wavel_resample: Optional[np.ndarray] = None, magnitude: bool = False, - fast_rot_broad: bool = True, ext_filter: Optional[str] = None, **kwargs, ) -> ModelBox: @@ -628,14 +626,6 @@ def get_model( magnitude : bool Normalize the spectrum with a flux calibrated spectrum of Vega and return the magnitude instead of flux density. - fast_rot_broad : bool - Apply fast algorithm for the rotational broadening if set - to ``True``, otherwise a slow but more accurate broadening - is applied if set to ``False``. The fast algorithm will - only provide an accurate broadening if the wavelength range - of the spectrum is somewhat narrow (e.g. only the :math:`K` - band). The argument is only used if the ``vsini`` parameter - is included in the ``model_param`` dictionary. ext_filter : str, None Filter that is associated with the (optional) extinction parameter, ``ism_ext``. When the argument of ``ext_filter`` @@ -879,51 +869,11 @@ def get_model( # Apply rotational broadening vsin(i) in km/s if "vsini" in model_param: - # fastRotBroad requires a linear wavelength sampling - # while pRT uses a logarithmic wavelength sampling, - # so change temporarily to a linear sampling - # with a factor 5 larger number of wavelengths - - wavel_linear = np.linspace( - model_box.wavelength[0], - model_box.wavelength[-1], - 5 * model_box.wavelength.size, - ) - - flux_linear = spectres_numba( - wavel_linear, - model_box.wavelength, - model_box.flux, - spec_errs=None, - fill=np.nan, - verbose=True, - ) - - if fast_rot_broad: - flux_broad = fastRotBroad( - wvl=wavel_linear, - flux=flux_linear, - epsilon=0.0, - vsini=model_param["vsini"], - effWvl=None, - ) - - else: - flux_broad = rotBroad( - wvl=wavel_linear, - flux=flux_linear, - epsilon=0.0, - vsini=model_param["vsini"], - edgeHandling="firstlast", - ) - - model_box.flux = spectres_numba( - model_box.wavelength, - wavel_linear, - flux_broad, - spec_errs=None, - fill=np.nan, - verbose=True, + model_box.flux = rot_int_cmj( + wavel=model_box.wavelength, + flux=model_box.flux, + vsini=model_param["vsini"], + eps=0.0, ) # Apply veiling @@ -1398,6 +1348,31 @@ def get_data( spec_res=spec_res, ) + # Apply rotational broadening vsin(i) in km/s + + if "vsini" in model_param: + model_box.flux = rot_int_cmj( + wavel=model_box.wavelength, + flux=model_box.flux, + vsini=model_param["vsini"], + eps=0.0, + ) + + # Apply veiling + + if ( + "veil_a" in model_param + and "veil_b" in model_param + and "veil_ref" in model_param + ): + lambda_ref = 0.5 # (um) + + veil_flux = model_param["veil_ref"] + model_param["veil_b"] * ( + model_box.wavelength - lambda_ref + ) + + model_box.flux = model_param["veil_a"] * model_box.flux + veil_flux + # Apply extinction if ( @@ -2007,25 +1982,26 @@ def integrate_spectrum(self, model_param: Dict[str, float]) -> float: or self.wavel_range[1] != wavel_points[-1] ): warnings.warn( - "The 'wavel_range' is not set to the " - "maximum available range. To maximize the " - "accuracy when calculating the bolometric " - "luminosity, it is recommended to set " - "'wavel_range=None'." + "The 'wavel_range' is not set to the maximum " + "available range. To maximize the accuracy when " + "calculating the bolometric luminosity, it is " + "recommended to set 'wavel_range=None'." ) - if "parallax" in model_param: - del model_param["parallax"] + param_copy = model_param.copy() - if "distance" in model_param: - del model_param["distance"] + if "parallax" in param_copy: + del param_copy["parallax"] - model_box = self.get_model(model_param) + if "distance" in param_copy: + del param_copy["distance"] + + model_box = self.get_model(param_copy) bol_lum = ( 4.0 * np.pi - * (model_param["radius"] * constants.R_JUP) ** 2 + * (param_copy["radius"] * constants.R_JUP) ** 2 * simpson(y=model_box.flux, x=model_box.wavelength) ) diff --git a/species/read/read_radtrans.py b/species/read/read_radtrans.py index 510113dd..2d0334c7 100644 --- a/species/read/read_radtrans.py +++ b/species/read/read_radtrans.py @@ -15,7 +15,6 @@ import spectres from matplotlib.ticker import MultipleLocator -from PyAstronomy.pyasl import fastRotBroad from scipy.interpolate import interp1d from spectres.spectral_resampling_numba import spectres_numba from typeguard import typechecked @@ -26,6 +25,7 @@ from species.read.read_filter import ReadFilter from species.util.convert_util import logg_to_mass from species.util.dust_util import apply_ism_ext +from species.util.model_util import rot_int_cmj from species.util.retrieval_util import ( calc_metal_ratio, calc_spectrum_clear, @@ -445,11 +445,9 @@ def get_model( - Rotational broadening can be applied by adding the ``vsini`` parameter, which is the projected spin velocity (km/s), :math:`v\\sin{i}`. The broadening - is applied with the ``fastRotBroad`` function from - ``PyAstronomy`` (see for details the `documentation - `_). + is applied with the function from `Carvalho & + Johns-Krull (2023) `_. quenching : str, None Quenching type for CO/CH$_4$/H$_2$O abundances. Either @@ -1163,46 +1161,11 @@ def get_model( # Convolve with a broadening kernel for vsin(i) if "vsini" in model_param: - # fastRotBroad requires a linear wavelength sampling - # while pRT uses a logarithmic wavelength sampling - # so change temporarily to a linear sampling - # with a factor 10 larger number of wavelengths - - wavel_linear = np.linspace( - np.amin(wavelength), np.amax(wavelength), wavelength.size * 10 - ) - - flux_linear = spectres_numba( - wavel_linear, - wavelength, - flux, - spec_errs=None, - fill=np.nan, - verbose=True, - ) - - # Apply the rotational broadening - # The rotBroad function is much slower than - # fastRotBroad when tested on a large array - - flux_broad = fastRotBroad( - wvl=wavel_linear, - flux=flux_linear, - epsilon=1.0, + flux = rot_int_cmj( + wavel=wavelength, + flux=flux, vsini=model_param["vsini"], - effWvl=None, - ) - - # And change back to the original (logarithmic) - # wavelength sampling, with constant R - - flux = spectres_numba( - wavelength, - wavel_linear, - flux_broad, - spec_errs=None, - fill=np.nan, - verbose=True, + eps=0.0, ) # Convolve the spectrum with a Gaussian LSF diff --git a/species/util/box_util.py b/species/util/box_util.py index 64dbabfa..cd0a0754 100644 --- a/species/util/box_util.py +++ b/species/util/box_util.py @@ -94,15 +94,25 @@ def update_objectbox( spec_tmp = value[0] if f"scaling_{key}" in model_param: - # Scale the flux of the spectrum + # Scale the fluxes and uncertainties of the spectrum scaling = model_param[f"scaling_{key}"] - print( - f"Scaling the flux of {key} by: {scaling:.2f}...", - end="", - flush=True, - ) + if scaling < 0.01: + print( + f"Scaling the flux of {key} by: {scaling:.2e}...", + end="", + flush=True, + ) + + else: + print( + f"Scaling the flux of {key} by: {scaling:.2f}...", + end="", + flush=True, + ) + spec_tmp[:, 1] *= model_param[f"scaling_{key}"] + spec_tmp[:, 2] *= model_param[f"scaling_{key}"] print(" [DONE]") if f"error_{key}" in model_param: diff --git a/species/util/data_util.py b/species/util/data_util.py index c9b2dcb5..37f03323 100644 --- a/species/util/data_util.py +++ b/species/util/data_util.py @@ -416,7 +416,7 @@ def add_missing( count_interp = 0 count_missing = 0 - if np.isinf(np.sum(flux)): + if not np.isfinite(np.sum(flux)): print("\nFix missing grid points with a linear interpolation:") if len(parameters) == 1: diff --git a/species/util/fit_util.py b/species/util/fit_util.py index 0ff43891..0c0a21e4 100644 --- a/species/util/fit_util.py +++ b/species/util/fit_util.py @@ -236,7 +236,7 @@ def get_residuals( objectbox : ObjectBox Box with the photometry and/or spectra of an object. A scaling and/or error inflation of the spectra should be applied with - :func:`~species.util.read_util.update_objectbox` beforehand. + :func:`~species.util.box_util.update_objectbox` beforehand. inc_phot : bool, list(str) Include photometric data in the fit. If a boolean, either all (``True``) or none (``False``) of the data are selected. If a diff --git a/species/util/model_util.py b/species/util/model_util.py index a5dceb7b..8ba956cc 100644 --- a/species/util/model_util.py +++ b/species/util/model_util.py @@ -12,7 +12,8 @@ import numpy as np from astropy import units as u -from PyAstronomy.pyasl import fastRotBroad + +# from PyAstronomy.pyasl import fastRotBroad from scipy.interpolate import RegularGridInterpolator from spectres.spectral_resampling_numba import spectres_numba from typeguard import typechecked @@ -344,44 +345,11 @@ def apply_obs( # Apply rotational broadening if rot_broad is not None: - # The fastRotBroad requires constant wavelength steps - # Upsample by a factor of 4 to not lose spectral information - - wavel_new = np.linspace( - model_wavel[0], - model_wavel[-1], - 4 * model_wavel.size, - ) - - flux_new = spectres_numba( - wavel_new, - model_wavel, - model_flux, - spec_errs=None, - fill=np.nan, - verbose=True, - ) - - # Apply fast rotational broadening - # Only to be used on a limited wavelength range - - flux_broad = fastRotBroad( - wvl=wavel_new, - flux=flux_new, - epsilon=0.0, + model_flux = rot_int_cmj( + wavel=model_wavel, + flux=model_flux, vsini=rot_broad, - effWvl=None, - ) - - # Interpolate back to the original wavelength sampling - - model_flux = spectres_numba( - model_wavel, - wavel_new, - flux_broad, - spec_errs=None, - fill=np.nan, - verbose=True, + eps=0.0, ) # Apply radial velocity shift @@ -518,8 +486,8 @@ def apply_obs( if data_wavel is not None: # The 'fill' should not happen because model_wavel is # 20 wavelength points broader than data_wavel, but - # spectres sometimes set the outer fluxes to 'fill' - # when a smaller margin than 20 wavelengths was used. + # spectres sometimes sets the outer fluxes to 'fill' + # depending on the spectral resolution of the data model_flux = spectres_numba( data_wavel, @@ -536,3 +504,87 @@ def apply_obs( # model_flux += model_param["flux_offset"] return model_flux + + +@typechecked +def rot_int_cmj( + wavel: np.ndarray, + flux: np.ndarray, + vsini: float, + eps: float = 0.6, + nr: int = 10, + ntheta: int = 100, + dif: float = 0.0, +): + """ + A routine to quickly rotationally broaden a spectrum in linear time. + This function has been adopted from `Carvalho & Johns-Krull (2023) + `_. + + Parameters + ---------- + wavel : np.ndarray + Array with the wavelengths. + flux : np.ndarray + Array with the fluxes + vsini : float + Projected rotational velocity (km s-1). + eps : float + Coefficient of the limb darkening law (default: 0.0). + nr : int + Number of radial bins on the projected disk (default: 10). + ntheta : int + Number of azimuthal bins in the largest radial annulus + (default: 100). Note: the number of bins at each r is + int(r*ntheta) where r < 1. + dif : float + Differential rotation coefficient (default = 0.0), applied + according to the law Omeg(th)/Omeg(eq) = (1 - dif/2 - (dif/2) + cos(2 th)). Dif = 0.675 reproduces the law proposed by Smith, + 1994, A&A, Vol. 287, p. 523-534, to unify WTTS and CTTS. + Dif = 0.23 is similar to observed solar differential rotation. + Note: the th in the expression above is the stellar co-latitude, + which is not the same as the integration variable used in the + function. This is a disk integration routine. + + Returns + ------- + np.ndarray + Array with the rotationally broadened spectrum. + """ + + ns = np.zeros(flux.shape) + tarea = 0.0 + dr = 1.0 / nr + + for j in range(nr): + r = dr / 2.0 + j * dr + area = ( + ((r + dr / 2.0) ** 2 - (r - dr / 2.0) ** 2) + / int(ntheta * r) + * (1.0 - eps + eps * np.cos(np.arcsin(r))) + ) + + for k in range(int(ntheta * r)): + th = np.pi / int(ntheta * r) + k * 2.0 * np.pi / int(ntheta * r) + + if dif != 0: + vl = ( + vsini + * r + * np.sin(th) + * ( + 1.0 + - dif / 2.0 + - dif / 2.0 * np.cos(2.0 * np.arccos(r * np.cos(th))) + ) + ) + else: + vl = r * vsini * np.sin(th) + + ns += area * np.interp( + wavel + wavel * vl / (1e-3 * constants.LIGHT), wavel, flux + ) + tarea += area + + return ns / tarea