#!pip install -I "phoebe>=2.4,<2.5"

import phoebe
from phoebe import u # units
import numpy as np
import matplotlib.pyplot as plt

logger = phoebe.logger()

"""
Explanation: Contact Binary Hierarchy
Setup
Let's first make sure we have the latest version of PHOEBE 2.4 installed (uncomment this line if running in an online notebook session such as colab).
End of explanation
"""


b_cb = phoebe.default_binary(contact_binary=True)

"""
Explanation: Here we'll initialize a default binary, but ask for it to be created as a contact system.
End of explanation
"""


b_detached = phoebe.default_binary()

"""
Explanation: We'll compare this to the default detached binary
End of explanation
"""


print(b_detached.hierarchy)

print(b_cb.hierarchy)

"""
Explanation: Hierarchy
Let's first look at the hierarchy of the default detached binary, and then compare that to the hierarchy of the overcontact system
End of explanation
"""


print(b_cb.filter(component='contact_envelope', kind='envelope', context='component'))

print(b_cb.filter(component='primary', kind='star', context='component'))

b_cb['requiv@primary'] = 1.5

b_cb['pot@contact_envelope@component']

b_cb['fillout_factor@contact_envelope@component']

b_cb['requiv@secondary@component']

"""
Explanation: As you can see, the overcontact system has an additional "component" with method "envelope" and component label "contact_envelope".
Next let's look at the parameters in the envelope and star components. You can see that most of parameters in the envelope class are constrained, while the equivalent radius of the primary is unconstrained. The value of primary equivalent radius constrains the potential and fillout factor of the envelope, as well as the equivalent radius of the secondary.
End of explanation
"""


b_detached.add_component('envelope', component='contact_envelope')

hier = phoebe.hierarchy.binaryorbit(b_detached['binary'], b_detached['primary'], b_detached['secondary'], b_detached['contact_envelope'])
print(hier)

b_detached.filter(context='constraint',constraint_func='pitch',component='primary')

b_detached.set_hierarchy(hier)

print(b_detached.hierarchy)

"""
Explanation: Now, of course, if we didn't originally know we wanted a contact binary and built the default detached system, we could still turn it into an contact binary just by changing the hierarchy.
End of explanation
"""


print(b_detached.run_checks())

"""
Explanation: However, since our system was detached, the system is not overflowing, and therefore doesn't pass system checks
End of explanation
"""


b_detached['pot@component']

b_detached['requiv@secondary@component']

"""
Explanation: And because of this, the potential and requiv@secondary constraints cannot be updated from the constraints.
End of explanation
"""


hier = phoebe.hierarchy.binaryorbit(b_detached['binary'], b_detached['primary'], b_detached['secondary'])
print(hier)

b_detached.set_hierarchy(hier)

print(b_detached.hierarchy)

"""
Explanation: Likewise, we can make a contact system detached again simply by removing the envelope from the hierarchy.  The parameters themselves will still exist (unless you remove them), so you can always just change the hierarchy again to change back to an overcontact system.
End of explanation
"""


b_detached['requiv@secondary@component'] = 1.0

"""
Explanation: Although the constraints have been removed, PHOEBE has lost the original value of the secondary radius (because of the failed contact constraints), so we'll have to reset that here as well.
End of explanation
"""


b_cb.add_dataset('mesh', compute_times=[0], dataset='mesh01')

b_cb.add_dataset('orb', compute_times=np.linspace(0,1,201), dataset='orb01')

b_cb.add_dataset('lc', times=np.linspace(0,1,21), dataset='lc01')

b_cb.add_dataset('rv', times=np.linspace(0,1,21), dataset='rv01')

"""
Explanation: Adding Datasets
End of explanation
"""


b_detached.add_dataset('mesh', compute_times=[0], dataset='mesh01')

b_detached.add_dataset('orb', compute_times=np.linspace(0,1,201), dataset='orb01')

b_detached.add_dataset('lc', times=np.linspace(0,1,21), dataset='lc01')

b_detached.add_dataset('rv', times=np.linspace(0,1,21), dataset='rv01')

"""
Explanation: For comparison, we'll do the same to our detached system
End of explanation
"""


b_cb.run_compute(irrad_method='none')

b_detached.run_compute(irrad_method='none')

"""
Explanation: Running Compute
End of explanation
"""


print(b_cb['mesh01@model'].components)

print(b_detached['mesh01@model'].components)

"""
Explanation: Synthetics
To ensure compatibility with computing synthetics in detached and semi-detached systems in Phoebe, the synthetic meshes for our overcontact system are attached to each component separetely, instead of the contact envelope.
End of explanation
"""


afig, mplfig = b_cb['mesh01@model'].plot(x='ws', show=True)

afig, mplfig = b_detached['mesh01@model'].plot(x='ws', show=True)

"""
Explanation: Plotting
Meshes
End of explanation
"""


afig, mplfig = b_cb['orb01@model'].plot(x='ws',show=True)

afig, mplfig = b_detached['orb01@model'].plot(x='ws',show=True)

"""
Explanation: Orbits
End of explanation
"""


afig, mplfig = b_cb['lc01@model'].plot(show=True)

afig, mplfig = b_detached['lc01@model'].plot(show=True)

"""
Explanation: Light Curves
End of explanation
"""


afig, mplfig = b_cb['rv01@model'].plot(show=True)

afig, mplfig = b_detached['rv01@model'].plot(show=True)

"""
Explanation: RVs
End of explanation
"""


