<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Passband Luminosity
Setup
Let's first make sure we have the latest version of PHOEBE 2.1 installed. (You can comment out this line if you don't use pip for your installation or don't want to update to the latest release).
Step1:  As always, let's do imports and initialize a logger and a new bundle.  See Building a System for more details.
Step2:  And we'll add a single light curve dataset so that we can see how passband luminosities affect the resulting synthetic light curve model.
Step3:  Lastly, just to make things a bit easier, we'll turn off limb-darkening and irradiation (reflection) and use blackbody atmospheres.
Step4:  Relevant Parameters
The 'pblum_ref' parameter exists for each component-dataset pair and it determines how the intensities for that star in that passband should be scaled, i.e. by the pblum provided by that component ('self') or coupled to the pblum provided by another component.
By default the passband luminosities are coupled (see below for explanations of coupled vs decoupled), with the passband luminosity being defined by the primary component in the system.
Step5:  The 'pblum' parameter is only relevant for each component-dataset pair in which pblum_ref=='self'.  This component will then have its intensities scaled such that they match the value provided by pblum.  In general, a pblum of 4pi will result in an out-of-eclipse flux of ~1.
Step6:  NOTE
Step7:  Now note that only a single pblum parameter is visible.
Step8:  Let's see how changing the value of pblum affects the computed light curve.  By default, pblum is set to be 4 pi, giving a total flux for the primary star of ~1.
Since the secondary star in the default binary is identical to the primary star, we'd expect an out-of-eclipse flux of the binary to be ~2.
Step9:  If we now set pblum to be only 2 pi, we should expect the entire light curve to be scaled in half.
Step10:  And if we halve the temperature of the secondary star - the resulting light curve changes to the new sum of fluxes, where the primary star dominates since the secondary star flux is reduced by a factor of 16, so we expect a total out-of-eclipse flux of ~0.5 + ~0.5/16 = ~0.53.
Step11:  Let us undo our changes before we look at decoupled luminosities.
Step12:  Decoupled Luminosities
The luminosities are decoupled when pblums are provided for the individual components.  To accomplish this, all 'pblum_ref' parameters should be set to 'self'.
Step13:  Now we see that both pblums are available and can have different values.
Step14:  If we set these to 4pi, then we'd expect each star to contribute 1.0 in flux units, meaning the baseline of the light curve should be at approximately 2.0
Step15:  Now let's make a significant temperature-ratio by making a very cool secondary star.  Since the luminosities are decoupled - this temperature change won't affect the resulting light curve very much (compare this to the case above with coupled luminosities).  What is happening here is that even though the secondary star is cooler, its luminosity is being rescaled to the same value as the primary star, so the eclipse depth doesn't change (you would see a similar lack-of-effect if you changed the radii).
Step16:  In most cases you will not want decoupled luminosities as they can easily break the self-consistency of your model.
Now we'll just undo our changes before we look at accessing model luminosities.
Step17:  Accessing Model Luminosities
NEW IN PHOEBE 2.1
Step18:  By default this exposes pblums for all component-dataset pairs in the form of a dictionary.  Alternatively, you can pass a label or list of labels to component and/or dataset.
Step19:  Note that this same logic is applied (at t0) to initialize all passband luminosities within the backend, so does not need to be called before run_compute.
In order to access passband luminosities at times other than t0, you can add a mesh dataset and request the pblum column to be exposed.  For stars that have pblum defined (as opposed to coupled to another star in the system), this value should be equivalent to the value of the parameter (at t0, and in simple circular cases will probably be equivalent at all times).
Let's create a mesh dataset at a few times and then access the synthetic luminosities.
Step20:  Since the luminosities are passband-dependent, they are stored with the same dataset as the light curve (or RV), but with the mesh method, and are available at each of the times at which a mesh was stored.
Step21:  Now let's compare the value of the synthetic luminosities to those of the input pblum
Step22:  In this case, since our two stars are identical, the synthetic luminosity of the secondary star should be the same as the primary (and the same as pblum@primary).
Step23:  However, if we change the temperature of the secondary star again, since the pblums are coupled, we'd expect the synthetic luminosity of the primary to remain fixed but the secondary to decrease.
Step24:  Now, we'll just undo our changes before continuing
Step25:  Role of Pblum
Let's now look at the intensities in the mesh to see how they're being scaled under-the-hood.
Step26:  'abs_normal_intensities' are the intensities per triangle in absolute units, i.e. W/m^3.
Step27:  The values of 'normal_intensities', however, are significantly samller (in this case).  These are the intensities in relative units which will eventually be integrated to give us flux for a light curve.
Step28:  'normal_intensities' are scaled from 'abs_normal_intensities' so that the computed luminosity matches the prescribed luminosity (pblum).
Here we compute the luminosity by summing over each triangle's intensity in the normal direction, and multiply it by pi to account for blackbody intensity emitted in all directions in the solid angle, and by the area of that triangle.
<ASSISTANT_TASK:>
Python Code:
!pip install -I "phoebe>=2.1,<2.2"
Explanation: Passband Luminosity
Setup
Let's first make sure we have the latest version of PHOEBE 2.1 installed. (You can comment out this line if you don't use pip for your installation or don't want to update to the latest release).
End of explanation
%matplotlib inline
import phoebe
from phoebe import u # units
import numpy as np
import matplotlib.pyplot as plt
logger = phoebe.logger()
b = phoebe.default_binary()
Explanation: As always, let's do imports and initialize a logger and a new bundle.  See Building a System for more details.
End of explanation
b.add_dataset('lc', times=np.linspace(0,1,101), dataset='lc01')
Explanation: And we'll add a single light curve dataset so that we can see how passband luminosities affect the resulting synthetic light curve model.
End of explanation
b.set_value_all('ld_func', 'logarithmic')
b.set_value_all('ld_coeffs', [0,0])
b.set_value_all('atm', 'blackbody')
b.set_value('irrad_method', 'none')
Explanation: Lastly, just to make things a bit easier, we'll turn off limb-darkening and irradiation (reflection) and use blackbody atmospheres.
End of explanation
print b['pblum_ref']
print b['pblum_ref@primary']
Explanation: Relevant Parameters
The 'pblum_ref' parameter exists for each component-dataset pair and it determines how the intensities for that star in that passband should be scaled, i.e. by the pblum provided by that component ('self') or coupled to the pblum provided by another component.
By default the passband luminosities are coupled (see below for explanations of coupled vs decoupled), with the passband luminosity being defined by the primary component in the system.
End of explanation
print b['pblum']
Explanation: The 'pblum' parameter is only relevant for each component-dataset pair in which pblum_ref=='self'.  This component will then have its intensities scaled such that they match the value provided by pblum.  In general, a pblum of 4pi will result in an out-of-eclipse flux of ~1.
End of explanation
b['pblum_ref@primary'] = 'self'
b['pblum_ref@secondary'] = 'primary'
Explanation: NOTE: other parameters also affect flux-levels, including limb darkening and distance
Coupled Luminosities
Passband luminosities are considered coupled when a single pblum value is provided, while the passband luminosity of the other component(s) is scaled by the same factor.  To accomplish this, ONE pblum_ref in the system must be set as 'self' and ALL OTHER pbscales must refer to that component. This is the default case, set explicitly by:
End of explanation
print b['pblum']
Explanation: Now note that only a single pblum parameter is visible.
End of explanation
b.run_compute()
afig, mplfig = b.plot(show=True)
Explanation: Let's see how changing the value of pblum affects the computed light curve.  By default, pblum is set to be 4 pi, giving a total flux for the primary star of ~1.
Since the secondary star in the default binary is identical to the primary star, we'd expect an out-of-eclipse flux of the binary to be ~2.
End of explanation
b['pblum@primary'] = 2 * np.pi
b.run_compute()
afig, mplfig = b.plot(show=True)
Explanation: If we now set pblum to be only 2 pi, we should expect the entire light curve to be scaled in half.
End of explanation
b['teff@secondary'] = 0.5 * b.get_value('teff@primary')
print b['teff']
b.run_compute()
afig, mplfig = b.plot(show=True)
Explanation: And if we halve the temperature of the secondary star - the resulting light curve changes to the new sum of fluxes, where the primary star dominates since the secondary star flux is reduced by a factor of 16, so we expect a total out-of-eclipse flux of ~0.5 + ~0.5/16 = ~0.53.
End of explanation
b.set_value_all('teff', 6000)
b.set_value_all('pblum', 4*np.pi)
Explanation: Let us undo our changes before we look at decoupled luminosities.
End of explanation
b.set_value_all('pblum_ref', 'self')
Explanation: Decoupled Luminosities
The luminosities are decoupled when pblums are provided for the individual components.  To accomplish this, all 'pblum_ref' parameters should be set to 'self'.
End of explanation
print b['pblum']
Explanation: Now we see that both pblums are available and can have different values.
End of explanation
b.set_value_all('pblum', 4*np.pi)
b.run_compute()
afig, mplfig = b.plot(show=True)
Explanation: If we set these to 4pi, then we'd expect each star to contribute 1.0 in flux units, meaning the baseline of the light curve should be at approximately 2.0
End of explanation
print b['teff']
b['teff@secondary'] = 3000
b.run_compute()
afig, mplfig = b.plot(show=True)
Explanation: Now let's make a significant temperature-ratio by making a very cool secondary star.  Since the luminosities are decoupled - this temperature change won't affect the resulting light curve very much (compare this to the case above with coupled luminosities).  What is happening here is that even though the secondary star is cooler, its luminosity is being rescaled to the same value as the primary star, so the eclipse depth doesn't change (you would see a similar lack-of-effect if you changed the radii).
End of explanation
b.set_value_all('teff', 6000)
b.set_value_all('pblum', 4*np.pi)
b['pblum_ref@primary'] = 'self'
b['pblum_ref@secondary'] = 'primary'
Explanation: In most cases you will not want decoupled luminosities as they can easily break the self-consistency of your model.
Now we'll just undo our changes before we look at accessing model luminosities.
End of explanation
print b.compute_pblums()
Explanation: Accessing Model Luminosities
NEW IN PHOEBE 2.1: Passband luminosities at t0@system per-star (including following all coupling logic) can be computed and exposed on the fly by calling compute_pblums.
End of explanation
print b.compute_pblums(dataset='lc01', component='primary')
Explanation: By default this exposes pblums for all component-dataset pairs in the form of a dictionary.  Alternatively, you can pass a label or list of labels to component and/or dataset.
End of explanation
b.add_dataset('mesh', times=np.linspace(0,1,5), dataset='mesh01', columns=['areas', 'pblum@lc01', 'ldint@lc01', 'ptfarea@lc01', 'abs_normal_intensities@lc01', 'normal_intensities@lc01'])
b.run_compute()
Explanation: Note that this same logic is applied (at t0) to initialize all passband luminosities within the backend, so does not need to be called before run_compute.
In order to access passband luminosities at times other than t0, you can add a mesh dataset and request the pblum column to be exposed.  For stars that have pblum defined (as opposed to coupled to another star in the system), this value should be equivalent to the value of the parameter (at t0, and in simple circular cases will probably be equivalent at all times).
Let's create a mesh dataset at a few times and then access the synthetic luminosities.
End of explanation
print b.filter(qualifier='pblum', context='model').twigs
Explanation: Since the luminosities are passband-dependent, they are stored with the same dataset as the light curve (or RV), but with the mesh method, and are available at each of the times at which a mesh was stored.
End of explanation
t0 = b.get_value('t0@system')
print b.get_value(qualifier='pblum', time=t0, component='primary', kind='mesh', context='model')
print b.get_value('pblum@primary@dataset')
Explanation: Now let's compare the value of the synthetic luminosities to those of the input pblum
End of explanation
print b.get_value(qualifier='pblum', time=t0, component='primary', kind='mesh', context='model')
print b.get_value(qualifier='pblum', time=t0, component='secondary', kind='mesh', context='model')
Explanation: In this case, since our two stars are identical, the synthetic luminosity of the secondary star should be the same as the primary (and the same as pblum@primary).
End of explanation
b['teff@secondary@component'] = 3000
print b.compute_pblums()
b.run_compute()
print b.get_value(qualifier='pblum', time=t0, component='primary', kind='mesh', context='model')
print b.get_value(qualifier='pblum', time=t0, component='secondary', kind='mesh', context='model')
Explanation: However, if we change the temperature of the secondary star again, since the pblums are coupled, we'd expect the synthetic luminosity of the primary to remain fixed but the secondary to decrease.
End of explanation
b.set_value_all('teff@component', 6000)
Explanation: Now, we'll just undo our changes before continuing
End of explanation
areas = b.get_value(qualifier='areas', dataset='mesh01', time=t0, component='primary', unit='m^2')
ldint = b.get_value(qualifier='ldint', component='primary', time=t0)
ptfarea = b.get_value(qualifier='ptfarea', component='primary', time=t0)
abs_normal_intensities = b.get_value(qualifier='abs_normal_intensities', dataset='lc01', time=t0, component='primary')
normal_intensities = b.get_value(qualifier='normal_intensities', dataset='lc01', time=t0, component='primary')
Explanation: Role of Pblum
Let's now look at the intensities in the mesh to see how they're being scaled under-the-hood.
End of explanation
np.median(abs_normal_intensities)
Explanation: 'abs_normal_intensities' are the intensities per triangle in absolute units, i.e. W/m^3.
End of explanation
np.median(normal_intensities)
Explanation: The values of 'normal_intensities', however, are significantly samller (in this case).  These are the intensities in relative units which will eventually be integrated to give us flux for a light curve.
End of explanation
pblum = b.get_value(qualifier='pblum', component='primary', context='dataset')
print np.sum(normal_intensities * ldint * np.pi * areas) * ptfarea, pblum
Explanation: 'normal_intensities' are scaled from 'abs_normal_intensities' so that the computed luminosity matches the prescribed luminosity (pblum).
Here we compute the luminosity by summing over each triangle's intensity in the normal direction, and multiply it by pi to account for blackbody intensity emitted in all directions in the solid angle, and by the area of that triangle.
End of explanation
<END_TASK>