Is it possible to expose individual component fluxes (light curves)?

[jsinkbaek]

Hi everyone!

I have a niche case where I would like to see the light curve split into two. Basically, the contribution of the primary and the secondary, as two separate outputs.
Is that something that is doable with the way PHOEBE is designed? And could you point me to either a ressource, or the code-snippet that handles the "output" light curve?

All the best,
Jeppe

[kecnry]

Hi Jeppe - this has been on the wishlist for a while, but would require a significant rewrite of some code to do in one single step and I don't see that being added anytime in the near future unless someone wants to take a fairly big project on. But this was a fun challenge to get some examples with what currently is supported.

If you do NOT want to see eclipses

You can tell one of the stars to skip meshing (emit no light, do NOT eclipse, but still have mass) with distortion_method='none'. This effectively gives you the contribution from one star with the other star treated as a point mass (so eclipses will disappear). Unfortunately this isn't currently dataset-dependent, so the only way to do it is with multiple calls to run_compute or multiple datasets and compute options.

Multiple calls to run_compute:

import phoebe
b = phoebe.default_binary()
b.add_dataset('lc', compute_phases=phoebe.linspace(0,1,101))
# b.set_value('pblum_mode', 'decoupled') # see note below
b.set_value('distortion_method', component='secondary', value='none')
b.run_compute(model='primary_only')
b.set_value('distortion_method', component='secondary', value='roche')
b.set_value('distortion_method', component='primary', value='none')
b.run_compute(model='secondary_only')
b.plot(legend=True, show=True)

Putting together this example actually exposed an interesting behavior here that I'm not sure if it should be considered a bug or just a consequence of the parameters (I'd be curious to hear your, and others', opinion). Since ld_mode='component-coupled' and the primary pblum is defining the system, but distortion_method@primary='none' is effectively saying there is no light coming from the primary star, the secondary is scaled to zero. If you uncomment the line to decouple the luminosities (you probably don't want to do this in practice as is done in this script), you get the following instead:

This makes sense with what you'd expect. In a real system, you could do this by using passband-coupling with both components on, manually compute the two passbands, and then set them after decoupling:

import phoebe
b = phoebe.default_binary()
b.add_dataset('lc', compute_phases=phoebe.linspace(0,1,101))
# set any system parameters, only call the following with FIXED parameters
pblums = b.compute_pblums()
b.set_value('pblum_mode', 'decoupled')
for twig, quantity in pblums:
    b.set_value(twig, quantity)
...

In theory, I would have also expected the following to work for a single call to run_compute with multiple datasets (ignoring the passband issue discussed above), but this fails during plotting (I created #840 to track this bug - but I doubt many people will run into this so probably isn't a high priority to track down and fix):

import phoebe
b = phoebe.default_binary()
b.add_dataset('lc', compute_phases=phoebe.linspace(0,1,101), dataset='lc_primary_only')
b.add_dataset('lc', compute_phases=phoebe.linspace(0,1,101), dataset='lc_secondary_only')
b.add_compute('phoebe', 
              distortion_method={'primary': 'roche', 'secondary': 'none'},
              enabled={'lc_primary_only': True, 'lc_secondary_only': False},
              compute='comp_primary_only')
b.add_compute('phoebe',
              distortion_method={'primary': 'none', 'secondary': 'roche'},
              enabled={'lc_primary_only': False, 'lc_secondary_only': True},
              compute='comp_secondary_only')
b.run_compute(compute=['comp_primary_only', 'comp_secondary_only'])
print(b.filter(context='model')  # does show the arrays are populated correctly 
b.plot(legend=True, show=True)  # raises an error currently

If you DO want to see eclipses

Now things are really inconvenient, and the only current way to get what you want is by populated and exposing the mesh at each time and integrating over a single star to build your own light curve.

import phoebe
import numpy as np
import matplotlib.pyplot as plt

times = np.linspace(0,1,11)

b = phoebe.default_binary()
b.add_dataset('lc', compute_times=times, dataset='lc01')
b.add_dataset('mesh', include_times='lc01', columns=['areas', 'mus', 'visibilities', 'intensities@lc01', 'ptfarea@lc01'])
b.run_compute()

lcs = {'primary': np.zeros_like(times), 'secondary': np.zeros_like(times)}
for i, time in enumerate(times):
    for component in ['primary', 'secondary']:
        mesh_ps = b.filter(time=time, component=component, context='model')
        intensities = mesh_ps.get_value('intensities@lc01')
        areas = mesh_ps.get_value('areas')
        mus = mesh_ps.get_value('mus')
        visibilities = mesh_ps.get_value('visibilities')
        ptfarea = mesh_ps.get_value('ptfarea@lc01')
    
        flux = np.nansum(intensities * areas * mus * visibilities) * ptfarea
        lcs[component][i] = flux

plt.plot(times, lcs['primary'], 'b.-')
plt.plot(times, lcs['secondary'], 'r.-')
plt.show()

[jsinkbaek]

Hi Kyle,

Thank you very much for the extensive details in the response, this will make everything a lot easier!

Unfortunately, I think I want the eclipses as well. I'm trying to estimate what the amplitude-modulation from the binarity has on the power spectrum with oscillations originating from just one of the stellar components.

I think, wrt. the comment you had on point-mass and component-coupled. To me it sounds like a pretty natural consequence of the parameters. If the luminosity of one is 0, and you define the second component luminosity as a fraction of the first, then it'll be 0.

In terms of wishlist. My two cents is that it would be nice to have this thing output. Either as flux1 and flux2, or as an array of light contribution values of each component going from 0 to 1. Both are equivalent.

Again, thank you a lot for your time!

All the best,
Jeppe