How to use the phoebe.BinaryRocheStar function in phoebe

wd_input_file = legacy_basename + '-test.lcin'
    ps,lc,rv = wd.lcin_to_ps(wd_input_file,version='wdphoebe')
    lc['jdstrt'] = time[1]
    lc['jdend'] = time[-1]+time[-1]-time[-2]
    lc['jdinc'] = time[-1]-time[-2]
    lc['indep_type'] = 'time (hjd)'
   
    #-- then create a BodyBag from WD: we want to make sure the output here
    #   is the same as before
    comp1,comp2,binary = wd.wd_to_phoebe(ps,lc,rv)
    star1,lcdep1,rvdep1 = comp1
    star2,lcdep2,rvdep2 = comp2
    mesh1 = phoebe.ParameterSet(frame='phoebe',context='mesh:marching',delta=0.2,alg='c')
    mesh2 = phoebe.ParameterSet(frame='phoebe',context='mesh:marching',delta=0.2,alg='c')
    star1 = phoebe.BinaryRocheStar(star1,binary,mesh1,pbdep=[lcdep1,rvdep1])
    star2 = phoebe.BinaryRocheStar(star2,binary,mesh2,pbdep=[lcdep2,rvdep2])
    wd_bbag = phoebe.BodyBag([star1,star2])
    
    #-- so write a file to compare the two (that's up to you...)
    curve,params = wd.lc(ps,request='curve',light_curve=lc,rv_curve=rv)
    serial_legacy = universe.serialize(system,color=False)
    serial_wildev = universe.serialize(wd_bbag,color=False)
    with open(phoebe_out_file,'w') as ff:
        for line1,line2 in zip(serial_legacy.split('\n'),serial_wildev.split('\n')):
            ff.write('PH:'+line1+'\n')
            ff.write('WD:'+line2+'\n')
    #============ COMPUTATIONS ===============================================
    #-- get mpi-stuff and details, but only if we're not profiling the code.
    if 'cProfile' in globals():
        mpi = None
    else:
        mpi = get_mpi_parameters()

ld_func='claret', ld_coeffs='kurucz_p00')


# Create observations
times = np.linspace(0, orbit['period'],100)
lcobs = phoebe.LCDataSet(time=times, ref=lcdep['ref'])
rvobs1 = phoebe.RVDataSet(time=times, ref=rvdep1['ref'])
rvobs2 = phoebe.RVDataSet(time=times, ref=rvdep2['ref'])

# Body setup
#--------------

# We need two BinaryRocheStars and put them in an orbit:
orbit['ecc'] = 0.1
star1 = phoebe.BinaryRocheStar(comp1, mesh=mesh, pbdep=[lcdep,rvdep1], obs=[rvobs1], orbit=orbit)
star2 = phoebe.BinaryRocheStar(comp2, mesh=mesh, pbdep=[lcdep,rvdep2], obs=[rvobs2], orbit=orbit)

system = phoebe.BodyBag([star1, star2], obs=[lcobs])


# Computation of observables
#-------------------------------

# Let us compute a whole orbit, sampled with 100 points in time, and compute
# the light curve and radial velocity curve. We use the binary eclipse algorithm
# to make computations go faster.

phoebe.observatory.compute(system, eclipse_alg='binary')

# Analysis of results:
#-----------------------

mesh2 = phoebe.ParameterSet(context='mesh:marching', delta=0.05)

lcdep1 = phoebe.PS('lcdep', ld_func='linear', ld_coeffs=[0.5], alb=0, ref='mylc')
lcdep2 = phoebe.PS('lcdep', ld_func='linear', ld_coeffs=[0.5], alb=1, ref='mylc')

# Add scattering properties to the light curve. We'll use the Henyey-Greenstein
# scattering phase function which only has one parameter (the asymmetry parameter)
tools.add_scattering(lcdep2, 'henyey')

# Let's compute the light curve in 200 points
time = np.linspace(0.05, 0.95*orbit['period'], 200)
lcobs = phoebe.LCDataSet(time=time, ref='mylc')

# Create the system consisting of two BinaryRocheStars
star1 = phoebe.BinaryRocheStar(comp1, orbit, mesh1, pbdep=[lcdep1])
star2 = phoebe.BinaryRocheStar(comp2, orbit, mesh2, pbdep=[lcdep2])
system = phoebe.BodyBag([star1, star2], obs=[lcobs])

# Cycle over the different asymmetry values and compute the light curve

for asymmetry in [-0.8,0.0,0.8]:
    
    # Set asymmetry factor
    lcdep2['asymmetry'] = asymmetry

    # Clear previous results and compute the LC
    system.reset_and_clear()
    system.compute()

    # Plot the LC
    phoebe.plotting.plot_lcsyn(system, 'o-', label='Asymmyetry = {}'.format(asymmetry))

