From b4786d517aba703c67a9b9f13cfb15561b9c3a9a Mon Sep 17 00:00:00 2001 From: Kyle Conroy Date: Fri, 4 Oct 2024 11:14:26 -0400 Subject: [PATCH 01/10] hacky support for contact temperature smoothing --- phoebe/backend/backends.py | 35 +++- phoebe/backend/contacts_smoothing.py | 174 ++++++++++++++++++ phoebe/backend/universe.py | 257 +++++++++++++++++++++++++++ phoebe/parameters/component.py | 2 + 4 files changed, 467 insertions(+), 1 deletion(-) create mode 100644 phoebe/backend/contacts_smoothing.py diff --git a/phoebe/backend/backends.py b/phoebe/backend/backends.py index 5ba840245..944c61231 100644 --- a/phoebe/backend/backends.py +++ b/phoebe/backend/backends.py @@ -14,7 +14,7 @@ from phoebe.parameters import StringParameter, DictParameter, ArrayParameter, ParameterSet from phoebe.parameters.parameters import _extract_index_from_string from phoebe import dynamics -from phoebe.backend import universe, horizon_analytic +from phoebe.backend import universe, horizon_analytic, contacts_smoothing from phoebe.atmospheres import passbands from phoebe.distortions import roche from phoebe.frontend import io @@ -974,6 +974,23 @@ def _compute_intrinsic_system_at_t0(self, b, compute, # kinds = [b.get_dataset(dataset=ds).kind for ds in datasets] system.update_positions(t0, x0, y0, z0, vx0, vy0, vz0, etheta0, elongan0, eincl0, ignore_effects=True) + + # TEMPERATURE SMOOTHING FOR CONTACTS + smoothing_enabled = b.get_value(qualifier='smoothing_enabled', context='component', **_skip_filter_checks) + if smoothing_enabled: + smoothing_factor = b.get_value(qualifier='smoothing_factor', context='component', **_skip_filter_checks) + primary_mesh, secondary_mesh = system.bodies[0].meshes.values() + coords1 = primary_mesh.roche_coords_for_computations + teffs1 = primary_mesh.teffs + coords2 = secondary_mesh.roche_coords_for_computations + teffs2 = secondary_mesh.teffs + + new_teffs1, new_teffs2 = contacts_smoothing.smooth_teffs(np.array(coords1), np.array(teffs1), + np.array(coords2), np.array(teffs2), + w=smoothing_factor, cutoff=0.1) + primary_mesh.update_columns(teffs=new_teffs1) + secondary_mesh.update_columns(teffs=new_teffs2) + system.populate_observables(t0, ['lc' for dataset in datasets], datasets, ignore_effects=True) @@ -1070,6 +1087,22 @@ def _run_single_time(self, b, i, time, infolist, **kwargs): logger.debug("rank:{}/{} PhoebeBackend._run_single_time: calling system.update_positions at time={}".format(mpi.myrank, mpi.nprocs, time)) system.update_positions(time, xi, yi, zi, vxi, vyi, vzi, ethetai, elongani, eincli, ds=di, Fs=Fi) + # TEMPERATURE SMOOTHING FOR CONTACTS + if i == 0: + smoothing_enabled = b.get_value(qualifier='smoothing_enabled', context='component', **_skip_filter_checks) + if smoothing_enabled: + smoothing_factor = b.get_value(qualifier='smoothing_factor', context='component', **_skip_filter_checks) + primary_mesh, secondary_mesh = system.bodies[0].meshes.values() + coords1 = primary_mesh.roche_coords_for_computations + teffs1 = primary_mesh.teffs + coords2 = secondary_mesh.roche_coords_for_computations + teffs2 = secondary_mesh.teffs + + new_teffs1, new_teffs2 = contacts_smoothing.smooth_teffs(np.array(coords1), np.array(teffs1), + np.array(coords2), np.array(teffs2), + w=smoothing_factor, cutoff=0.) + primary_mesh.update_columns(teffs=new_teffs1) + secondary_mesh.update_columns(teffs=new_teffs2) # Now we need to determine which triangles are visible and handle subdivision # NOTE: this should come after populate_observables so that each subdivided triangle # will have identical local quantities. The only downside to this is that we can't diff --git a/phoebe/backend/contacts_smoothing.py b/phoebe/backend/contacts_smoothing.py new file mode 100644 index 000000000..d12d30346 --- /dev/null +++ b/phoebe/backend/contacts_smoothing.py @@ -0,0 +1,174 @@ +import numpy as np +import matplotlib.pyplot as plt +from itertools import combinations + + +def _isolate_neck(coords_all, teffs_all, cutoff = 0.,component=1, plot=False): + if component == 1: + cond = coords_all[:,0] >= 0+cutoff + elif component == 2: + cond = coords_all[:,0] <= 1-cutoff + else: + raise ValueError + + if plot: + fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(10,5)) + axes[0].scatter(coords_all[cond][:,0], coords_all[cond][:,1]) + axes[1].scatter(coords_all[cond][:,0], teffs_all[cond]) + axes[0].set_xlabel('x') + axes[0].set_ylabel('y') + axes[1].set_xlabel('x') + axes[1].set_ylabel('teff') + fig.tight_layout() + plt.show() + + return coords_all[cond], teffs_all[cond], np.argwhere(cond).flatten() + +def _dist(p1, p2): + (x1, y1, z1), (x2, y2, z2) = p1, p2 + return ((x2 - x1)**2 + (y2 - y1)**2 + (z2 - z1)**2)**0.5 + + +def _isolate_sigma_fitting_regions(coords_neck, teffs_neck, direction='x', cutoff=0., component=1, plot=False): + + distances = [_dist(p1, p2) for p1, p2 in combinations(coords_neck, 2)] + min_dist = np.min(distances[distances != 0]) + + if direction == 'x': + cond = (coords_neck[:,1] >= -0.2*min_dist) & (coords_neck[:,1] <= 0.2*min_dist) + + elif direction == 'y': + + if component == 1: + cond = coords_neck[:,0] <= 0+cutoff+0.15*min_dist + + elif component == 2: + cond = coords_neck[:,0] >= 1-cutoff-0.15*min_dist + + else: + raise ValueError + + else: + raise ValueError + + if plot: + if direction == 'x': + plt.scatter(coords_neck[cond][:,0], teffs_neck[cond]) + plt.show() + elif direction == 'y': + plt.scatter(coords_neck[cond][:,1], teffs_neck[cond]) + plt.show() + + + return coords_neck[cond], teffs_neck[cond] + +def _compute_new_teff_at_neck(coords1, teffs1, coords2, teffs2, w=0.5, offset=0.): + + distances1 = [_dist(p1, p2) for p1, p2 in combinations(coords1, 2)] + min_dist1 = np.min(distances1[distances1 != 0]) + distances2 = [_dist(p1, p2) for p1, p2 in combinations(coords2, 2)] + min_dist2 = np.min(distances2[distances2 != 0]) + + x_neck = np.average((coords1[:,0].max(), coords2[:,0].min())) + amplitude_1 = teffs1.max() + amplitude_2 = teffs2.max() + + teffs_neck1 = teffs1[coords1[:,0] >= coords1[:,0].max() - 0.25*min_dist1] + teffs_neck2 = teffs2[coords2[:,0] <= coords2[:,0].min() + 0.25*min_dist2] + + teff1 = np.average(teffs_neck1) + teff2 = np.average(teffs_neck2) + tavg = w*teff1 + (1-w)*teff2 + if tavg > teffs2.max(): + print('Warning: Tavg > Teff2, setting new temperature to 1 percent of Teff2 max. %i > %i' % (int(tavg), int(teffs2.max()))) + tavg = teffs2.max() - 0.01*teffs2.max() + + return x_neck, tavg + + +def _compute_sigmax(Tavg, x, x0, offset, amplitude): + return ((-1)*(x-x0)**2/np.log((Tavg-offset)/amplitude))**0.5 + + +def _fit_sigma_amplitude(coords, teffs, offset=0., cutoff=0, direction='y', component=1, plot=False): + + def gaussian_1d(x, sigma): + a = 1./sigma**2 + g = offset + amplitude * np.exp(- (a * ((x - x0) ** 2))) + return g + + from scipy.optimize import curve_fit + + if direction == 'y': + coord_ind = 1 + x0 = 0. + elif direction == 'x': + coord_ind = 0 + if component == 1: + x0 = 0+cutoff + elif component == 2: + x0 = 1-cutoff + else: + raise ValueError + else: + raise ValueError + + amplitude = teffs.max() - offset + sigma_0 = 0.5 + result = curve_fit(gaussian_1d, + xdata=coords[:,coord_ind], + ydata=teffs, + p0=(0.5,), + bounds=[0.01,1000]) + + sigma = result[0] + model = gaussian_1d(coords[:,coord_ind], sigma) + + if plot: + plt.scatter(coords[:,coord_ind], teffs) + plt.scatter(coords[:,coord_ind], model) + plt.show() + + return sigma, amplitude, model + + +def _compute_twoD_Gaussian(coords, sigma_x, sigma_y, amplitude, cutoff=0, offset=0., component=1): + y0 = 0. + x0 = 0+cutoff if component == 1 else 1-cutoff + a = 1. / sigma_x ** 2 + b = 1. / sigma_y ** 2 + return offset + amplitude * np.exp(- (a * ((coords[:, 0] - x0) ** 2))) * np.exp(- (b * ((coords[:, 1] - y0) ** 2))) + + + +def smooth_teffs(xyz1, teffs1, xyz2, teffs2, w=0.5, cutoff=0., offset=0.): + + coords_neck_1, teffs_neck_1, cond1 = _isolate_neck(xyz1, teffs1, cutoff = cutoff, component=1, plot=False) + coords_neck_2, teffs_neck_2, cond2 = _isolate_neck(xyz2, teffs2, cutoff = cutoff, component=2, plot=False) + + x_neck, Tavg = _compute_new_teff_at_neck(coords_neck_1, teffs_neck_1, coords_neck_2, teffs_neck_2, w=w, offset=offset) + + sigma_x1 = _compute_sigmax(Tavg, x_neck, x0=0+cutoff, offset=offset, amplitude=teffs_neck_1.max()) + sigma_x2 = _compute_sigmax(Tavg, x_neck, x0=1-cutoff, offset=offset, amplitude=teffs_neck_2.max()) + + coords_fit_y1, teffs_fit_y1 = _isolate_sigma_fitting_regions(coords_neck_1, teffs_neck_1, direction='y', cutoff=cutoff, + component=1, plot=False) + coords_fit_y2, teffs_fit_y2 = _isolate_sigma_fitting_regions(coords_neck_2, teffs_neck_2, direction='y', cutoff=cutoff, + component=2, plot=False) + + sigma_y1, amplitude_y1, model_y1 = _fit_sigma_amplitude(coords_fit_y1, teffs_fit_y1, offset, direction='y', + component=1, plot=False) + sigma_y2, amplitude_y2, model_y2 = _fit_sigma_amplitude(coords_fit_y2, teffs_fit_y2, offset, direction='y', + component=2, plot=False) + + new_teffs1 = _compute_twoD_Gaussian(coords_neck_1, sigma_x1, sigma_y1, teffs_neck_1.max(), + cutoff=cutoff, offset=offset, component=1) + new_teffs2 = _compute_twoD_Gaussian(coords_neck_2, sigma_x2, sigma_y2, teffs_neck_2.max(), + cutoff=cutoff, offset=offset, component=2) + + # print(cond1, len(cond1), len(new_teffs1)) + # print(cond2, len(cond2), len(new_teffs2)) + teffs1[cond1] = new_teffs1 + teffs2[cond2] = new_teffs2 + + return teffs1, teffs2 \ No newline at end of file diff --git a/phoebe/backend/universe.py b/phoebe/backend/universe.py index 6424eac15..776f9c54a 100644 --- a/phoebe/backend/universe.py +++ b/phoebe/backend/universe.py @@ -3326,3 +3326,260 @@ def process_teffs(self, teffs, coords, s=np.array([0., 0., 1.]), t=None): return teffs raise NotImplementedError("teffext=True not yet supported for pulsations") + + +# def dist(p1, p2): +# (x1, y1, z1), (x2, y2, z2) = p1, p2 +# return ((x2 - x1)**2 + (y2 - y1)**2 + (z2 - z1)**2)**0.5 + +# class ContactSmoothing(Feature): +# def __init__(self, w=0.5, cutoff=0, **kwargs): +# self._w = w +# self._cutoff = cutoff + +# @classmethod +# def from_bundle(cls, b, feature): +# """ +# Initialize a Spot feature from the bundle. +# """ + +# # feature_ps = b.get_feature(feature=feature, **_skip_filter_checks) + +# # colat = feature_ps.get_value(qualifier='colat', unit=u.rad, **_skip_filter_checks) +# # longitude = feature_ps.get_value(qualifier='long', unit=u.rad, **_skip_filter_checks) + +# # if len(b.hierarchy.get_stars())>=2: +# # star_ps = b.get_component(component=feature_ps.component, **_skip_filter_checks) +# # orbit_ps = b.get_component(component=b.hierarchy.get_parent_of(feature_ps.component), **_skip_filter_checks) +# # # TODO: how should this handle dpdt? + +# # # we won't use syncpar directly because that is defined wrt sidereal period and we want to make sure +# # # this translated to roche longitude correctly. In the non-apsidal motion case +# # # syncpar = period_anom_orb / period_star +# # period_anom_orb = orbit_ps.get_value(qualifier='period_anom', unit=u.d, **_skip_filter_checks) +# # period_star = star_ps.get_value(qualifier='period', unit=u.d, **_skip_filter_checks) +# # dlongdt = 2*pi * (period_anom_orb/period_star - 1) / period_anom_orb +# # else: +# # star_ps = b.get_component(component=feature_ps.component, **_skip_filter_checks) +# # dlongdt = star_ps.get_value(qualifier='freq', unit=u.rad/u.d, **_skip_filter_checks) +# # longitude += np.pi/2 + +# # radius = feature_ps.get_value(qualifier='radius', unit=u.rad, **_skip_filter_checks) +# # relteff = feature_ps.get_value(qualifier='relteff', unit=u.dimensionless_unscaled, **_skip_filter_checks) + +# # t0 = b.get_value(qualifier='t0', context='system', unit=u.d, **_skip_filter_checks) + +# # return cls(colat, longitude, dlongdt, radius, relteff, t0) + +# @staticmethod +# def _isolate_neck(coords_all, teffs_all, cutoff = 0.,component=1, plot=False): +# if component == 1: +# cond = coords_all[:,0] >= 0+cutoff +# elif component == 2: +# cond = coords_all[:,0] <= 1-cutoff +# else: +# raise ValueError + +# if plot: +# import matplotlib.pyplot as plt +# fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(10,5)) +# axes[0].scatter(coords_all[cond][:,0], coords_all[cond][:,1]) +# axes[1].scatter(coords_all[cond][:,0], teffs_all[cond]) +# axes[0].set_xlabel('x') +# axes[0].set_ylabel('y') +# axes[1].set_xlabel('x') +# axes[1].set_ylabel('teff') +# fig.tight_layout() +# plt.show() + +# return coords_all[cond], teffs_all[cond] + + +# @staticmethod +# def _isolate_sigma_fitting_regions(coords_neck, teffs_neck, direction='x', cutoff=0, component=1, plot=False): +# from itertools import combinations + +# distances = [dist(p1, p2) for p1, p2 in combinations(coords_neck, 2)] +# min_dist = np.min(distances[distances != 0]) + +# if direction == 'x': +# cond = (coords_neck[:,1] >= -0.2*min_dist) & (coords_neck[:,1] <= 0.2*min_dist) + +# elif direction == 'y': + +# if component == 1: +# cond = coords_neck[:,0] <= 0+cutoff+0.15*min_dist + +# elif component == 2: +# cond = coords_neck[:,0] >= 1-cutoff-0.15*min_dist + +# else: +# raise ValueError + +# else: +# raise ValueError + +# if plot: +# import matplotlib.pyplot as plt +# if direction == 'x': +# plt.scatter(coords_neck[cond][:,0], teffs_neck[cond]) +# plt.show() +# elif direction == 'y': +# plt.scatter(coords_neck[cond][:,1], teffs_neck[cond]) +# plt.show() + + +# return coords_neck[cond], teffs_neck[cond] + +# @staticmethod +# def _compute_new_teff_at_neck(coords1, teffs1, coords2, teffs2, w=0.5, offset=0.): + +# from itertools import combinations + +# distances1 = [dist(p1, p2) for p1, p2 in combinations(coords1, 2)] +# min_dist1 = np.min(distances1[distances1 != 0]) +# distances2 = [dist(p1, p2) for p1, p2 in combinations(coords2, 2)] +# min_dist2 = np.min(distances2[distances2 != 0]) + +# x_neck = np.average((coords1[:,0].max(), coords2[:,0].min())) +# amplitude_1 = teffs1.max() +# amplitude_2 = teffs2.max() + +# teffs_neck1 = teffs1[coords1[:,0] >= coords1[:,0].max() - 0.25*min_dist1] +# teffs_neck2 = teffs2[coords2[:,0] <= coords2[:,0].min() + 0.25*min_dist2] + +# teff1 = np.average(teffs_neck1) +# teff2 = np.average(teffs_neck2) +# tavg = w*teff1 + (1-w)*teff2 +# if tavg > teffs2.max(): +# print('Warning: Tavg > Teff2, setting new temperature to 1 percent of Teff2 max. %i > %i' % (int(tavg), int(teffs2.max()))) +# tavg = teffs2.max() - 0.01*teffs2.max() + +# return x_neck, tavg + + +# @staticmethod +# def _compute_sigmax(Tavg, x, x0, offset, amplitude): +# return ((-1)*(x-x0)**2/np.log((Tavg-offset)/amplitude))**0.5 + + +# @staticmethod +# def _fit_sigma_amplitude(coords, teffs, offset, cutoff=0, direction='y', component=1, plot=False): + +# def gaussian_1d(x, sigma): +# a = 1./sigma**2 +# g = offset + amplitude * np.exp(- (a * ((x - x0) ** 2))) +# return g + +# from scipy.optimize import curve_fit + +# if direction == 'y': +# coord_ind = 1 +# x0 = 0. +# elif direction == 'x': +# coord_ind = 0 +# if component == 1: +# x0 = 0+cutoff +# elif component == 2: +# x0 = 1-cutoff +# else: +# raise ValueError +# else: +# raise ValueError + +# amplitude = teffs.max() - offset +# sigma_0 = 0.5 +# result = curve_fit(gaussian_1d, +# xdata=coords[:,coord_ind], +# ydata=teffs, +# p0=(0.5,), +# bounds=[0.01,1000]) + +# sigma = result[0] +# model = gaussian_1d(coords[:,coord_ind], sigma) + +# if plot: +# import matplotlib.pyplot as plt +# plt.scatter(coords[:,coord_ind], teffs) +# plt.scatter(coords[:,coord_ind], model) +# plt.show() + +# return sigma, amplitude, model + + +# @staticmethod +# def _compute_twoD_Gaussian(coords, sigma_x, sigma_y, amplitude, cutoff=0, offset=0, component=1): +# y0 = 0. +# x0 = 0+cutoff if component == 1 else 1-cutoff +# a = 1. / sigma_x ** 2 +# b = 1. / sigma_y ** 2 +# return offset + amplitude * np.exp(- (a * ((coords[:, 0] - x0) ** 2))) * np.exp(- (b * ((coords[:, 1] - y0) ** 2))) + + +# def smooth_teffs(self, xyz1, teffs1, xyz2, teffs2, offset=0.): + +# coords_neck_1, teffs_neck_1 = self._isolate_neck(xyz1, teffs1, cutoff = self._cutoff, component=1, plot=False) +# coords_neck_2, teffs_neck_2 = self._isolate_neck(xyz2, teffs2, cutoff = self._cutoff, component=2, plot=False) + +# x_neck, Tavg = self._compute_new_teff_at_neck(coords_neck_1, teffs_neck_1, coords_neck_2, teffs_neck_2, w=self._w, offset=offset) + +# sigma_x1 = self._compute_sigmax(Tavg, x_neck, x0=0+self._cutoff, offset=offset, amplitude=teffs_neck_1.max()) +# sigma_x2 = self._compute_sigmax(Tavg, x_neck, x0=1-self._cutoff, offset=offset, amplitude=teffs_neck_2.max()) + +# coords_fit_y1, teffs_fit_y1 = self._isolate_sigma_fitting_regions(coords_neck_1, teffs_neck_1, direction='y', cutoff=self._cutoff, +# component=1, plot=False) +# coords_fit_y2, teffs_fit_y2 = self._isolate_sigma_fitting_regions(coords_neck_2, teffs_neck_2, direction='y', cutoff=self._cutoff, +# component=2, plot=False) + +# sigma_y1, amplitude_y1, model_y1 = self._fit_sigma_amplitude(coords_fit_y1, teffs_fit_y1, offset, direction='y', +# component=1, plot=False) +# sigma_y2, amplitude_y2, model_y2 = self._fit_sigma_amplitude(coords_fit_y2, teffs_fit_y2, offset, direction='y', +# component=2, plot=False) + +# new_teffs1 = self. _compute_twoD_Gaussian(coords_neck_1, sigma_x1, sigma_y1, teffs_neck_1.max(), +# cutoff=self._cutoff, offset=offset, component=1) +# new_teffs2 = self._compute_twoD_Gaussian(coords_neck_2, sigma_x2, sigma_y2, teffs_neck_2.max(), +# cutoff=self._cutoff, offset=offset, component=2) + +# return new_teffs1, new_teffs2 + + +# @property +# def proto_coords(self): +# """ +# """ +# return True + + +# def process_teffs(self, teffs, coords, s=np.array([0., 0., 1.]), t=None): +# """ +# Change the local effective temperatures for any values within the +# "cone" defined by the spot. Any teff within the spot will have its +# current value multiplied by the "relteff" factor + +# :parameter array teffs: array of teffs for computations +# :parameter array coords: array of coords for computations +# :t float: current time +# """ + +# # THE COORDINATES NEED TO BE ROCHE! + +# xyz1 = b['value@xyz_elements@mesh@primary'][:,0] +# teffs1 = b['value@teffs@mesh@primary'] +# xyz2 = b['value@xyz_elements@mesh@secondary'][:,0] +# teffs2 = b['value@teffs@mesh@secondary'] + +# # if t is None: +# # # then assume at t0 +# # t = self._t0 + +# # pointing_vector = self.pointing_vector(s,t) +# # logger.debug("spot.process_teffs at t={} with pointing_vector={} and radius={}".format(t, pointing_vector, self._radius)) + +# # cos_alpha_coords = np.dot(coords, pointing_vector) / np.linalg.norm(coords, axis=1) +# # cos_alpha_spot = np.cos(self._radius) + +# # filter_ = cos_alpha_coords > cos_alpha_spot +# # teffs[filter_] = teffs[filter_] * self._relteff + +# # return teffs \ No newline at end of file diff --git a/phoebe/parameters/component.py b/phoebe/parameters/component.py index 257208974..07896f03f 100644 --- a/phoebe/parameters/component.py +++ b/phoebe/parameters/component.py @@ -302,6 +302,8 @@ def envelope(component, **kwargs): params += [FloatParameter(qualifier='pot', latexfmt=r'\Omega_\mathrm{{ {component} }}', value=kwargs.get('pot', 3.5), default_unit=u.dimensionless_unscaled, limits=(0.0,None), description='Potential of the envelope (from the primary component\'s reference)')] params += [FloatParameter(qualifier='pot_min', latexfmt=r'\Omega_\mathrm{{ min, {component} }}', value=kwargs.get('pot_min', 3.5), default_unit=u.dimensionless_unscaled, limits=(0.0,None), description='Critical (minimum) value of the potential to remain a contact')] params += [FloatParameter(qualifier='pot_max', latexfmt=r'\Omega_\mathrm{{ max, {component} }}', value=kwargs.get('pot_max', 3.5), default_unit=u.dimensionless_unscaled, limits=(0.0,None), description='Critical (maximum) value of the potential to remain a contact')] + params += [BoolParameter(qualifier='smoothing_enabled', value=kwargs.get('smoothing_enabled', True), description='Whether to smooth the envelope temperature distribution across the neck')] + params += [FloatParameter(qualifier='smoothing_factor', latexfmt=r'\mathrm{{w}}_\mathrm{{ {component}} }}', value=kwargs.get('smoothing_factor', 0.5), default_unit=u.dimensionless_unscaled, limits=(0.0,1.0), description='Factor for temperature smoothing of the envelope')] # params += [FloatParameter(qualifier='intens_coeff1', value=kwargs.get('intens_coeff1', 1.0), default_unit=u.dimensionless_unscaled, description='')] # params += [FloatParameter(qualifier='intens_coeff2', value=kwargs.get('intens_coeff2', 1.0), default_unit=u.dimensionless_unscaled, description='')] # params += [FloatParameter(qualifier='intens_coeff3', value=kwargs.get('intens_coeff3', 1.0), default_unit=u.dimensionless_unscaled, description='')] From f0fe6057c0d4f08261236a76e227bff270c55c22 Mon Sep 17 00:00:00 2001 From: Angela Pond Date: Thu, 12 Nov 2020 18:20:45 -0500 Subject: [PATCH 02/10] added .save_as_standard_mesh after smoothing --- phoebe/backend/backends.py | 36 ++++++++++++++++++++---------------- 1 file changed, 20 insertions(+), 16 deletions(-) diff --git a/phoebe/backend/backends.py b/phoebe/backend/backends.py index 944c61231..f03ed0e50 100644 --- a/phoebe/backend/backends.py +++ b/phoebe/backend/backends.py @@ -991,6 +991,9 @@ def _compute_intrinsic_system_at_t0(self, b, compute, primary_mesh.update_columns(teffs=new_teffs1) secondary_mesh.update_columns(teffs=new_teffs2) + system.bodies[0]._halves[0].save_as_standard_mesh(primary_mesh) + system.bodies[0]._halves[1].save_as_standard_mesh(secondary_mesh) + system.populate_observables(t0, ['lc' for dataset in datasets], datasets, ignore_effects=True) @@ -1087,22 +1090,23 @@ def _run_single_time(self, b, i, time, infolist, **kwargs): logger.debug("rank:{}/{} PhoebeBackend._run_single_time: calling system.update_positions at time={}".format(mpi.myrank, mpi.nprocs, time)) system.update_positions(time, xi, yi, zi, vxi, vyi, vzi, ethetai, elongani, eincli, ds=di, Fs=Fi) - # TEMPERATURE SMOOTHING FOR CONTACTS - if i == 0: - smoothing_enabled = b.get_value(qualifier='smoothing_enabled', context='component', **_skip_filter_checks) - if smoothing_enabled: - smoothing_factor = b.get_value(qualifier='smoothing_factor', context='component', **_skip_filter_checks) - primary_mesh, secondary_mesh = system.bodies[0].meshes.values() - coords1 = primary_mesh.roche_coords_for_computations - teffs1 = primary_mesh.teffs - coords2 = secondary_mesh.roche_coords_for_computations - teffs2 = secondary_mesh.teffs - - new_teffs1, new_teffs2 = contacts_smoothing.smooth_teffs(np.array(coords1), np.array(teffs1), - np.array(coords2), np.array(teffs2), - w=smoothing_factor, cutoff=0.) - primary_mesh.update_columns(teffs=new_teffs1) - secondary_mesh.update_columns(teffs=new_teffs2) + #TEMPERATURE SMOOTHING FOR CONTACTS + # if i == 0: + # smoothing_enabled = b.get_value(qualifier='smoothing_enabled', context='component', **_skip_filter_checks) + # if smoothing_enabled: + # smoothing_factor = b.get_value(qualifier='smoothing_factor', context='component', **_skip_filter_checks) + # primary_mesh, secondary_mesh = system.bodies[0].meshes.values() + # coords1 = primary_mesh.roche_coords_for_computations + # teffs1 = primary_mesh.teffs + # coords2 = secondary_mesh.roche_coords_for_computations + # teffs2 = secondary_mesh.teffs + + # new_teffs1, new_teffs2 = contacts_smoothing.smooth_teffs(np.array(coords1), np.array(teffs1), + # np.array(coords2), np.array(teffs2), + # w=smoothing_factor, cutoff=0.1) + # primary_mesh.update_columns(teffs=new_teffs1) + # secondary_mesh.update_columns(teffs=new_teffs2) + # Now we need to determine which triangles are visible and handle subdivision # NOTE: this should come after populate_observables so that each subdivided triangle # will have identical local quantities. The only downside to this is that we can't From 4720723ef471327931b625f3b083b1fed31b85fb Mon Sep 17 00:00:00 2001 From: Kyle Date: Thu, 12 Nov 2020 19:34:08 -0500 Subject: [PATCH 03/10] cache "smoothed_teffs" --- phoebe/backend/backends.py | 75 ++++++++++++++++------------- phoebe/backend/universe.py | 99 ++++++++++++++++++++------------------ 2 files changed, 93 insertions(+), 81 deletions(-) diff --git a/phoebe/backend/backends.py b/phoebe/backend/backends.py index f03ed0e50..a0570aaba 100644 --- a/phoebe/backend/backends.py +++ b/phoebe/backend/backends.py @@ -974,25 +974,27 @@ def _compute_intrinsic_system_at_t0(self, b, compute, # kinds = [b.get_dataset(dataset=ds).kind for ds in datasets] system.update_positions(t0, x0, y0, z0, vx0, vy0, vz0, etheta0, elongan0, eincl0, ignore_effects=True) - + # TEMPERATURE SMOOTHING FOR CONTACTS - smoothing_enabled = b.get_value(qualifier='smoothing_enabled', context='component', **_skip_filter_checks) - if smoothing_enabled: - smoothing_factor = b.get_value(qualifier='smoothing_factor', context='component', **_skip_filter_checks) - primary_mesh, secondary_mesh = system.bodies[0].meshes.values() - coords1 = primary_mesh.roche_coords_for_computations - teffs1 = primary_mesh.teffs - coords2 = secondary_mesh.roche_coords_for_computations - teffs2 = secondary_mesh.teffs - - new_teffs1, new_teffs2 = contacts_smoothing.smooth_teffs(np.array(coords1), np.array(teffs1), - np.array(coords2), np.array(teffs2), - w=smoothing_factor, cutoff=0.1) - primary_mesh.update_columns(teffs=new_teffs1) - secondary_mesh.update_columns(teffs=new_teffs2) - - system.bodies[0]._halves[0].save_as_standard_mesh(primary_mesh) - system.bodies[0]._halves[1].save_as_standard_mesh(secondary_mesh) + # smoothing_enabled = b.get_value(qualifier='smoothing_enabled', context='component', **_skip_filter_checks) + # if smoothing_enabled: + # smoothing_factor = b.get_value(qualifier='smoothing_factor', context='component', **_skip_filter_checks) + # primary_mesh, secondary_mesh = system.bodies[0].meshes.values() + # coords1 = primary_mesh.roche_coords_for_computations + # teffs1 = primary_mesh.teffs + # coords2 = secondary_mesh.roche_coords_for_computations + # teffs2 = secondary_mesh.teffs + # + # new_teffs1, new_teffs2 = contacts_smoothing.smooth_teffs(np.array(coords1), np.array(teffs1), + # np.array(coords2), np.array(teffs2), + # w=smoothing_factor, cutoff=0.1) + # # primary_mesh.update_columns(teffs=new_teffs1) + # # secondary_mesh.update_columns(teffs=new_teffs2) + # + # # system.bodies[0]._halves[0].save_as_standard_mesh(primary_mesh) + # # system.bodies[0]._halves[1].save_as_standard_mesh(secondary_mesh) + # system.bodies[0]._halves[0].smoothed_teffs = new_teffs1 + # system.bodies[0]._halves[1].smoothed_teffs = new_teffs2 system.populate_observables(t0, ['lc' for dataset in datasets], datasets, ignore_effects=True) @@ -1091,22 +1093,27 @@ def _run_single_time(self, b, i, time, infolist, **kwargs): system.update_positions(time, xi, yi, zi, vxi, vyi, vzi, ethetai, elongani, eincli, ds=di, Fs=Fi) #TEMPERATURE SMOOTHING FOR CONTACTS - # if i == 0: - # smoothing_enabled = b.get_value(qualifier='smoothing_enabled', context='component', **_skip_filter_checks) - # if smoothing_enabled: - # smoothing_factor = b.get_value(qualifier='smoothing_factor', context='component', **_skip_filter_checks) - # primary_mesh, secondary_mesh = system.bodies[0].meshes.values() - # coords1 = primary_mesh.roche_coords_for_computations - # teffs1 = primary_mesh.teffs - # coords2 = secondary_mesh.roche_coords_for_computations - # teffs2 = secondary_mesh.teffs - - # new_teffs1, new_teffs2 = contacts_smoothing.smooth_teffs(np.array(coords1), np.array(teffs1), - # np.array(coords2), np.array(teffs2), - # w=smoothing_factor, cutoff=0.1) - # primary_mesh.update_columns(teffs=new_teffs1) - # secondary_mesh.update_columns(teffs=new_teffs2) - + smoothing_enabled = b.get_value(qualifier='smoothing_enabled', context='component', **_skip_filter_checks) + if smoothing_enabled and i==0: + smoothing_factor = b.get_value(qualifier='smoothing_factor', context='component', **_skip_filter_checks) + primary_mesh, secondary_mesh = system.bodies[0].meshes.values() + coords1 = primary_mesh.roche_coords_for_computations + teffs1 = primary_mesh.teffs + coords2 = secondary_mesh.roche_coords_for_computations + teffs2 = secondary_mesh.teffs + + new_teffs1, new_teffs2 = contacts_smoothing.smooth_teffs(np.array(coords1), np.array(teffs1), + np.array(coords2), np.array(teffs2), + w=smoothing_factor, cutoff=0.1) + # primary_mesh.update_columns(teffs=new_teffs1) + # secondary_mesh.update_columns(teffs=new_teffs2) + + # print("new_teffs1.median", np.median(new_teffs1)) + # print("new_teffs2.median", np.median(new_teffs2)) + system.bodies[0]._halves[0].smoothed_teffs = new_teffs1 + system.bodies[0]._halves[1].smoothed_teffs = new_teffs2 + + # Now we need to determine which triangles are visible and handle subdivision # NOTE: this should come after populate_observables so that each subdivided triangle # will have identical local quantities. The only downside to this is that we can't diff --git a/phoebe/backend/universe.py b/phoebe/backend/universe.py index 776f9c54a..c9cf9a80f 100644 --- a/phoebe/backend/universe.py +++ b/phoebe/backend/universe.py @@ -1581,22 +1581,27 @@ def _fill_teffs(self, mesh=None, ignore_effects=False, **kwargs): if mesh is None: mesh = self.mesh - # Now we can compute the local temperatures. - # see PHOEBE Legacy scientific reference eq 5.23 - teffs = self.instantaneous_tpole*mesh.gravs.for_computations**0.25 + if hasattr(self, 'smoothed_teffs'): + teffs = self.smoothed_teffs + # print("teffs.median", np.median(teffs)) + else: - if not ignore_effects: - for feature in self.features: - if feature.proto_coords: + # Now we can compute the local temperatures. + # see PHOEBE Legacy scientific reference eq 5.23 + teffs = self.instantaneous_tpole*mesh.gravs.for_computations**0.25 - if self.__class__.__name__ == 'Star_roche_envelope_half' and self.ind_self != self.ind_self_vel: - # then this is the secondary half of a contact envelope - roche_coords_for_computations = np.array([1.0, 0.0, 0.0]) - mesh.roche_coords_for_computations + if not ignore_effects: + for feature in self.features: + if feature.proto_coords: + + if self.__class__.__name__ == 'Star_roche_envelope_half' and self.ind_self != self.ind_self_vel: + # then this is the secondary half of a contact envelope + roche_coords_for_computations = np.array([1.0, 0.0, 0.0]) - mesh.roche_coords_for_computations + else: + roche_coords_for_computations = mesh.roche_coords_for_computations + teffs = feature.process_teffs(teffs, roche_coords_for_computations, s=self.polar_direction_xyz, t=self.time) else: - roche_coords_for_computations = mesh.roche_coords_for_computations - teffs = feature.process_teffs(teffs, roche_coords_for_computations, s=self.polar_direction_xyz, t=self.time) - else: - teffs = feature.process_teffs(teffs, mesh.coords_for_computations, s=self.polar_direction_xyz, t=self.time) + teffs = feature.process_teffs(teffs, mesh.coords_for_computations, s=self.polar_direction_xyz, t=self.time) mesh.update_columns(teffs=teffs) @@ -3379,7 +3384,7 @@ def process_teffs(self, teffs, coords, s=np.array([0., 0., 1.]), t=None): # cond = coords_all[:,0] <= 1-cutoff # else: # raise ValueError - + # if plot: # import matplotlib.pyplot as plt # fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(10,5)) @@ -3391,7 +3396,7 @@ def process_teffs(self, teffs, coords, s=np.array([0., 0., 1.]), t=None): # axes[1].set_ylabel('teff') # fig.tight_layout() # plt.show() - + # return coords_all[cond], teffs_all[cond] @@ -3401,10 +3406,10 @@ def process_teffs(self, teffs, coords, s=np.array([0., 0., 1.]), t=None): # distances = [dist(p1, p2) for p1, p2 in combinations(coords_neck, 2)] # min_dist = np.min(distances[distances != 0]) - + # if direction == 'x': -# cond = (coords_neck[:,1] >= -0.2*min_dist) & (coords_neck[:,1] <= 0.2*min_dist) - +# cond = (coords_neck[:,1] >= -0.2*min_dist) & (coords_neck[:,1] <= 0.2*min_dist) + # elif direction == 'y': # if component == 1: @@ -3412,13 +3417,13 @@ def process_teffs(self, teffs, coords, s=np.array([0., 0., 1.]), t=None): # elif component == 2: # cond = coords_neck[:,0] >= 1-cutoff-0.15*min_dist - + # else: # raise ValueError - + # else: # raise ValueError - + # if plot: # import matplotlib.pyplot as plt # if direction == 'x': @@ -3427,7 +3432,7 @@ def process_teffs(self, teffs, coords, s=np.array([0., 0., 1.]), t=None): # elif direction == 'y': # plt.scatter(coords_neck[cond][:,1], teffs_neck[cond]) # plt.show() - + # return coords_neck[cond], teffs_neck[cond] @@ -3440,39 +3445,39 @@ def process_teffs(self, teffs, coords, s=np.array([0., 0., 1.]), t=None): # min_dist1 = np.min(distances1[distances1 != 0]) # distances2 = [dist(p1, p2) for p1, p2 in combinations(coords2, 2)] # min_dist2 = np.min(distances2[distances2 != 0]) - + # x_neck = np.average((coords1[:,0].max(), coords2[:,0].min())) # amplitude_1 = teffs1.max() # amplitude_2 = teffs2.max() - + # teffs_neck1 = teffs1[coords1[:,0] >= coords1[:,0].max() - 0.25*min_dist1] # teffs_neck2 = teffs2[coords2[:,0] <= coords2[:,0].min() + 0.25*min_dist2] - + # teff1 = np.average(teffs_neck1) # teff2 = np.average(teffs_neck2) # tavg = w*teff1 + (1-w)*teff2 # if tavg > teffs2.max(): # print('Warning: Tavg > Teff2, setting new temperature to 1 percent of Teff2 max. %i > %i' % (int(tavg), int(teffs2.max()))) # tavg = teffs2.max() - 0.01*teffs2.max() - + # return x_neck, tavg - - + + # @staticmethod # def _compute_sigmax(Tavg, x, x0, offset, amplitude): # return ((-1)*(x-x0)**2/np.log((Tavg-offset)/amplitude))**0.5 - + # @staticmethod # def _fit_sigma_amplitude(coords, teffs, offset, cutoff=0, direction='y', component=1, plot=False): - + # def gaussian_1d(x, sigma): # a = 1./sigma**2 # g = offset + amplitude * np.exp(- (a * ((x - x0) ** 2))) # return g # from scipy.optimize import curve_fit - + # if direction == 'y': # coord_ind = 1 # x0 = 0. @@ -3486,26 +3491,26 @@ def process_teffs(self, teffs, coords, s=np.array([0., 0., 1.]), t=None): # raise ValueError # else: # raise ValueError - + # amplitude = teffs.max() - offset # sigma_0 = 0.5 -# result = curve_fit(gaussian_1d, -# xdata=coords[:,coord_ind], -# ydata=teffs, -# p0=(0.5,), +# result = curve_fit(gaussian_1d, +# xdata=coords[:,coord_ind], +# ydata=teffs, +# p0=(0.5,), # bounds=[0.01,1000]) - + # sigma = result[0] # model = gaussian_1d(coords[:,coord_ind], sigma) - + # if plot: # import matplotlib.pyplot as plt # plt.scatter(coords[:,coord_ind], teffs) # plt.scatter(coords[:,coord_ind], model) # plt.show() - + # return sigma, amplitude, model - + # @staticmethod # def _compute_twoD_Gaussian(coords, sigma_x, sigma_y, amplitude, cutoff=0, offset=0, component=1): @@ -3526,24 +3531,24 @@ def process_teffs(self, teffs, coords, s=np.array([0., 0., 1.]), t=None): # sigma_x1 = self._compute_sigmax(Tavg, x_neck, x0=0+self._cutoff, offset=offset, amplitude=teffs_neck_1.max()) # sigma_x2 = self._compute_sigmax(Tavg, x_neck, x0=1-self._cutoff, offset=offset, amplitude=teffs_neck_2.max()) -# coords_fit_y1, teffs_fit_y1 = self._isolate_sigma_fitting_regions(coords_neck_1, teffs_neck_1, direction='y', cutoff=self._cutoff, +# coords_fit_y1, teffs_fit_y1 = self._isolate_sigma_fitting_regions(coords_neck_1, teffs_neck_1, direction='y', cutoff=self._cutoff, # component=1, plot=False) # coords_fit_y2, teffs_fit_y2 = self._isolate_sigma_fitting_regions(coords_neck_2, teffs_neck_2, direction='y', cutoff=self._cutoff, # component=2, plot=False) -# sigma_y1, amplitude_y1, model_y1 = self._fit_sigma_amplitude(coords_fit_y1, teffs_fit_y1, offset, direction='y', +# sigma_y1, amplitude_y1, model_y1 = self._fit_sigma_amplitude(coords_fit_y1, teffs_fit_y1, offset, direction='y', # component=1, plot=False) -# sigma_y2, amplitude_y2, model_y2 = self._fit_sigma_amplitude(coords_fit_y2, teffs_fit_y2, offset, direction='y', +# sigma_y2, amplitude_y2, model_y2 = self._fit_sigma_amplitude(coords_fit_y2, teffs_fit_y2, offset, direction='y', # component=2, plot=False) -# new_teffs1 = self. _compute_twoD_Gaussian(coords_neck_1, sigma_x1, sigma_y1, teffs_neck_1.max(), +# new_teffs1 = self. _compute_twoD_Gaussian(coords_neck_1, sigma_x1, sigma_y1, teffs_neck_1.max(), # cutoff=self._cutoff, offset=offset, component=1) -# new_teffs2 = self._compute_twoD_Gaussian(coords_neck_2, sigma_x2, sigma_y2, teffs_neck_2.max(), +# new_teffs2 = self._compute_twoD_Gaussian(coords_neck_2, sigma_x2, sigma_y2, teffs_neck_2.max(), # cutoff=self._cutoff, offset=offset, component=2) # return new_teffs1, new_teffs2 - + # @property # def proto_coords(self): # """ @@ -3582,4 +3587,4 @@ def process_teffs(self, teffs, coords, s=np.array([0., 0., 1.]), t=None): # # filter_ = cos_alpha_coords > cos_alpha_spot # # teffs[filter_] = teffs[filter_] * self._relteff -# # return teffs \ No newline at end of file +# # return teffs From a9c363f6f8287d7b1aa7774a3b7002022cd77a22 Mon Sep 17 00:00:00 2001 From: Angela Pond Date: Fri, 13 Nov 2020 14:40:03 -0500 Subject: [PATCH 04/10] update weighing to use highest teff2 instead of teff1_neck --- phoebe/backend/backends.py | 20 ++++++++++---------- phoebe/backend/contacts_smoothing.py | 2 +- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/phoebe/backend/backends.py b/phoebe/backend/backends.py index a0570aaba..ce9096269 100644 --- a/phoebe/backend/backends.py +++ b/phoebe/backend/backends.py @@ -984,17 +984,17 @@ def _compute_intrinsic_system_at_t0(self, b, compute, # teffs1 = primary_mesh.teffs # coords2 = secondary_mesh.roche_coords_for_computations # teffs2 = secondary_mesh.teffs - # + # new_teffs1, new_teffs2 = contacts_smoothing.smooth_teffs(np.array(coords1), np.array(teffs1), # np.array(coords2), np.array(teffs2), - # w=smoothing_factor, cutoff=0.1) - # # primary_mesh.update_columns(teffs=new_teffs1) - # # secondary_mesh.update_columns(teffs=new_teffs2) - # - # # system.bodies[0]._halves[0].save_as_standard_mesh(primary_mesh) - # # system.bodies[0]._halves[1].save_as_standard_mesh(secondary_mesh) - # system.bodies[0]._halves[0].smoothed_teffs = new_teffs1 - # system.bodies[0]._halves[1].smoothed_teffs = new_teffs2 + # w=smoothing_factor, cutoff=0.0) + # primary_mesh.update_columns(teffs=new_teffs1) + # secondary_mesh.update_columns(teffs=new_teffs2) + + # system.bodies[0]._halves[0].save_as_standard_mesh(primary_mesh) + # system.bodies[0]._halves[1].save_as_standard_mesh(secondary_mesh) + # system.bodies[0]._halves[0].smoothed_teffs = new_teffs1 + # system.bodies[0]._halves[1].smoothed_teffs = new_teffs2 system.populate_observables(t0, ['lc' for dataset in datasets], datasets, ignore_effects=True) @@ -1104,7 +1104,7 @@ def _run_single_time(self, b, i, time, infolist, **kwargs): new_teffs1, new_teffs2 = contacts_smoothing.smooth_teffs(np.array(coords1), np.array(teffs1), np.array(coords2), np.array(teffs2), - w=smoothing_factor, cutoff=0.1) + w=smoothing_factor, cutoff=0.) # primary_mesh.update_columns(teffs=new_teffs1) # secondary_mesh.update_columns(teffs=new_teffs2) diff --git a/phoebe/backend/contacts_smoothing.py b/phoebe/backend/contacts_smoothing.py index d12d30346..da66c3898 100644 --- a/phoebe/backend/contacts_smoothing.py +++ b/phoebe/backend/contacts_smoothing.py @@ -76,7 +76,7 @@ def _compute_new_teff_at_neck(coords1, teffs1, coords2, teffs2, w=0.5, offset=0. teffs_neck1 = teffs1[coords1[:,0] >= coords1[:,0].max() - 0.25*min_dist1] teffs_neck2 = teffs2[coords2[:,0] <= coords2[:,0].min() + 0.25*min_dist2] - teff1 = np.average(teffs_neck1) + teff1 = np.max(teffs2)#np.average(teffs_neck1) teff2 = np.average(teffs_neck2) tavg = w*teff1 + (1-w)*teff2 if tavg > teffs2.max(): From 7e429c75399ebac0a5968beb2c712c1f53a3cd16 Mon Sep 17 00:00:00 2001 From: Kyle Conroy Date: Fri, 4 Oct 2024 11:18:07 -0400 Subject: [PATCH 05/10] added lateral, isotropic and spotty mixing --- phoebe/backend/backends.py | 42 +++++++++++++---------- phoebe/backend/contacts_smoothing.py | 51 ++++++++++++++++++++++++++-- phoebe/parameters/component.py | 5 +-- 3 files changed, 75 insertions(+), 23 deletions(-) diff --git a/phoebe/backend/backends.py b/phoebe/backend/backends.py index ce9096269..07f44dbae 100644 --- a/phoebe/backend/backends.py +++ b/phoebe/backend/backends.py @@ -975,26 +975,29 @@ def _compute_intrinsic_system_at_t0(self, b, compute, system.update_positions(t0, x0, y0, z0, vx0, vy0, vz0, etheta0, elongan0, eincl0, ignore_effects=True) - # TEMPERATURE SMOOTHING FOR CONTACTS - # smoothing_enabled = b.get_value(qualifier='smoothing_enabled', context='component', **_skip_filter_checks) - # if smoothing_enabled: - # smoothing_factor = b.get_value(qualifier='smoothing_factor', context='component', **_skip_filter_checks) - # primary_mesh, secondary_mesh = system.bodies[0].meshes.values() - # coords1 = primary_mesh.roche_coords_for_computations - # teffs1 = primary_mesh.teffs - # coords2 = secondary_mesh.roche_coords_for_computations - # teffs2 = secondary_mesh.teffs + #TEMPERATURE SMOOTHING FOR CONTACTS + smoothing_enabled = b.get_value(qualifier='smoothing_enabled', context='component', **_skip_filter_checks) - # new_teffs1, new_teffs2 = contacts_smoothing.smooth_teffs(np.array(coords1), np.array(teffs1), - # np.array(coords2), np.array(teffs2), - # w=smoothing_factor, cutoff=0.0) + if smoothing_enabled: + smoothing_method = b.get_value(qualifier='smoothing_method', context='component', **_skip_filter_checks) + smoothing_factor = b.get_value(qualifier='smoothing_factor', context='component', **_skip_filter_checks) + primary_mesh, secondary_mesh = system.bodies[0].meshes.values() + coords1 = primary_mesh.roche_coords_for_computations + teffs1 = primary_mesh.teffs + coords2 = secondary_mesh.roche_coords_for_computations + teffs2 = secondary_mesh.teffs + + new_teffs1, new_teffs2 = contacts_smoothing.smooth_teffs(np.array(coords1), np.array(teffs1), + np.array(coords2), np.array(teffs2), + smoothing_method=smoothing_method) + # w=smoothing_factor, cutoff=0.) # primary_mesh.update_columns(teffs=new_teffs1) # secondary_mesh.update_columns(teffs=new_teffs2) - - # system.bodies[0]._halves[0].save_as_standard_mesh(primary_mesh) - # system.bodies[0]._halves[1].save_as_standard_mesh(secondary_mesh) - # system.bodies[0]._halves[0].smoothed_teffs = new_teffs1 - # system.bodies[0]._halves[1].smoothed_teffs = new_teffs2 + + # print("new_teffs1.median", np.median(new_teffs1)) + # print("new_teffs2.median", np.median(new_teffs2)) + system.bodies[0]._halves[0].smoothed_teffs = new_teffs1 + system.bodies[0]._halves[1].smoothed_teffs = new_teffs2 system.populate_observables(t0, ['lc' for dataset in datasets], datasets, ignore_effects=True) @@ -1094,7 +1097,9 @@ def _run_single_time(self, b, i, time, infolist, **kwargs): #TEMPERATURE SMOOTHING FOR CONTACTS smoothing_enabled = b.get_value(qualifier='smoothing_enabled', context='component', **_skip_filter_checks) + if smoothing_enabled and i==0: + smoothing_method = b.get_value(qualifier='smoothing_method', context='component', **_skip_filter_checks) smoothing_factor = b.get_value(qualifier='smoothing_factor', context='component', **_skip_filter_checks) primary_mesh, secondary_mesh = system.bodies[0].meshes.values() coords1 = primary_mesh.roche_coords_for_computations @@ -1104,7 +1109,8 @@ def _run_single_time(self, b, i, time, infolist, **kwargs): new_teffs1, new_teffs2 = contacts_smoothing.smooth_teffs(np.array(coords1), np.array(teffs1), np.array(coords2), np.array(teffs2), - w=smoothing_factor, cutoff=0.) + smoothing_method=smoothing_method) + # w=smoothing_factor, cutoff=0.) # primary_mesh.update_columns(teffs=new_teffs1) # secondary_mesh.update_columns(teffs=new_teffs2) diff --git a/phoebe/backend/contacts_smoothing.py b/phoebe/backend/contacts_smoothing.py index da66c3898..f0b16e45e 100644 --- a/phoebe/backend/contacts_smoothing.py +++ b/phoebe/backend/contacts_smoothing.py @@ -140,9 +140,7 @@ def _compute_twoD_Gaussian(coords, sigma_x, sigma_y, amplitude, cutoff=0, offset return offset + amplitude * np.exp(- (a * ((coords[:, 0] - x0) ** 2))) * np.exp(- (b * ((coords[:, 1] - y0) ** 2))) - -def smooth_teffs(xyz1, teffs1, xyz2, teffs2, w=0.5, cutoff=0., offset=0.): - +def gaussian_smoothing(xyz1, teffs1, xyz2, teffs2, w=0.5, cutoff=0., offset=0.): coords_neck_1, teffs_neck_1, cond1 = _isolate_neck(xyz1, teffs1, cutoff = cutoff, component=1, plot=False) coords_neck_2, teffs_neck_2, cond2 = _isolate_neck(xyz2, teffs2, cutoff = cutoff, component=2, plot=False) @@ -171,4 +169,51 @@ def smooth_teffs(xyz1, teffs1, xyz2, teffs2, w=0.5, cutoff=0., offset=0.): teffs1[cond1] = new_teffs1 teffs2[cond2] = new_teffs2 + return teffs1, teffs2 + + +def lateral_transfer(t2s, teffs2, lat=0.22, power=0.25): + z2s = t2s[:,2] + lat = 0.15 + power = 0.75 + filt = (z2s>-lat) & (z2s Date: Sat, 14 Nov 2020 13:02:57 -0500 Subject: [PATCH 06/10] parameterization of mixing methods with user-provided power and fixed teffratio from bundle --- phoebe/backend/backends.py | 57 +++++++++++++++------------- phoebe/backend/contacts_smoothing.py | 34 ++++++++++------- phoebe/parameters/component.py | 6 +-- 3 files changed, 55 insertions(+), 42 deletions(-) diff --git a/phoebe/backend/backends.py b/phoebe/backend/backends.py index 07f44dbae..64a9c292e 100644 --- a/phoebe/backend/backends.py +++ b/phoebe/backend/backends.py @@ -976,28 +976,28 @@ def _compute_intrinsic_system_at_t0(self, b, compute, system.update_positions(t0, x0, y0, z0, vx0, vy0, vz0, etheta0, elongan0, eincl0, ignore_effects=True) #TEMPERATURE SMOOTHING FOR CONTACTS - smoothing_enabled = b.get_value(qualifier='smoothing_enabled', context='component', **_skip_filter_checks) + # smoothing_enabled = b.get_value(qualifier='smoothing_enabled', context='component', **_skip_filter_checks) - if smoothing_enabled: - smoothing_method = b.get_value(qualifier='smoothing_method', context='component', **_skip_filter_checks) - smoothing_factor = b.get_value(qualifier='smoothing_factor', context='component', **_skip_filter_checks) - primary_mesh, secondary_mesh = system.bodies[0].meshes.values() - coords1 = primary_mesh.roche_coords_for_computations - teffs1 = primary_mesh.teffs - coords2 = secondary_mesh.roche_coords_for_computations - teffs2 = secondary_mesh.teffs - - new_teffs1, new_teffs2 = contacts_smoothing.smooth_teffs(np.array(coords1), np.array(teffs1), - np.array(coords2), np.array(teffs2), - smoothing_method=smoothing_method) - # w=smoothing_factor, cutoff=0.) - # primary_mesh.update_columns(teffs=new_teffs1) - # secondary_mesh.update_columns(teffs=new_teffs2) - - # print("new_teffs1.median", np.median(new_teffs1)) - # print("new_teffs2.median", np.median(new_teffs2)) - system.bodies[0]._halves[0].smoothed_teffs = new_teffs1 - system.bodies[0]._halves[1].smoothed_teffs = new_teffs2 + # if smoothing_enabled: + # smoothing_method = b.get_value(qualifier='smoothing_method', context='component', **_skip_filter_checks) + # smoothing_factor = b.get_value(qualifier='smoothing_factor', context='component', **_skip_filter_checks) + # primary_mesh, secondary_mesh = system.bodies[0].meshes.values() + # coords1 = primary_mesh.roche_coords_for_computations + # teffs1 = primary_mesh.teffs + # coords2 = secondary_mesh.roche_coords_for_computations + # teffs2 = secondary_mesh.teffs + + # new_teffs1, new_teffs2 = contacts_smoothing.smooth_teffs(np.array(coords1), np.array(teffs1), + # np.array(coords2), np.array(teffs2), + # smoothing_method=smoothing_method) + # # w=smoothing_factor, cutoff=0.) + # # primary_mesh.update_columns(teffs=new_teffs1) + # # secondary_mesh.update_columns(teffs=new_teffs2) + + # # print("new_teffs1.median", np.median(new_teffs1)) + # # print("new_teffs2.median", np.median(new_teffs2)) + # system.bodies[0]._halves[0].smoothed_teffs = new_teffs1 + # system.bodies[0]._halves[1].smoothed_teffs = new_teffs2 system.populate_observables(t0, ['lc' for dataset in datasets], datasets, ignore_effects=True) @@ -1096,11 +1096,15 @@ def _run_single_time(self, b, i, time, infolist, **kwargs): system.update_positions(time, xi, yi, zi, vxi, vyi, vzi, ethetai, elongani, eincli, ds=di, Fs=Fi) #TEMPERATURE SMOOTHING FOR CONTACTS - smoothing_enabled = b.get_value(qualifier='smoothing_enabled', context='component', **_skip_filter_checks) + mixing_enabled = b.get_value(qualifier='mixing_enabled', context='component', **_skip_filter_checks) - if smoothing_enabled and i==0: - smoothing_method = b.get_value(qualifier='smoothing_method', context='component', **_skip_filter_checks) - smoothing_factor = b.get_value(qualifier='smoothing_factor', context='component', **_skip_filter_checks) + if mixing_enabled and i==0: + mixing_method = b.get_value(qualifier='mixing_method', context='component', **_skip_filter_checks) + mixing_power = b.get_value(qualifier='mixing_power', context='component', **_skip_filter_checks) + secondary_teff = b.get_value(qualifier='teff', context='component', component='secondary', **_skip_filter_checks) + primary_teff = b.get_value(qualifier='teff', context='component', component='primary', **_skip_filter_checks) + teff_factor = 1.-secondary_teff/primary_teff + primary_mesh, secondary_mesh = system.bodies[0].meshes.values() coords1 = primary_mesh.roche_coords_for_computations teffs1 = primary_mesh.teffs @@ -1109,7 +1113,8 @@ def _run_single_time(self, b, i, time, infolist, **kwargs): new_teffs1, new_teffs2 = contacts_smoothing.smooth_teffs(np.array(coords1), np.array(teffs1), np.array(coords2), np.array(teffs2), - smoothing_method=smoothing_method) + mixing_method=mixing_method, + mixing_power=mixing_power, teff_factor=teff_factor) # w=smoothing_factor, cutoff=0.) # primary_mesh.update_columns(teffs=new_teffs1) # secondary_mesh.update_columns(teffs=new_teffs2) diff --git a/phoebe/backend/contacts_smoothing.py b/phoebe/backend/contacts_smoothing.py index f0b16e45e..c7a7c0f0a 100644 --- a/phoebe/backend/contacts_smoothing.py +++ b/phoebe/backend/contacts_smoothing.py @@ -172,19 +172,27 @@ def gaussian_smoothing(xyz1, teffs1, xyz2, teffs2, w=0.5, cutoff=0., offset=0.): return teffs1, teffs2 -def lateral_transfer(t2s, teffs2, lat=0.22, power=0.25): +def lateral_transfer(t2s, teffs2, mixing_power=0.5, teff_factor=1.): + x2s = t2s[:,0] + y2s = t2s[:,1] z2s = t2s[:,2] - lat = 0.15 - power = 0.75 + + y2s_neck = y2s[x2s<1] + z2s_neck = z2s[x2s<1] + rs_neck = (y2s_neck**2+z2s_neck**2)**0.5 + lat = np.min(rs_neck) + print('mixing latitude = ', lat) + # lat = 0.15 + # power = 0.75 filt = (z2s>-lat) & (z2s Date: Sat, 14 Nov 2020 14:12:35 -0500 Subject: [PATCH 07/10] updated limits for mixing power to None --- phoebe/parameters/component.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/phoebe/parameters/component.py b/phoebe/parameters/component.py index 460ca5dfd..7f893c54b 100644 --- a/phoebe/parameters/component.py +++ b/phoebe/parameters/component.py @@ -299,7 +299,7 @@ def envelope(component, **kwargs): params += [BoolParameter(qualifier='mixing_enabled', latexfmt=r'\mathrm{{wenabled}}_\mathrm{{ {component} }}', value=kwargs.get('mixing_enabled', True), description='Whether to smooth the envelope temperature distribution across the neck')] - params += [FloatParameter(qualifier='mixing_power', latexfmt=r'\mathrm{{w}}_\mathrm{{ {component} }}', value=kwargs.get('mixing_power', 0.5), default_unit=u.dimensionless_unscaled, limits=(0.0,1.0), description='Power of the thermal mixing of the envelope')] + params += [FloatParameter(qualifier='mixing_power', latexfmt=r'\mathrm{{w}}_\mathrm{{ {component} }}', value=kwargs.get('mixing_power', 0.5), default_unit=u.dimensionless_unscaled, limits=(None,None), description='Power of the thermal mixing of the envelope')] params += [ChoiceParameter(qualifier='mixing_method', latexfmt=r'\mathrm{{smethod}}_\mathrm{{ {component} }}', value=kwargs.get('mixing_method', 'lateral'), choices=['lateral', 'isotropic'], description='Method for thermal mixing of the envelope')] params += [FloatParameter(qualifier='fillout_factor', latexfmt=r'\mathrm{{FF}}_\mathrm{{ {component} }}', value=kwargs.get('fillout_factor', 0.5), default_unit=u.dimensionless_unscaled, limits=(0.0,1.0), description='Fillout-factor of the envelope')] params += [FloatParameter(qualifier='pot', latexfmt=r'\Omega_\mathrm{{ {component} }}', value=kwargs.get('pot', 3.5), default_unit=u.dimensionless_unscaled, limits=(0.0,None), description='Potential of the envelope (from the primary component\'s reference)')] From 1530d586b6e23b6ab012ce5e255415e661722a8a Mon Sep 17 00:00:00 2001 From: Andrej Prsa Date: Sun, 15 Nov 2020 18:12:05 -0500 Subject: [PATCH 08/10] Added "spotted" option to mixing_method. --- phoebe/parameters/component.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/phoebe/parameters/component.py b/phoebe/parameters/component.py index 7f893c54b..7e6c8e1ea 100644 --- a/phoebe/parameters/component.py +++ b/phoebe/parameters/component.py @@ -300,7 +300,7 @@ def envelope(component, **kwargs): params += [BoolParameter(qualifier='mixing_enabled', latexfmt=r'\mathrm{{wenabled}}_\mathrm{{ {component} }}', value=kwargs.get('mixing_enabled', True), description='Whether to smooth the envelope temperature distribution across the neck')] params += [FloatParameter(qualifier='mixing_power', latexfmt=r'\mathrm{{w}}_\mathrm{{ {component} }}', value=kwargs.get('mixing_power', 0.5), default_unit=u.dimensionless_unscaled, limits=(None,None), description='Power of the thermal mixing of the envelope')] - params += [ChoiceParameter(qualifier='mixing_method', latexfmt=r'\mathrm{{smethod}}_\mathrm{{ {component} }}', value=kwargs.get('mixing_method', 'lateral'), choices=['lateral', 'isotropic'], description='Method for thermal mixing of the envelope')] + params += [ChoiceParameter(qualifier='mixing_method', latexfmt=r'\mathrm{{smethod}}_\mathrm{{ {component} }}', value=kwargs.get('mixing_method', 'lateral'), choices=['lateral', 'isotropic', 'spotty'], description='Method for thermal mixing of the envelope')] params += [FloatParameter(qualifier='fillout_factor', latexfmt=r'\mathrm{{FF}}_\mathrm{{ {component} }}', value=kwargs.get('fillout_factor', 0.5), default_unit=u.dimensionless_unscaled, limits=(0.0,1.0), description='Fillout-factor of the envelope')] params += [FloatParameter(qualifier='pot', latexfmt=r'\Omega_\mathrm{{ {component} }}', value=kwargs.get('pot', 3.5), default_unit=u.dimensionless_unscaled, limits=(0.0,None), description='Potential of the envelope (from the primary component\'s reference)')] params += [FloatParameter(qualifier='pot_min', latexfmt=r'\Omega_\mathrm{{ min, {component} }}', value=kwargs.get('pot_min', 3.5), default_unit=u.dimensionless_unscaled, limits=(0.0,None), description='Critical (minimum) value of the potential to remain a contact')] From 9b7116d007de61c5ae8ec52a837cfae63ca551a4 Mon Sep 17 00:00:00 2001 From: Andrej Prsa Date: Sun, 15 Nov 2020 20:33:26 -0500 Subject: [PATCH 09/10] Commented out a print statement. --- phoebe/backend/contacts_smoothing.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/phoebe/backend/contacts_smoothing.py b/phoebe/backend/contacts_smoothing.py index c7a7c0f0a..a4a134c74 100644 --- a/phoebe/backend/contacts_smoothing.py +++ b/phoebe/backend/contacts_smoothing.py @@ -181,7 +181,7 @@ def lateral_transfer(t2s, teffs2, mixing_power=0.5, teff_factor=1.): z2s_neck = z2s[x2s<1] rs_neck = (y2s_neck**2+z2s_neck**2)**0.5 lat = np.min(rs_neck) - print('mixing latitude = ', lat) + # print('mixing latitude = ', lat) # lat = 0.15 # power = 0.75 filt = (z2s>-lat) & (z2s Date: Fri, 4 Oct 2024 11:19:30 -0400 Subject: [PATCH 10/10] updates to distl and stuff to allow plotting --- phoebe/dependencies/distl/distl.py | 5 ++++- phoebe/frontend/bundle.py | 2 +- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/phoebe/dependencies/distl/distl.py b/phoebe/dependencies/distl/distl.py index 032e062fa..82b74d6ae 100644 --- a/phoebe/dependencies/distl/distl.py +++ b/phoebe/dependencies/distl/distl.py @@ -3380,6 +3380,7 @@ def plot_sample(self, size=1e5, **kwargs): plot_uncertainties = kwargs.pop('plot_uncertainties', True) + labels = kwargs.pop('labels', [self._xlabel(dim) for dim in range(self.ndimensions)]) if plot_uncertainties: if plot_uncertainties is True: plot_uncertainties = [1, 2, 3] @@ -3393,9 +3394,10 @@ def plot_sample(self, size=1e5, **kwargs): kwargs.setdefault('levels', [1-_np.exp(-s**2 / 2.) for s in (1,2,3)]) fig = corner.corner(self.sample(size=int(size), cache_sample=False), - labels=[self._xlabel(dim) for dim in range(self.ndimensions)], + labels=labels, quantiles=kwargs.pop('quantiles', None), levels=kwargs.pop('levels', None), + show_titles=False, **kwargs) @@ -4365,6 +4367,7 @@ def _range(dist): range=kwargs.pop('range', [_range(dist) for dist in self.dists]), quantiles=kwargs.pop('quantiles', None), levels=kwargs.pop('levels', None), + show_titles=False **kwargs) if plot_uncertainties: diff --git a/phoebe/frontend/bundle.py b/phoebe/frontend/bundle.py index 92a8f928e..373a1c4c2 100644 --- a/phoebe/frontend/bundle.py +++ b/phoebe/frontend/bundle.py @@ -9082,7 +9082,7 @@ def plot_distribution_collection(self, twig=None, """ plot_kwargs = {} for k in list(kwargs.keys()): - if k in ['plot_uncertainties', 'label', 'xlabel']: + if k in ['plot_uncertainties', 'label', 'labels', 'xlabel']: plot_kwargs[k] = kwargs.pop(k) elif k == 'sample_size': plot_kwargs['size'] = kwargs.pop('sample_size')