<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Comparing PHOEBE 2 vs PHOEBE Legacy
NOTE
Step1:  As always, let's do imports and initialize a logger and a new bundle.
Step2:  Adding Datasets and Compute Options
Step3:  Let's add compute options for phoebe using both the new (marching) method for creating meshes as well as the WD method which imitates the format of the mesh used within legacy.
Step4:  Now we add compute options for the 'legacy' backend.
Step5:  And set the two RV datasets to use the correct methods (for both compute options)
Step6:  Let's use the external atmospheres available for both phoebe1 and phoebe2
Step7:  Let's make sure both 'phoebe1' and 'phoebe2wd' use the same value for gridsize
Step8:  Let's also disable other special effect such as heating, gravity, and light-time effects.
Step9:  Finally, let's compute all of our models
Step10:  Plotting
Light Curve
Step11:  Now let's plot the residuals between these two models
Step12:  Dynamical RVs
Since the dynamical RVs don't depend on the mesh, there should be no difference between the 'phoebe2marching' and 'phoebe2wd' synthetic models.  Here we'll just choose one to plot.
Step13:  And also plot the residuals of both the primary and secondary RVs (notice the scale on the y-axis)
Step14:  Numerical (flux-weighted) RVs
<ASSISTANT_TASK:>
Python Code:
#!pip install -I "phoebe>=2.3,<2.4"
Explanation: Comparing PHOEBE 2 vs PHOEBE Legacy
NOTE: PHOEBE 1.0 legacy is an alternate backend and is not installed with PHOEBE 2.  In order to run this backend, you'll need to have PHOEBE 1.0 installed and manually build the python bindings in the phoebe-py directory.
Setup
Let's first make sure we have the latest version of PHOEBE 2.3 installed (uncomment this line if running in an online notebook session such as colab).
End of explanation
import phoebe
from phoebe import u
import numpy as np
import matplotlib.pyplot as plt
phoebe.devel_on()  # needed to use WD-style meshing, which isn't fully supported yet
logger = phoebe.logger()
b = phoebe.default_binary()
b['q'] = 0.7
b['requiv@secondary'] = 0.7
Explanation: As always, let's do imports and initialize a logger and a new bundle.
End of explanation
b.add_dataset('lc', times=np.linspace(0,1,101), dataset='lc01')
b.add_dataset('rv', times=np.linspace(0,1,101), dataset='rvdyn')
b.add_dataset('rv', times=np.linspace(0,1,101), dataset='rvnum')
Explanation: Adding Datasets and Compute Options
End of explanation
b.add_compute(compute='phoebe2marching', irrad_method='none', mesh_method='marching')
b.add_compute(compute='phoebe2wd', irrad_method='none', mesh_method='wd', eclipse_method='graham')
Explanation: Let's add compute options for phoebe using both the new (marching) method for creating meshes as well as the WD method which imitates the format of the mesh used within legacy.
End of explanation
b.add_compute('legacy', compute='phoebe1', irrad_method='none')
Explanation: Now we add compute options for the 'legacy' backend.
End of explanation
b.set_value_all('rv_method', dataset='rvdyn', value='dynamical')
b.set_value_all('rv_method', dataset='rvnum', value='flux-weighted')
Explanation: And set the two RV datasets to use the correct methods (for both compute options)
End of explanation
b.set_value_all('atm', 'extern_planckint')
Explanation: Let's use the external atmospheres available for both phoebe1 and phoebe2
End of explanation
b.set_value_all('gridsize', 30)
Explanation: Let's make sure both 'phoebe1' and 'phoebe2wd' use the same value for gridsize
End of explanation
b.set_value_all('ld_mode', 'manual')
b.set_value_all('ld_func', 'logarithmic')
b.set_value_all('ld_coeffs', [0.,0.])
b.set_value_all('rv_grav', False)
b.set_value_all('ltte', False)
Explanation: Let's also disable other special effect such as heating, gravity, and light-time effects.
End of explanation
b.run_compute(compute='phoebe2marching', model='phoebe2marchingmodel')
b.run_compute(compute='phoebe2wd', model='phoebe2wdmodel')
b.run_compute(compute='phoebe1', model='phoebe1model')
Explanation: Finally, let's compute all of our models
End of explanation
colors = {'phoebe2marchingmodel': 'g', 'phoebe2wdmodel': 'b', 'phoebe1model': 'r'}
afig, mplfig = b['lc01'].plot(c=colors, legend=True, show=True)
Explanation: Plotting
Light Curve
End of explanation
artist, = plt.plot(b.get_value('fluxes@lc01@phoebe2marchingmodel') - b.get_value('fluxes@lc01@phoebe1model'), 'g-')
artist, = plt.plot(b.get_value('fluxes@lc01@phoebe2wdmodel') - b.get_value('fluxes@lc01@phoebe1model'), 'b-')
artist = plt.axhline(0.0, linestyle='dashed', color='k')
ylim = plt.ylim(-0.003, 0.003)
Explanation: Now let's plot the residuals between these two models
End of explanation
afig, mplfig = b.filter(dataset='rvdyn', model=['phoebe2wdmodel', 'phoebe1model']).plot(c=colors, legend=True, show=True)
Explanation: Dynamical RVs
Since the dynamical RVs don't depend on the mesh, there should be no difference between the 'phoebe2marching' and 'phoebe2wd' synthetic models.  Here we'll just choose one to plot.
End of explanation
artist, = plt.plot(b.get_value('rvs@rvdyn@primary@phoebe2wdmodel') - b.get_value('rvs@rvdyn@primary@phoebe1model'), color='b', ls=':')
artist, = plt.plot(b.get_value('rvs@rvdyn@secondary@phoebe2wdmodel') - b.get_value('rvs@rvdyn@secondary@phoebe1model'), color='b', ls='-.')
artist = plt.axhline(0.0, linestyle='dashed', color='k')
ylim = plt.ylim(-1.5e-12, 1.5e-12)
Explanation: And also plot the residuals of both the primary and secondary RVs (notice the scale on the y-axis)
End of explanation
afig, mplfig = b.filter(dataset='rvnum').plot(c=colors, show=True)
artist, = plt.plot(b.get_value('rvs@rvnum@primary@phoebe2marchingmodel', ) - b.get_value('rvs@rvnum@primary@phoebe1model'), color='g', ls=':')
artist, = plt.plot(b.get_value('rvs@rvnum@secondary@phoebe2marchingmodel') - b.get_value('rvs@rvnum@secondary@phoebe1model'), color='g', ls='-.')
artist, = plt.plot(b.get_value('rvs@rvnum@primary@phoebe2wdmodel', ) - b.get_value('rvs@rvnum@primary@phoebe1model'), color='b', ls=':')
artist, = plt.plot(b.get_value('rvs@rvnum@secondary@phoebe2wdmodel') - b.get_value('rvs@rvnum@secondary@phoebe1model'), color='b', ls='-.')
artist = plt.axhline(0.0, linestyle='dashed', color='k')
ylim = plt.ylim(-1e-2, 1e-2)
Explanation: Numerical (flux-weighted) RVs
End of explanation
<END_TASK>