How to use phoebe - 10 common examples

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

b.add_solver('estimator.rv_periodogram')
    b.add_solver('estimator.lc_geometry')
    b.add_solver('estimator.rv_geometry')
    b.add_solver('estimator.ebai')
    b.add_solver('optimizer.nelder_mead')
    b.add_solver('optimizer.differential_evolution')
    b.add_solver('sampler.emcee')
    b.add_solver('sampler.dynesty')


    # TODO: include constraint_func?  Shouldn't matter since they're not in twigs
    should_be_forbidden = b.qualifiers + b.contexts + b.kinds + [c.split('@')[0] for c in b.get_parameter('columns').choices]

    if verbose:
        for l in should_be_forbidden:
            if l not in phoebe.parameters.parameters._forbidden_labels:
                print(l)

    for l in should_be_forbidden:
        assert(l in phoebe.parameters.parameters._forbidden_labels)

    phoebe.reset_settings()
    return b

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

twigs_atm = mybundle.search('atm')
    for atm in twigs_atm:
        if atm.split('@')[0] == 'value':
            continue
        if mybundle[atm] in ['blackbody','kurucz']:
            mybundle[atm] = 'blackbody_uniform_none_teff.fits'
            passband_twig = 'passband@{}'.format("@".join(atm.split('@')[1:]))
            if passband_twig in mybundle and mybundle[passband_twig] == 'JOHNSON.V':
                mybundle[passband_twig] = 'johnson_v.ptf'
    
    mybundle['pblum@secondary'] = phb.getpar('phoebe_plum2')
    mybundle.run_compute(label='from_legacy', irradiation_alg='point_source')
    #mybundle.get_system().compute(animate=True)
    lc_ph2 = mybundle['flux@{}@lcsyn'.format(dataref)]
    
    U = phb2.units.conversions.Unit
    R1 = U(mybundle.get_system()[0].params['component'].request_value('r_pole'), 'm')
    T1 = U(mybundle['teff@primary'], 'K')
    R2 = U(mybundle.get_system()[1].params['component'].request_value('r_pole'), 'm')
    T2 = U(mybundle['teff@secondary'], 'K')
    sigma = U('sigma')
    L1 = (4*np.pi*R1**2*sigma*T1**4).convert('Lsol')
    L2 = (4*np.pi*R2**2*sigma*T2**4).convert('Lsol')
    print("Numerical bolometric luminosity (primary) = {} Lsol".format(phb2.convert('erg/s', 'Lsol',mybundle['primary'].luminosity())))
    print("Numerical bolometric luminosity (secondary) = {} Lsol".format(phb2.convert('erg/s', 'Lsol',mybundle['secondary'].luminosity())))
    print("Eq. sphere bolometric luminosity (primary) = {}".format(L1))
    print("Eq. sphere bolometric luminosity (secondary) = {}".format(L2))
    print("Numerical passband luminosity (primary) = {} Lsol".format(phb2.convert('erg/s', 'Lsol',mybundle['primary'].luminosity(ref='LC'))))
    print("Numerical passband luminosity (secondary) = {} Lsol".format(phb2.convert('erg/s', 'Lsol',mybundle['secondary'].luminosity(ref='LC'))))
    
    
    # Passband luminosities:

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

    #-- compute the system if the light curves haven't been computed before
    if not os.path.isfile(phoebe_lc_file) or recompute:
        #-- compute the system
        #time = time[::5]

ref='mylc')
    
    
    # create bodies
    starA = universe.BinaryRocheStar(Apars, mesh=meshpars, orbit=orbitpars,
                                    pbdep=[pbdep_lc,pbdep_rv_a], obs=[obs_rv_a])
    starB = universe.BinaryRocheStar(Bpars, mesh=meshpars, orbit=orbitpars,
                                    pbdep=[pbdep_lc,pbdep_rv_b], obs=[obs_rv_b])
    system = universe.BinaryBag([starA,starB], orbit=orbitpars, label='V380Cyg',
                                position=globs,obs=[obs_lc])
    
    # compute
    #obs_lc.set_adjust('offset',True)
    #obs_lc.set_adjust('scale',True)
    tools.group([obs_rv_a, obs_rv_b], 'rv', scale=False, offset=True)
    compute_options = parameters.ParameterSet(context='compute')
    system.compute(params=compute_options, mpi=None)

    lcobs = system.get_obs(category='lc', ref='mylc').asarray()
    lcsyn = system.get_synthetic(category='lc', ref='mylc').asarray()

    rvobs1 = system[0].get_obs(category='rv', ref='myrv_comp1').asarray()
    rvsyn1 = system[0].get_synthetic(category='rv', ref='myrv_comp1').asarray()

    rvobs2 = system[1].get_obs(category='rv', ref='myrv_comp2').asarray()
    rvsyn2 = system[1].get_synthetic(category='rv', ref='myrv_comp2').asarray()

    if not debug:
        assert(np.mean((lcobs['flux']-(lcsyn['flux']*lcobs['scale'] - lcobs['offset']))**2/lcobs['sigma']**2)<0.00017)
        assert(np.mean((rvobs1['rv']-(rvsyn1['rv'] - rvobs1['offset']))**2/rvobs1['sigma']**2)<2.14)
        assert(np.mean((rvobs2['rv']-(rvsyn2['rv'] - rvobs2['offset']))**2/rvobs2['sigma']**2)<2.41)
    else:

