Sinusoidal Intensities feature

phoebe/backend/universe.py

@@ -1804,9 +1804,6 @@ def _populate_lc(self, dataset, ignore_effects=False, **kwargs):
         else:
             raise NotImplementedError
 
-
-
-
         logger.debug("ld_func={}, ld_coeffs={}, atm={}, ldatm={}".format(ld_func, ld_coeffs, atm, ldatm))
 
         pblum = kwargs.get('pblum', 4*np.pi)
@@ -1931,6 +1928,18 @@ def _populate_lc(self, dataset, ignore_effects=False, **kwargs):
         else:
             raise NotImplementedError("lc_method '{}' not recognized".format(lc_method))
 
+        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 = self.mesh.roche_coords_for_computations
+                    abs_normal_intensities, normal_intensities, abs_intensities, intensities = feature.process_intensities(abs_normal_intensities, normal_intensities, abs_intensities, intensities, roche_coords_for_computations, s=self.polar_direction_xyz, t=self.time)
+                else:
+                    abs_normal_intensities, normal_intensities, abs_intensities, intensities = feature.process_intensities(abs_normal_intensities, normal_intensities, abs_intensities, intensities, self.mesh.coords_for_computations, s=self.polar_direction_xyz, t=self.time)
+
         # TODO: do we really need to store all of these if store_mesh==False?
         # Can we optimize by only returning the essentials if we know we don't need them?
         return {'abs_normal_intensities': abs_normal_intensities,
@@ -3114,6 +3123,16 @@ def process_teffs(self, teffs, coords, s=np.array([0., 0., 1.]), t=None):
         """
         return teffs
 
+    def process_intensities(self, abs_normal_intensities, normal_intensities, abs_intensities, intensities,
+                            coords, s=np.array([0., 0., 1.]), t=None):
+        """
+        Method for a feature to process/modify the intensities.
+
+        Features that affect intensities should override this method
+        """
+        return abs_normal_intensities, normal_intensities, abs_intensities, intensities
+
+
 class Spot(Feature):
     def __init__(self, colat, longitude, dlongdt, radius, relteff, t0, **kwargs):
         """
@@ -3328,3 +3347,44 @@ 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")
+
+class Sinusoidal_Intensities(Feature):
+    def __init__(self, amplitude, frequency, phase0, t0, **kwargs):
+        """
+        Initialize a Spot feature
+        """
+        super(Sinusoidal_Intensities, self).__init__(**kwargs)
+        self._amplitude = amplitude
+        self._frequency = frequency
+        self._phase0 = phase0
+        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)
+        amplitude = feature_ps.get_value(qualifier='amplitude', unit=u.dimensionless_unscaled, **_skip_filter_checks)
+        period = feature_ps.get_value(qualifier='period', unit=u.d, **_skip_filter_checks)
+        frequency = 2 * np.pi / period
+        phase0 = feature_ps.get_value(qualifier='phase0', unit=u.cycle, **_skip_filter_checks)
+        t0 = b.get_value(qualifier='t0', context='system', unit=u.d, **_skip_filter_checks)
+
+        return cls(amplitude, frequency, phase0, t0)
+
+    @property
+    def _remeshing_required(self):
+        return False
+
+    def process_intensities(self, abs_normal_intensities, normal_intensities, abs_intensities, intensities,
+                            coords, s=np.array([0., 0., 1.]), t=None):
+        """
+        """
+        if t is None:
+            # then assume at t0
+            t = self._t0
+
+        factor = 1 + self._amplitude * np.sin(self._frequency * (t - self._t0) + 2 * np.pi * self._phase0)
+        return abs_normal_intensities * factor, normal_intensities * factor, abs_intensities * factor, intensities * factor

phoebe/parameters/feature.py

@@ -18,6 +18,7 @@ def _component_allowed_for_feature(feature_kind, component_kind):
     _allowed['gp_sklearn'] = [None]
     _allowed['gp_celerite2'] = [None]
     _allowed['gaussian_process'] = [None]  # deprecated: remove in 2.5
+    _allowed['sinusoidal_intensities'] = ['star', 'envelope']
 
     return component_kind in _allowed[feature_kind]
 
@@ -28,6 +29,7 @@ def _dataset_allowed_for_feature(feature_kind, dataset_kind):
     _allowed['gp_sklearn'] = ['lc', 'rv', 'lp']
     _allowed['gp_celerite2'] = ['lc', 'rv', 'lp']
     _allowed['gaussian_process'] = ['lc', 'rv', 'lp']  # deprecated: remove in 2.5
+    _allowed['sinusoidal_intensities'] = [None]
 
     return dataset_kind in _allowed[feature_kind]
 
@@ -253,3 +255,19 @@ def gaussian_process(feature, **kwargs):
     """
     logger.warning("gaussian_process is deprecated.  Use gp_celerite2 instead.")
     return gp_celerite2(feature, **kwargs)
+
+def sinusoidal_intensities(feature, **kwargs):
+    """
+    Modify the per-element intensities of a mesh with a sinusoidal factor.
+    """
+    params = []
+
+    params += [FloatParameter(qualifier='amplitude', value=kwargs.get('amplitude', 0.01), default_unit=u.dimensionless_unscaled, description='Amplitude of the sinusoidal intensity variation')]
+    params += [FloatParameter(qualifier='period', value=kwargs.get('period', 1.0), default_unit=u.d, description='Period of the sinusoidal intensity variation')]
+    params += [FloatParameter(qualifier='phase0', value=kwargs.get('phase0', 0.0), default_unit=u.cycle, description='Phase of the sinusoidal intensity variation at time t0@system')]
+
+    # TODO: add frequency (and constraints) and t0 (constrained to t0@system by default)
+
+    constraints = []
+
+    return ParameterSet(params), constraints