pblums in contact systems are susceptible to neck splitting

[kecnry]

For contact systems, we define the passband luminosities of the individual "stars" that make up the binary. By default these pblums are coupled and provided for the primary component (but could also be decoupled by setting pblum_mode to 'decoupled'). However, since these are susceptible to the "neck" splitting, this can result in somewhat unexpected results for the overall luminosity of the system. This can most easily be seen for a system with identical sub-components:

import phoebe
b=phoebe.default_contact_binary()
b.add_dataset('lc')

print(b.filter(qualifier='pblum*'))
# ParameterSet: 3 parameters
#         pblum_mode@lc01@dataset: component-coupled
#         pblum_component@lc01@dataset: primary
#         pblum@primary@lc01@dataset: 12.5663706144 W

print(b.compute_pblums(pblum_ext=False))
# {'pblum@contact_envelope@lc01': <Quantity 25.17564236 W>, 'pblum@secondary@lc01': <Quantity 12.60927174 W>, 'pblum@primary@lc01': <Quantity 12.56637061 W>}

note that the computed pblum for the primary star matches the input value of 4pi, but since the secondary isn't truly identical at the neck (as far as the mesh is concerned), it has a slightly different computed pblum, and therefore the resulting luminosity of the entire envelope is not exactly 2*pblum@primary@dataset as you might expect.

Increasing the number of triangles seems to cause this discrepancy to converge, which makes sense:

b.set_value('ntriangles', 6000)
print(b.compute_pblums(pblum_ext=False))
# {'pblum@contact_envelope@lc01': <Quantity 25.11051153 W>, 'pblum@secondary@lc01': <Quantity 12.54414091 W>, 'pblum@primary@lc01': <Quantity 12.56637061 W>}
b.set_value('ntriangles', 9000)
print(b.compute_pblums(pblum_ext=False))
# {'pblum@contact_envelope@lc01': <Quantity 25.14510989 W>, 'pblum@secondary@lc01': <Quantity 12.57873928 W>, 'pblum@primary@lc01': <Quantity 12.56637061 W>}

I'm not sure if this is something that we can "fix", if we should reconsider the parameterization, or if we should just clarify this in a tutorial somewhere (perhaps in http://phoebe-project.org/docs/development/tutorials/pblum or in a standalone contact binary pblum tutorial). I have already updated the API docstring for compute_pblums to add a note briefly explaining the effect.

[gecheline]

I think the only case in which this matters is what you mentioned above: if you assume the two components are identical, you expect identical pblums. We could in principle have an "envelope pblum" that in the case of identical components is divided by two, but in all other cases we depend on the mesh to estimate the individual pblums, so, in my experience, this behavior is expected. We should definitely explain this in a tutorial and I think either the existing pblums or contacts tutorials will do (no need for a separate contact binary pblum tutorial).