import time
import numpy as np
from matplotlib import pyplot as plt
import phoebe

logger = phoebe.get_basic_logger()

c0 = time.time()

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

# Create a :ref:`star ParameterSet ` with parameters
# matching the Sun, but make a fine-enough mesh. Also, the rotation period is
# set to almost 90% of the Sun's critical rotation period.
star = phoebe.ParameterSet(context='star')
star['atm'] = 'blackbody'
star['ld_func'] = 'linear'
star['ld_coeffs'] = [0.6]
star['rotperiod'] = 0.24,'d'
star['shape'] = 'sphere'
star['teff'] = 10000.
star['radius'] = 1.0,'Rsol'

mesh = phoebe.ParameterSet(context='mesh:marching')
mesh['delta'] = 0.05

# These are the parameters for the calculation of the :ref:`spectrum `. We assume the
# spectral line is in the visual region, and use a linear limb darkening law.
spdep1 = phoebe.ParameterSet(context='spdep')
spdep1['ld_func'] = 'linear'
spdep1['atm'] = 'blackbody'

spdep1 = phoebe.ParameterSet(context='spdep')
    spdep1['ld_func'] = 'linear'
    spdep1['atm'] = 'blackbody'
    spdep1['ld_coeffs'] = [0.6]
    spdep1['passband'] = 'JOHNSON.V'
    spdep1['method'] = 'numerical'
    spdep1['ref'] = 'Numerical'

    spdep2 = spdep1.copy()
    spdep2['method'] = 'analytical'
    spdep2['ref'] = 'Via convolution'

    wavelengths = np.linspace(399.7, 400.3, 1000)
    spobs1 = phoebe.ParameterSet(context='spobs', ref=spdep1['ref'], wavelength=wavelengths)
    spobs2 = phoebe.ParameterSet(context='spobs', ref=spdep2['ref'], wavelength=wavelengths)

    mesh1 = phoebe.Star(star, mesh, pbdep=[spdep1, spdep2])

    mesh1.set_time(0)

    mesh1.sp(obs=spobs1)
    mesh1.sp(obs=spobs2)

    result1 = mesh1.get_synthetic(category='sp',ref=0)
    result2 = mesh1.get_synthetic(category='sp',ref=1)

    flux1 = np.array(result1['flux'][0])/np.array(result1['continuum'][0])
    flux2 = np.array(result2['flux'][0])/np.array(result2['continuum'][0])

    if not from_main:
        assert (np.all(np.abs((flux1-flux2)/flux1)<=0.00061))

#-- setup MPI stuff: check which host the script is running on. If it is
    #   clusty and we can readily find a hosts file, we fill in the MPI parameter
    #   set.
    hostname = socket.gethostname()
    if hostname=='clusty':
        hostfile = os.path.expanduser('~/mpi/hosts')
        print("Running on Clusty, trying to load hostfile {}".format(hostfile)) 
        if not os.path.isfile(hostfile):
            print("Cannot load hostfile {} (file does not exist), falling back to defaults".format(hostfile))
        else:
            mpi = phoebe.ParameterSet(context='mpi',np=24,hostfile=hostfile,
                                        byslot=True,python='python2.7')
    #-- If the hostname isn't clusty, we're probably running on a normal computer
    #   in that case, just take the number of available processes.
    else:
        mpi = phoebe.ParameterSet(context='mpi',np=multiprocessing.cpu_count())
    return mpi

The priors are set for each `real` parameter: we set the inclination and
potential values to have a normal distribution, and the eccentricity to have
a uniform distribution.
"""
pset.get_parameter('incl').set_prior(distribution='normal',mu=90.,sigma=5.)
pset.get_parameter('ecc').set_prior(distribution='uniform',lower=0.0,upper=1.0)
pset.get_parameter('pot1').set_prior(distribution='normal',mu=6.2,sigma=0.2)

# You might wonder what the purpose is for setting priors for parameters if we
# want to use a regular nonlinear fitting routine. The reason is twofold: one
# is to tell the code we want to fit it, and the other is that the nonlinear
# fitting routines accept boundaries for the values. If you don't want to use
# these boundaries (they may influence the error determination considerably), you
# need to set ``bounded=False`` in the parameterSet controlling the properties
# of the fit:
fitparams = phoebe.ParameterSet(context='fitting:lmfit',method='leastsq',bounded=True, iters=1)

"""

Printing this parameterSet reveals the following parameters::
    
    print(fitparams)
        method leastsq                               --   phoebe Nonlinear fitting method
         iters 0                                     --   phoebe Number of iterations
         label c8b61573-a686-4886-aef9-d2d0a3846c93  --   phoebe Fit run name
    compute_ci False                                 --   phoebe Compute detailed confidence intervals
       bounded True                                  --   phoebe Include boundaries in fit
      feedback {:}                                   --   phoebe Results from fitting procedure

The type of nonlinear fitting algorithm is given by ``method``, and you can
give an additional label for easy reference. If you set ``iters=0``, the algorithm
will be run just once starting from the current values (i.e. the once we've