# We can easily convert pyWD parameters to pyphoebe parameters:
comp1,comp2,binary,globals = wd.wd_to_phoebe(ps,lc,rv)
star1,lcdep1,rvdep1 = comp1
star2,lcdep2,rvdep2 = comp2

# Set the fineness of the mesh manually, there is no conversion possible between
# a WD-type and pyphoebe-type mesh.
mesh1 = phoebe.ParameterSet(frame='phoebe',context='mesh:marching',delta=0.07,alg='c')
mesh2 = phoebe.ParameterSet(frame='phoebe',context='mesh:marching',delta=0.07,alg='c')

# Body setup
# ----------
star1 = phoebe.BinaryRocheStar(star1,binary,mesh1,pbdep=[lcdep1,rvdep1])
star2 = phoebe.BinaryRocheStar(star2,binary,mesh2,pbdep=[lcdep2,rvdep2])
c1 = time.time()

# Computation of observables
# --------------------------
# Compute WD light curve
curve,params = wd.lc(ps,request='curve',light_curve=lc,rv_curve=rv)
image,params = wd.lc(ps,request='image',light_curve=lc,rv_curve=rv)
c2 = time.time()

mpi = None#phoebe.ParameterSet(context='mpi',np=4)

# Compute pyphoebe light curve
times = curve['indeps']

