Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Dataset Features #891

Draft
wants to merge 3 commits into
base: release-2.5
Choose a base branch
from
Draft

Dataset Features #891

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
09e83f5
refactor features and add hook for dataset features
kecnry Jun 11, 2024
5a88d21
build features
kecnry Jun 11, 2024
ed40baf
use class-defined allowed kinds
kecnry Jun 11, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions phoebe/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -384,6 +384,7 @@ def progressbars(self):
from . import dynamics as dynamics
from . import distortions as distortions
from . import algorithms as algorithms
from . import features as features
import libphoebe

# Shortcut to building logger
Expand Down
305 changes: 3 additions & 302 deletions phoebe/backend/universe.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@
from phoebe.distortions import roche, rotstar
from phoebe.backend import eclipse, oc_geometry, mesh, mesh_wd
from phoebe.utils import _bytes
from phoebe.features import component_features
import libphoebe

from phoebe import u
Expand Down Expand Up @@ -1249,7 +1250,7 @@ def from_bundle(cls, b, component, compute=None,
feature_ps = b.get_feature(feature=feature, **_skip_filter_checks)
if feature_ps.component != component:
continue
feature_cls = globals()[feature_ps.kind.title()]
feature_cls = getattr(component_features, feature_ps.kind.title())
features.append(feature_cls.from_bundle(b, feature))

if conf.devel:
Expand Down Expand Up @@ -3027,304 +3028,4 @@ def populate_observable(self, time, kind, dataset, **kwargs):
"""

for half in self._halves:
half.populate_observable(time, kind, dataset, **kwargs)


################################################################################
################################################################################
################################################################################


class Feature(object):
"""
Note that for all features, each of the methods below will be called. So
changing the coordinates WILL affect the original/intrinsic loggs which
will then be used as input for that method call.

In other words, its probably safest if each feature only overrides a
SINGLE one of the methods. Overriding multiple methods should be done
with great care.

Each feature may or may not require recomputing a mesh, depending on the
kind of change it exacts to the mesh. For example, pulsations will require
recomputing a mesh while spots will not. By default, the mesh will be
recomputed (set in this superclass' `__init__()` method) but inherited
classes should overload `self._remeshing_required`.
"""
def __init__(self, *args, **kwargs):
pass

@property
def _remeshing_required(self):
return True

@property
def proto_coords(self):
"""
Override this to True if all methods (except process_coords*... those
ALWAYS expect protomesh coordinates) are expecting coordinates
in the protomesh (star) frame-of-reference rather than the
current in-orbit system frame-of-reference.
"""
return False

def process_coords_for_computations(self, coords_for_computations, s, t):
"""
Method for a feature to process the coordinates. Coordinates are
processed AFTER scaling but BEFORE being placed in orbit.

NOTE: coords_for_computations affect physical properties only and
not geometric properties (areas, eclipse detection, etc). If you
want to override geometric properties, use the hook for
process_coords_for_observations as well.

Features that affect coordinates_for_computations should override
this method
"""
return coords_for_computations

def process_coords_for_observations(self, coords_for_computations, coords_for_observations, s, t):
"""
Method for a feature to process the coordinates. Coordinates are
processed AFTER scaling but BEFORE being placed in orbit.

NOTE: coords_for_observations affect the geometry only (areas of each
element and eclipse detection) but WILL NOT affect any physical
parameters (loggs, teffs, intensities). If you want to override
physical parameters, use the hook for process_coords_for_computations
as well.

Features that affect coordinates_for_observations should override this method.
"""
return coords_for_observations

def process_loggs(self, loggs, coords, s=np.array([0., 0., 1.]), t=None):
"""
Method for a feature to process the loggs.

Features that affect loggs should override this method
"""
return loggs

def process_teffs(self, teffs, coords, s=np.array([0., 0., 1.]), t=None):
"""
Method for a feature to process the teffs.

Features that affect teffs should override this method
"""
return teffs

class Spot(Feature):
def __init__(self, colat, longitude, dlongdt, radius, relteff, t0, **kwargs):
"""
Initialize a Spot feature
"""
super(Spot, self).__init__(**kwargs)
self._colat = colat
self._longitude = longitude
self._radius = radius
self._relteff = relteff
self._dlongdt = dlongdt
self._t0 = t0

@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)

@property
def _remeshing_required(self):
return False

@property
def proto_coords(self):
"""
"""
return True

def pointing_vector(self, s, time):
"""
s is the spin vector in roche coordinates
time is the current time
"""
t = time - self._t0
longitude = self._longitude + self._dlongdt * t

# define the basis vectors in the spin (primed) coordinates in terms of
# the Roche coordinates.
# ez' = s
# ex' = (ex - s(s.ex)) /|i - s(s.ex)|
# ey' = s x ex'
ex = np.array([1., 0., 0.])
ezp = s
exp = (ex - s*np.dot(s,ex))
eyp = np.cross(s, exp)

return np.sin(self._colat)*np.cos(longitude)*exp +\
np.sin(self._colat)*np.sin(longitude)*eyp +\
np.cos(self._colat)*ezp

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
"""
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

class Pulsation(Feature):
def __init__(self, radamp, freq, l=0, m=0, tanamp=0.0, teffext=False, **kwargs):
self._freq = freq
self._radamp = radamp
self._l = l
self._m = m
self._tanamp = tanamp

self._teffext = teffext

@classmethod
def from_bundle(cls, b, feature):
"""
Initialize a Pulsation feature from the bundle.
"""

feature_ps = b.get_feature(feature=feature, **_skip_filter_checks)
freq = feature_ps.get_value(qualifier='freq', unit=u.d**-1, **_skip_filter_checks)
radamp = feature_ps.get_value(qualifier='radamp', unit=u.dimensionless_unscaled, **_skip_filter_checks)
l = feature_ps.get_value(qualifier='l', unit=u.dimensionless_unscaled, **_skip_filter_checks)
m = feature_ps.get_value(qualifier='m', unit=u.dimensionless_unscaled, **_skip_filter_checks)
teffext = feature_ps.get_value(qualifier='teffext', **_skip_filter_checks)

GM = c.G.to('solRad3 / (solMass d2)').value*b.get_value(qualifier='mass', component=feature_ps.component, context='component', unit=u.solMass, **_skip_filter_checks)
R = b.get_value(qualifier='rpole', component=feature_ps.component, section='component', unit=u.solRad, **_skip_filter_checks)

tanamp = GM/R**3/freq**2

return cls(radamp, freq, l, m, tanamp, teffext)

@property
def proto_coords(self):
"""
"""
return True

def dYdtheta(self, m, l, theta, phi):
if abs(m) > l:
return 0

# TODO: just a quick hack
if abs(m+1) > l:
last_term = 0.0
else:
last_term = Y(m+1, l, theta, phi)

return m/np.tan(theta)*Y(m, l, theta, phi) + np.sqrt((l-m)*(l+m+1))*np.exp(-1j*phi)*last_term

def dYdphi(self, m, l, theta, phi):
return 1j*m*Y(m, l, theta, phi)

def process_coords_for_computations(self, coords_for_computations, s, t):
"""
"""
if self._teffext:
return coords_for_computations

x, y, z, r = coords_for_computations[:,0], coords_for_computations[:,1], coords_for_computations[:,2], np.sqrt((coords_for_computations**2).sum(axis=1))
theta = np.arccos(z/r)
phi = np.arctan2(y, x)

xi_r = self._radamp * Y(self._m, self._l, theta, phi) * np.exp(-1j*2*np.pi*self._freq*t)
xi_t = self._tanamp * self.dYdtheta(self._m, self._l, theta, phi) * np.exp(-1j*2*np.pi*self._freq*t)
xi_p = self._tanamp/np.sin(theta) * self.dYdphi(self._m, self._l, theta, phi) * np.exp(-1j*2*np.pi*self._freq*t)

new_coords = np.zeros(coords_for_computations.shape)
new_coords[:,0] = coords_for_computations[:,0] + xi_r * np.sin(theta) * np.cos(phi)
new_coords[:,1] = coords_for_computations[:,1] + xi_r * np.sin(theta) * np.sin(phi)
new_coords[:,2] = coords_for_computations[:,2] + xi_r * np.cos(theta)

return new_coords

def process_coords_for_observations(self, coords_for_computations, coords_for_observations, s, t):
"""
Displacement equations:

xi_r(r, theta, phi) = a(r) Y_lm (theta, phi) exp(-i*2*pi*f*t)
xi_theta(r, theta, phi) = b(r) dY_lm/dtheta (theta, phi) exp(-i*2*pi*f*t)
xi_phi(r, theta, phi) = b(r)/sin(theta) dY_lm/dphi (theta, phi) exp(-i*2*pi*f*t)

where:

b(r) = a(r) GM/(R^3*f^2)
"""
# TODO: we do want to displace the coords_for_observations, but the x,y,z,r below are from the ALSO displaced coords_for_computations
# if not self._teffext:
# return coords_for_observations

x, y, z, r = coords_for_computations[:,0], coords_for_computations[:,1], coords_for_computations[:,2], np.sqrt((coords_for_computations**2).sum(axis=1))
theta = np.arccos(z/r)
phi = np.arctan2(y, x)

xi_r = self._radamp * Y(self._m, self._l, theta, phi) * np.exp(-1j*2*np.pi*self._freq*t)
xi_t = self._tanamp * self.dYdtheta(self._m, self._l, theta, phi) * np.exp(-1j*2*np.pi*self._freq*t)
xi_p = self._tanamp/np.sin(theta) * self.dYdphi(self._m, self._l, theta, phi) * np.exp(-1j*2*np.pi*self._freq*t)

new_coords = np.zeros(coords_for_observations.shape)
new_coords[:,0] = coords_for_observations[:,0] + xi_r * np.sin(theta) * np.cos(phi)
new_coords[:,1] = coords_for_observations[:,1] + xi_r * np.sin(theta) * np.sin(phi)
new_coords[:,2] = coords_for_observations[:,2] + xi_r * np.cos(theta)

return new_coords

def process_teffs(self, teffs, coords, s=np.array([0., 0., 1.]), t=None):
"""
"""
if not self._teffext:
return teffs

raise NotImplementedError("teffext=True not yet supported for pulsations")
half.populate_observable(time, kind, dataset, **kwargs)
2 changes: 2 additions & 0 deletions phoebe/features/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
from .dataset_features import *
from .component_features import *
Loading
Loading