How to use the phoebe.PS function in phoebe

# First, import necessary modules
import matplotlib.pyplot as plt
import phoebe
from phoebe.parameters import tools
from phoebe.atmospheres import passbands

logger = phoebe.get_basic_logger()

# Parameter preparation
# ----------------------

# We load the parameters of Vega from the library, but add a few new ones
# which we like to fit.

star = phoebe.create.from_library('Vega', gravblaw='espinosa')
mesh = phoebe.PS(context='mesh:marching',alg='c',delta=0.1)

# The parameters that we want to fit but are not available, are not very
# exotic so they can be added via the :py:mod:`parameters.tools `
# module

tools.add_surfgrav(star, 4.0, derive='mass')
tools.add_rotfreqcrit(star, 0.5)
tools.add_teffpolar(star, 9000.)

# Next, we tell to code to consider these parameters to be fitted.
star.set_adjust(('teffpolar','radius','surfgrav','incl','rotfreqcrit'),True)

# To be sure the code doesn't go of the grid, we define priors, which will be
# interpreted as boundaries for the Levenberg-Marquardt algorithm.
star.get_parameter('teffpolar').set_prior(distribution='uniform',lower=4500,upper=15000.)
star.get_parameter('surfgrav').set_prior(distribution='uniform',lower=3.,upper=5.0)

# Create the star
    star = phoebe.create.from_library('vega', create_body=True)
    star.params['position']['epoch'] = 'J1991.25'
    star.params['position']['pmra'] = pmra
    star.params['position']['pmdec'] = pmdec
    star.params['position']['ra'] = ra
    star.params['position']['dec'] = dec
    star.params['position']['distance'] = distance
    
    # Read in the Hipparcos positions
    times, x_earth, y_earth, z_earth = np.genfromtxt(os.path.join(basedir,'hipparcos_position.dat'),
                     delimiter=',', skip_header=53, skip_footer=32,
                     usecols=(0,2,3,4), unpack=True)
    
    # Create the dep and obs
    amdep = phoebe.PS(context='amdep', ref='myam')
    amobs = phoebe.DataSet('amobs', ref='myam', time=times, eclx=x_earth,
                       ecly=y_earth, eclz=z_earth)[::100]

    # Add the dep and obs, and compute
    star.add_pbdeps([amdep])
    star.add_obs([amobs])
    star.init_mesh()
    star.compute()

    # Get the synthetic stuff
    output = star.get_synthetic(category='am',ref='myam').asarray()
    return output['delta_ra']*3600*1000, output['delta_dec']*3600*1000

# in case the Kurucz grid does not cover a wide enough parameter space.
# First the case with doppler boosting:
comp1, comp2, orbit = create.from_library('KPD1946+4340_bis')
mesh = phoebe.PS('mesh:marching')
lcdep1 = phoebe.PS('lcdep', atm=comp1['atm'], ld_coeffs=comp1['ld_coeffs'],
                  ld_func=comp1['ld_func'], passband='KEPLER.V', boosting=True)
lcdep2 = phoebe.PS('lcdep', atm=comp2['atm'], ld_coeffs=comp2['ld_coeffs'],
                  ld_func=comp2['ld_func'], passband='KEPLER.V', boosting=True)

star1 = phoebe.BinaryRocheStar(comp1, mesh=mesh, orbit=orbit, pbdep=[lcdep1])
star2 = phoebe.BinaryRocheStar(comp2, mesh=mesh, orbit=orbit, pbdep=[lcdep2])
system1 = phoebe.BodyBag([star1, star2])

# Then without doppler boosting:
comp1, comp2, orbit = create.from_library('KPD1946+4340_bis')
lcdep1 = phoebe.PS('lcdep', atm=comp1['atm'], ld_coeffs=comp1['ld_coeffs'],
                  ld_func=comp1['ld_func'], passband='KEPLER.V', boosting=False)
lcdep2 = phoebe.PS('lcdep', atm=comp2['atm'], ld_coeffs=comp2['ld_coeffs'],
                  ld_func=comp2['ld_func'], passband='KEPLER.V', boosting=False)
star1 = phoebe.BinaryRocheStar(comp1, mesh=mesh, orbit=orbit, pbdep=[lcdep1])
star2 = phoebe.BinaryRocheStar(comp2, mesh=mesh, orbit=orbit, pbdep=[lcdep2])
system2 = phoebe.BodyBag([star1, star2])

# Computation of observables
# --------------------------
# Generate a time series to simulate the light curve and images on, and compute
# them:
period = system1[0].params['orbit']['period']
t0 = system1[0].params['orbit']['t0']
times = np.linspace(t0,t0+period,150)

phoebe.observe(system1,times,lc=True,heating=False,refl=False,

# Technically speaking, we don't need light curve information. However we'll
# use it anyway, because then we can readily call the ``compute`` function,
# which will take care of calling ``set_time`` and the radiation budget 
# functionality in the right way (since these need to be computed to get the
# light curve anyway). We could do it manually, but this is easier.

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

# Choose a moderate mesh density for the Sun (not so important), but a finer
# grid for Mercury to resolve the latitudes better.

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

# Take care of output to the screen
logger = phoebe.get_basic_logger(clevel='INFO')

# Create the two stars and put them in a binary orbit
star1 = phoebe.PS('star', teff=30000, ld_func='linear', ld_coeffs=[0.5], radius=(1., 'Rsol'))
star2 = phoebe.PS('star', teff=1000, ld_func='linear', ld_coeffs=[0.5], radius=(2.,'Rsol'), alb=1)

comp1, comp2, orbit = phoebe.create.binary_from_stars(star1, star2, period=(1,'d'))
orbit['incl'] = 85.0

# Create the mesh and the light curves
mesh1 = phoebe.ParameterSet(context='mesh:marching', delta=0.075)
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

def do_mercury(redist):
    sun, position = phoebe.create.from_library('sun')
    sun['irradiator'] = True


    mercury = phoebe.PS('star', teff=100., atm=atm, radius=(0.3829,'Rearth'),
                    mass=(0.055,'Mearth'), ld_func='uniform', shape='sphere',
                    rotperiod=(58.646,'d'),alb=alb, label='mercury', redist=redist,
                    incl=90., long=0, irradiator=False)

    orbit = phoebe.PS('orbit', period=87.9691, t0=0, t0type='superior conjunction',
                    sma=(0.307499,'au'), q=mercury['mass']/sun['mass'], ecc=0.205630,
                    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])

rotperiod=(58.646,'d'),alb=alb, label='mercury', redist=redist,
                    incl=90., long=0, irradiator=False)

    orbit = phoebe.PS('orbit', period=87.9691, t0=0, t0type='superior conjunction',
                    sma=(0.307499,'au'), q=mercury['mass']/sun['mass'], ecc=0.205630,
                    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

def do_mercury(redist):
    sun, position = phoebe.create.from_library('sun')
    sun['irradiator'] = True


    mercury = phoebe.PS('star', teff=100., atm=atm, radius=(0.3829,'Rearth'),
                    mass=(0.055,'Mearth'), ld_func='uniform', shape='sphere',
                    rotperiod=(58.646,'d'),alb=alb, label='mercury', redist=redist,
                    incl=90., long=0, irradiator=False)

    orbit = phoebe.PS('orbit', period=87.9691, t0=0, t0type='superior conjunction',
                    sma=(0.307499,'au'), q=mercury['mass']/sun['mass'], ecc=0.205630,
                    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