EclipsingBinaryStar.evaluate

src/tdastro/sources/periodic_variable_star.py

@@ -21,8 +21,9 @@ class PeriodicVariableStar(PeriodicSource, ABC):
     """
 
     def __init__(self, period, epoch, **kwargs):
+        distance = kwargs.pop("distance", None)
         super().__init__(period, epoch, **kwargs)
-        self.add_parameter("distance", required=True, **kwargs)
+        self.add_parameter("distance", value=distance, required=True, **kwargs)
 
     def _evaluate_phases(self, phases, wavelengths, **kwargs):
         """Draw effect-free observations for this object, as a function of phase.
@@ -129,7 +130,7 @@ def _dl_dnu_domega_phases(self, phases, wavelengths, **kwargs):
             self.primary_temperature, self.primary_radius, wavelengths
         )
         secondary_lum = black_body_luminosity_density_per_solid(
-            self.primary_temperature, self.primary_radius, wavelengths
+            self.secondary_temperature, self.secondary_radius, wavelengths
         )
 
         # Distance between the centers of the stars on the plane of the sky
@@ -144,16 +145,14 @@ def _dl_dnu_domega_phases(self, phases, wavelengths, **kwargs):
         # For phases around the main minimum (0) primary star is eclipsed
         secondary_star_closer = np.logical_or(phases < 0.25, phases > 0.75)
         primary_star_overlap_ratio = np.where(
-            secondary_star_closer, overlap_area / (np.pi * self.primary_radius**2), 0
+            secondary_star_closer, overlap_area / (np.pi * self.primary_radius**2), 0.0
         )
-        # For phases around the secondary minimum (0.5) secondary star is eclipsed
-        primary_star_closer = np.logical_and(phases >= 0.25, phases <= 0.75)
         secondary_star_overlap_ratio = np.where(
-            primary_star_closer, overlap_area / (np.pi * self.secondary_radius**2), 0
+            ~secondary_star_closer, overlap_area / (np.pi * self.secondary_radius**2), 0.0
         )
 
-        primary_lum_eclipsed = primary_lum * (1 - primary_star_overlap_ratio)
-        secondary_lum_eclipsed = secondary_lum * (1 - secondary_star_overlap_ratio)
+        primary_lum_eclipsed = primary_lum * (1.0 - primary_star_overlap_ratio)
+        secondary_lum_eclipsed = secondary_lum * (1.0 - secondary_star_overlap_ratio)
 
         # Just in case, we clip negative values to zero
         total_lum = np.where(primary_lum_eclipsed >= 0, primary_lum_eclipsed, 0) + np.where(
@@ -183,8 +182,6 @@ def _norm_star_center_distance(phase_fraction, inclination_degree):
         inclination_radians = np.radians(inclination_degree)
         return np.hypot(np.cos(inclination_radians) * np.cos(phase_radians), np.sin(phase_radians))
 
-    # Math description can be found here:
-    # https://math.stackexchange.com/a/3543551/1348501
     @staticmethod
     def _circle_overlap_area(d, r1, r2):
         """Calculate the area of overlap between two circles.

tests/tdastro/sources/test_periodic_variable_star.py

@@ -1,3 +1,5 @@
+import astropy.constants as const
+import astropy.units as u
 import numpy as np
 import pytest
 from numpy.testing import assert_allclose
@@ -112,6 +114,63 @@ def test_norm_star_center_distance_with_spherical_geometry():
     x, y, z = rotation.apply(np.stack([x_prime, y_prime, z_prime], axis=-1)).T
     # x is looking to the observer, y and z are the observer's plane.
     expected_distance = np.hypot(y, z)
-
     actual_distance = EclipsingBinaryStar._norm_star_center_distance(phase, inclination)
     assert_allclose(actual_distance, expected_distance)
+
+
+def test_eclipsing_binary_star():
+    """Test EclipsingBinaryStar light curve satisfies astrophysical expectations"""
+    points_per_period = 1000
+
+    distance_pc = 10
+    primary_mass = 1.0 * u.M_sun
+    secondary_mass = 0.6 * u.M_sun
+    major_semiaxis = 10 * u.R_sun
+
+    # Mass-radius relation for main-sequence stars
+    primary_radius = 1.0 * u.R_sun * (primary_mass / const.M_sun) ** 0.8
+    secondary_radius = 1.0 * u.R_sun * (secondary_mass / const.M_sun) ** 0.8
+    assert primary_radius > secondary_radius, "Primary should be larger than secondary"
+    # Luminosity-mass relation for main-sequence stars
+    primary_lum = 1.0 * u.L_sun * (primary_mass / const.M_sun) ** 4
+    secondary_lum = 1.0 * u.L_sun * (secondary_mass / const.M_sun) ** 4
+    # Effective temperature from Stefan-Boltzmann law
+    primary_temperature = (primary_lum / (4 * np.pi * primary_radius**2 * const.sigma_sb)) ** 0.25
+    secondary_temperature = (secondary_lum / (4 * np.pi * secondary_radius**2 * const.sigma_sb)) ** 0.25
+    assert primary_temperature > secondary_temperature, "Primary should be hotter than secondary"
+
+    # Third law of Kepler, with gravitation constant
+    period = np.sqrt(4 * np.pi**2 * major_semiaxis**3 / (const.G * (primary_mass + secondary_mass)))
+
+    source = EclipsingBinaryStar(
+        distance=distance_pc,
+        period=period.to_value(u.day),
+        epoch=0.0,
+        major_semiaxis=major_semiaxis.cgs.value,
+        inclination=89.0,
+        primary_radius=primary_radius.cgs.value,
+        primary_temperature=primary_temperature.cgs.value,
+        secondary_radius=secondary_radius.cgs.value,
+        secondary_temperature=secondary_temperature.cgs.value,
+    )
+    # Times in days
+    times = np.linspace(0, 2.0 * period.to_value(u.day), 2 * points_per_period + 1)
+    # Wavelengths in cm
+    wavelengths = np.array([4500, 6000]) * 1e-8
+
+    fluxes = source.evaluate(times, wavelengths)
+
+    # Check the fluxes are positive
+    assert np.all(fluxes >= 0)
+    # Check the fluxes are periodic
+    assert_allclose(fluxes[0], fluxes[-1])
+    assert_allclose(fluxes[:points_per_period], fluxes[points_per_period : 2 * points_per_period])
+    # Check the fluxes are symmetric
+    assert_allclose(fluxes[: points_per_period + 1], fluxes[points_per_period::-1])
+
+    # Check the primary minimum is deeper than everything else
+    assert np.all(fluxes[0] <= fluxes)
+
+    # Check if primary minimum is redder than everything else
+    color = -2.5 * np.log10(fluxes[:, 0] / fluxes[:, 1])
+    assert np.all(color[0] >= color)