system = phoebe.BodyBag([star1,star2], globals=globals)
phoebe.observe(system,times,subdiv_num=0,eclipse_alg='binary',

logger.info('\n'+str(lc))
    logger.info('\n'+str(rv))
    
    #-- compute the light curve
    curve,params = wd.lc(ps,request='curve',light_curve=lc,rv_curve=rv)
    
    #-- compare with phoebe if needed
    if do_phoebe:
        comp1,comp2,binary = wd.wd_to_phoebe(ps,lc,rv)
        star1,lcdep1,rvdep1 = comp1
        star2,lcdep2,rvdep2 = comp2
        crit_times = tools.critical_times(binary)
        mesh1 = phoebe.ParameterSet(frame='phoebe',context='mesh:marching',delta=0.075,alg='c')
        mesh2 = phoebe.ParameterSet(frame='phoebe',context='mesh:marching',delta=0.075,alg='c')
        star1 = phoebe.BinaryRocheStar(star1,binary,mesh1,pbdep=[lcdep1,rvdep1])
        star2 = phoebe.BinaryRocheStar(star2,binary,mesh2,pbdep=[lcdep2,rvdep2])
        wd_bbag = phoebe.BodyBag([star1,star2])
        mpi = phoebe.ParameterSet('mpi', np=4)
        phoebe.universe.serialize(wd_bbag,filename='check.pars')
        phoebe.observe(wd_bbag, curve['indeps'], lc=True, rv=True, mpi=mpi)
                #extra_func=[phoebe.observatory.ef_binary_image],
                #extra_func_kwargs=[dict(select='teff',cmap=plt.cm.spectral)])
    
    #-- now do something with it
    plt.figure()
    plt.subplot(121)
    plt.plot(curve['indeps'],curve['lc']/curve['lc'].mean(),'ro-',lw=2,label='WD')
    plt.subplot(122)
    plt.plot(curve['indeps'],curve['rv1'],'ro-',lw=2,label='WD(a)')
    plt.plot(curve['indeps'],curve['rv2'],'ro--',lw=2,label='WD(b)')

apply:

    * The first object you give is the primary, the second one the secondary
    
    * You can set an object equal to the ``None`` variable, in which case
      that component will not be generated.
    
    * An object can be a Body, a BodyBag, a BinaryBag or a list of Bodies.
      In the latter case, the list will first be packed in a BodyBag itself.
      This makes it easy to do *hiearchical packing*.

When creating a Body which requires an orbit to initialize, it is allowed in
this case to set it to ``None``. But only in this case!
"""

star1 = phoebe.BinaryRocheStar(comp1,mesh=mesh,pbdep=[primary_lc1,primary_lc2,spdep])
star2 = phoebe.BinaryRocheStar(comp2,mesh=mesh,pbdep=[secondary_lc1,secondary_lc2,spdep])
system = phoebe.BinaryBag([star1,star2],orbit=orbit)

"""

Under the hood, the latter three calls are equivalent with the
three lines before that.

.. warning:: Difference between ``BodyBag`` and ``BinaryBag``
   
   There is subtle but important difference between ``BodyBag`` and
   ``BinaryBag`` when they are called the following way::
   
       systemA = universe.BodyBag([body1,body2],orbit=orbit)
       systemB = universe.BinaryBag([body1,body2],orbit=orbit)

"""
Now we need to tell the code that the light curve belongs to the whole system,
while the spectroscopic data only pertains to the separate components. You
could do that after initializing the Bodies, but since this is not advised
(you are recommended to initialize the Bodies after creating all ParameterSets
-- and I do mean *all* ParameterSets, also the DataSets), we do it upon
initialization here.

Every Body accepts the keyword ``data``. You should pass the data that belongs
to it. The spectra belong to each component, while the light curve belongs to
the combination of both (i.e. the BodyBag). Note that you are not required
to add data to all Bodies.

"""

star1 = phoebe.BinaryRocheStar(comp1,orbit=orbit,mesh=mesh,
                                 pbdep=[primary_lc1,primary_lc2,spdep],obs=[spdata1])
star2 = phoebe.BinaryRocheStar(comp2,orbit,mesh,
                                 pbdep=[secondary_lc1,secondary_lc2,spdep],obs=[spdata2])
system = phoebe.BodyBag([star1,star2],obs=[lcdata])

"""

Step 6: Computing observables (reprise)
---------------------------------------

The function ``compute_observables`` in ``observatory`` is meant as a
convenience function for unfirom, one-time calculations. In the case of
fitting data, this a bit cumbersome, because you light curves and spectra
are most of the time not taken simultaneously. Thus, you should iterate
over that function, for different time arrays for different observables.
Moreover, in that case of the data, you shouldn't need to manually specify

* The first object you give is the primary, the second one the secondary
    
    * You can set an object equal to the ``None`` variable, in which case
      that component will not be generated.
    
    * An object can be a Body, a BodyBag, a BinaryBag or a list of Bodies.
      In the latter case, the list will first be packed in a BodyBag itself.
      This makes it easy to do *hiearchical packing*.

When creating a Body which requires an orbit to initialize, it is allowed in
this case to set it to ``None``. But only in this case!
"""

star1 = phoebe.BinaryRocheStar(comp1,mesh=mesh,pbdep=[primary_lc1,primary_lc2,spdep])
star2 = phoebe.BinaryRocheStar(comp2,mesh=mesh,pbdep=[secondary_lc1,secondary_lc2,spdep])
system = phoebe.BinaryBag([star1,star2],orbit=orbit)

"""

Under the hood, the latter three calls are equivalent with the
three lines before that.

.. warning:: Difference between ``BodyBag`` and ``BinaryBag``
   
   There is subtle but important difference between ``BodyBag`` and
   ``BinaryBag`` when they are called the following way::
   
       systemA = universe.BodyBag([body1,body2],orbit=orbit)
       systemB = universe.BinaryBag([body1,body2],orbit=orbit)
    
   In the first call, ``systemA`` will be put as a whole in the orbit ``orbit``.

b. *With data*, you can use ``observatory.compute``. It requires
        no extra input, since all information is derived from the data. For
        each added dataset, the synthetics will be computed at the relevant
        time points (see Step 5).
    
The first option is covered in earlier and later tutorials, in this tutorial
the automatic computation option is explored.

We've fiddled a bit above with the creation of the Bodies. Just for clarity,
we start over quickly. We make a circular version of V380 Cyg:

"""

comp1,comp2,orbit = phoebe.create.from_library('V380_Cyg')
star1 = phoebe.BinaryRocheStar(comp1,orbit=orbit,mesh=mesh,pbdep=[primary_lc1,primary_lc2,spdep])
star2 = phoebe.BinaryRocheStar(comp2,orbit=orbit,mesh=mesh,pbdep=[secondary_lc1,secondary_lc2,spdep])
system = phoebe.BodyBag([star1,star2])
orbit['ecc'] = 0.


"""

.. admonition:: ParameterSets are passed by reference.

   You can (and should) always work with the ParameterSets you defined before
   creating the Bodies. They are passed by reference, so changing them outside
   of the Body will propogate into the Bodies too.


Let us define a list (or array for that matter) of times at which we want to
compute observables. The list of times does not need to be equidistant. For all I

Possibility 2: Via separate components
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

It is easy to create a BodyBag yourself, but of course you need to
create the two BinaryRocheStars first. To do that, you also need to specify
information on how to compute the mesh. This is the last parameterSet that
is required:

"""
mesh = phoebe.ParameterSet(context='mesh:marching')

# The ``BinaryRocheStars`` are then easily created and (optionally) put in
# a ``BodyBag``, after importing the universe:

star1 = phoebe.BinaryRocheStar(comp1, orbit=orbit, mesh=mesh, pbdep=[primary_lc1, primary_lc2, spdep])
star2 = phoebe.BinaryRocheStar(comp2, orbit=orbit, mesh=mesh, pbdep=[secondary_lc1, secondary_lc2, spdep])
system = phoebe.BodyBag([star1, star2])

"""

.. admonition:: Adding one or more pbdeps

    You can add only one pbdep (e.g. only ``primary_lc1`` and ``secondary_lc1``),
    but you need to add one pbdep for each component. Adding a pbdep basically
    tells the code *how* to compute something. These specifications can be
    different for each body, so you need to give each body this information!


Possibility 3: Via hierarchical packing
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

For the specific case of a binary, there is a class called ``BinaryBag``, which