How to use the phoebe.BodyBag function in phoebe

mesh1 = phoebe.PS('mesh:marching', delta=0.1)
mesh2 = phoebe.PS('mesh:marching', delta=0.03, maxpoints=40000)

# Put the system at about 1 AU:

globals = phoebe.PS('position', distance=(1,'au')) 

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

sun = phoebe.BinaryStar(sun, mesh=mesh1, orbit=orbit, pbdep=[lcdep1])
mercury = phoebe.BodyBag(phoebe.Star(mercury, mesh=mesh2, pbdep=[lcdep2]),
                         orbit=orbit)

system = phoebe.BodyBag([sun, mercury], obs=[obs], globals=globals)


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

system.compute(heating=True, refl=True, refl_num=1, eclipse_alg='graham')


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

# Let's make a temperature map, with the extremes being the temperature of the
# planet in absence of radiation, and the temperature of maximum heating.

out = system[1].plot2D(select='teff',vmin=100, vmax=700)
cbar = plt.colorbar(out[-1])

c1label=sun['label'], c2label=mercury['label'],incl=86.62,
                    per0=(29.124,'deg'))

    lcdep1 = phoebe.PS('lcdep', atm='kurucz', ld_coeffs='kurucz', ld_func='claret', ref='apparent')
    lcdep2 = phoebe.PS('lcdep', atm=atm, ld_func='uniform', ref='apparent', alb=alb)
    obs = phoebe.LCDataSet(time=np.array([orbit['period']*0.25]), columns=['time'], ref=lcdep1['ref'])

    mesh1 = phoebe.PS('mesh:marching', delta=0.1)
    mesh2 = phoebe.PS('mesh:marching', delta=0.04, maxpoints=22000)

    globals = phoebe.PS('position', distance=(1,'au')) 
    
    sun = phoebe.BinaryStar(sun, mesh=mesh1, orbit=orbit, pbdep=[lcdep1])
    mercury = phoebe.BinaryStar(mercury, mesh=mesh2, orbit=orbit, pbdep=[lcdep2])

    system = phoebe.BodyBag([sun, mercury], obs=[obs], position=globals)
    system.compute(heating=True, refl=True, refl_num=1, boosting_alg='none')
    return system

mesh = phoebe.ParameterSet(context='mesh:marching', delta=0.05, alg='c')
lcdep1 = create.dep_from_object(sun,context='lcdep', passband='JOHNSON.V',
                                ref='Visual')
lcdep2 = create.dep_from_object(venus,context='lcdep', passband='JOHNSON.V',
                                ref='Visual')

globals = phoebe.ParameterSet('position', distance=(4.84813681108e-06,'pc'))

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

# Then the Sun and Venus are easily created as ``BinaryStars``.
bsun = phoebe.BinaryStar(sun, orbit, mesh, pbdep=lcdep1)
bvenus = phoebe.BinaryStar(venus, orbit, mesh, pbdep=lcdep2)

system = phoebe.BodyBag([bsun, bvenus], globals=globals)

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

#observatory.compute(system,[0.35*orbit['period']],lc=True,reflection=True,heating=False,circular=False)
phoebe.observe(system,[0.41*orbit['period']],lc=True,refl=True,heating=False)

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

# We make some nice images of Venus and the SUn, and convert the computed
# intensity to Johnson magnitudes. For fun, we also check what the distance
# to the two bodies is.
bvenus.plot2D(ref='Visual',select='proj',savefig='venus_proj')
bvenus.plot2D(ref='Visual',select='teff',cmap='eye',savefig='venus_eye')
bsun.plot2D(ref='Visual',select='proj',savefig='sun_proj')

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))

# Define the parameters to compute light curve and the mesh

lcdep = phoebe.ParameterSet('lcdep',atm='blackbody',ld_coeffs=[0.5],ld_func='linear')
mesh = phoebe.ParameterSet('mesh:marching',delta=0.05,alg='python')

# We only consider the primary to be misaligned. For reference, we also make a
# version of the primary that *is* aligned.

star1a = phoebe.MisalignedBinaryRocheStar(comp1,mesh=mesh,orbit=orbit,pbdep=[lcdep])
star1b = phoebe.BinaryRocheStar(comp1,mesh=mesh,orbit=orbit,pbdep=[lcdep])
star2 = phoebe.BinaryRocheStar(comp2,mesh=mesh,orbit=orbit,pbdep=[lcdep])

# Thus we have two systems, the first with a misaligned primary, the second
# with two stars that are fully aligned.

system1 = phoebe.BodyBag([star1a,star2])
system2 = phoebe.BodyBag([star1b,star2])

# We want to observe the star during one full orbital period
times = np.linspace(0., orbit['period'], 100)[:-1]

# Observe the reference system, compute the light curve and make some check images.
phoebe.observe(system2,times,lc=True,extra_func=[phoebe.observatory.ef_binary_image,
                                                 phoebe.observatory.ef_binary_image],
                                     extra_func_kwargs=[dict(select='rv',name='rv00',),
                                                        dict(select='teff',name='teff00',vmin=10000.,vmax=11300.)])

# Keep track of the computed light curves in a nice figure:
colors = itertools.cycle([plt.cm.spectral(i) for i in np.linspace(0,1,5)])
plotting.plot_lcsyn(system2,'x-',color=colors.next(),ms=10,mew=2,scale=None,label='Aligned')
plt.xlabel("Time [d]")
plt.ylabel("Flux [erg/s/cm2/AA]")

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',
                    lc=True,rv=True,mpi=mpi)

system.set_time(0.)
L1 = phoebe.universe.luminosity(star1)
L2 = phoebe.universe.luminosity(star2)
print("Luminosity ratio Phoebe 2.0 = {:.6f}".format(L2/L1))
print("Luminosity ratio WD         = {:.6f}".format(10**(params['Mbol2']/-2.5)/10**(params['Mbol1']/-2.5)))
print("WD passband luminosity 1 = {:.6f}".format(params['L1']))
print("WD passband luminosity 2 = {:.6f}".format(params['L2']))

c3 = time.time()         

# Analysis of results
# -------------------
# Retrieve the computed fluxes and radial velocities

# 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:
#-----------------------

# Retrieve the results of the Kepler light curve computations, and convert all
# arrays in the parameterSet to arrays (instead of lists), for easy computation.

# Create the bodies.
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 and radial velocity curve and image
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()

# Compute pyphoebe light curve and radial velocity curve:
P = binary['period']
times = np.linspace(binary['t0'],binary['t0']+P,len(curve['indeps']))
bbag = phoebe.BodyBag([star1,star2], globals=globals)


phoebe.observe(bbag,times,subdiv_num=2,lc=True,rv=True,eclipse_alg='binary')#,extra_func=[phoebe.observatory.ef_binary_image],extra_func_kwargs=[dict(select='rv')])
c3 = time.time()
bbag.save('wdvspyphoebe.phoebe')

# Analysis of results
# -------------------
# Access the computed quantities
flux1 = np.array(star1.params['syn']['lcsyn'].values()[0]['flux'])
flux2 = np.array(star2.params['syn']['lcsyn'].values()[0]['flux'])
flux = flux1+flux2
rv1 = np.array(star1.params['syn']['rvsyn'].values()[0]['rv'])
rv2 = np.array(star2.params['syn']['rvsyn'].values()[0]['rv'])
rv1 = phoebe.convert('Rsol/d','km/s',rv1)
rv2 = phoebe.convert('Rsol/d','km/s',rv2)

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

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
is nothing more than a subclass of ``BodyBag``, but which is a little more