How to use the phoebe.ParameterSet function in phoebe

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

from matplotlib import pyplot as plt
import numpy as np
import phoebe
from phoebe.backend import observatory
from phoebe.utils import plotlib

logger = phoebe.get_basic_logger(clevel='INFO')
c0 = time.time()

# Parameter preparation
# ---------------------
incl = 75.,'deg'
long = 30.,'deg'
    
# Create a ParameterSet representing a pulsating star
pstar = phoebe.ParameterSet(frame='phoebe',context='star',add_constraints=True)
pstar['teff'] = 10000.
pstar['incl'] = incl
pstar['long'] = long
pstar['radius'] = 4.5,'Rsol'
pstar['mass'] = 3.0,'Msol'
pstar['rotperiod'] = 10.,'d'
pstar['shape'] = 'sphere'
pstar['atm'] = 'kurucz'
pstar['ld_coeffs'] = 'kurucz'
pstar['ld_func'] = 'claret'
pstar['label'] = 'Pulsating Star'
mesh1 = phoebe.ParameterSet(context='mesh:marching',alg='c')
mesh1['delta'] = 0.1

# Create a ParameterSet with Parameters for the pulsation mode
freq_pars1 = phoebe.ParameterSet(frame='phoebe',context='puls',add_constraints=True)

freq_pars1['ampl'] = 0.01
freq_pars1['k'] = 0.0
freq_pars1['l'] = sys.argv[1] if len(sys.argv)>1 else 4
freq_pars1['m'] = sys.argv[2] if len(sys.argv)>2 else 4
freq_pars1['ledoux_coeff'] = 0.1585
freq_pars1['amplteff'] = 0.01
freq_pars1['scheme'] = 'nonrotating'

# Create a ParameterSet with parameters for the light curve
lcdep1 = phoebe.ParameterSet(context='lcdep',ref='light curve')
lcdep1['ld_func'] = 'claret'
lcdep1['ld_coeffs'] = 'kurucz'
lcdep1['atm'] = 'kurucz'

# Create a parameterSet with parameters for the spectra
spdep1 = phoebe.ParameterSet(context='spdep',ref='line profile')
spdep1['ld_func'] = 'claret'
spdep1['ld_coeffs'] = 'kurucz'
spdep1['atm'] = 'kurucz'

name = ""


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

# Compute the light curve and spectra

name = 'pulsating_star_{0}rot_l{1:d}m{2:d}_{3}'.format((np.isinf(lac['rotperiod']) and 'no' or 'w'),freq_pars1['l'],freq_pars1['m'],name)

def extra_func(system,time,i_time):

# Create a ParameterSet with parameters matching Vega. The bolometric atmosphere
# and limb darkening coefficients are set the Kurucz and Claret models.
vega = phoebe.ParameterSet(frame='phoebe', context='star', label='True Vega',
                           add_constraints=True)
vega['teff'] = 8900, 'K'
vega['mass'] = 2.3, 'Msol'
vega['radius'] = 2.26, 'Rsol'
vega['rotperiod'] = 12.5, 'h'
vega['gravb'] = 1.
vega['incl'] = 4.7, 'deg'
vega['ld_func'] = 'claret'
vega['atm'] = 'kurucz'
vega['ld_coeffs'] = 'kurucz'

globs = phoebe.ParameterSet('position')
globs['distance'] = 7.756, 'pc'

# We want to compute interferometric visibilities in the K band:
ifdep1 = phoebe.ParameterSet(frame='phoebe', context='ifdep', ref='Fast rotator')
ifdep1['ld_func'] = 'claret'
ifdep1['ld_coeffs'] = 'kurucz'
ifdep1['atm']= 'kurucz'
ifdep1['passband'] = '2MASS.KS'

mesh = phoebe.ParameterSet(frame='phoebe', context='mesh:marching', alg='c')
mesh['delta'] = 0.1

# Next, wee'll make a Star similar to Vega but that is not a rapid rotator:
# We'll make uniform disk model and power-law limbdarkened disk on a
# non-rotating star. Note that we scale the radius of the stars according to the
# fitted angular diameter of Aufdenberg. That is, we add a constraint on the

star['long'] = 20, 'deg'
star['mass'] = 14.,'Msol'
star['radius'] = 6.4,'Rsol'
star['vgamma'] = 7.632287,'km/s'

# Just for fun, also add the parallax, the surface gravity and vsini (and define
# a mesh).
tools.add_parallax(star,parallax=4.76,unit='mas')
tools.add_surfgrav(star,3.70,derive='radius')
tools.add_vsini(star,32.359421,unit='km/s',derive='incl')

mesh = phoebe.PS('mesh:marching',delta=0.05, alg='c')


# For the parameters of the pulsation mode, we're inspired by [Telting1997]_:
freq1 = phoebe.ParameterSet('puls',freq=  5.24965427519,phase=0.122545,ampl=0.10525/50.,l=3,m=2,amplteff=0.01)

# For the parameters of the magnetic dipole, we take guesses from our own
# research:
mag_field = phoebe.PS('magnetic_field')
mag_field['Bpolar'] = 276.01,'G'
mag_field['beta'] = 61.29,'deg'
mag_field['phi0'] = 80.522,'deg'


# Fake some observations, and specify the passband dependent parameters.
plobs = datasets.PLDataSet(time=np.linspace(0, 0.5*star['rotperiod'], 500),
                           wavelength=np.linspace(399.92,400.08,500), ref='mylsd',
                           columns=['time','wavelength','flux','V','continuum'])

pldep = phoebe.ParameterSet('pldep',ref='mylsd', weak_field=False,
                            atm='kurucz',ld_func='claret',ld_coeffs='kurucz')

# Next, create the parameters of the uniform source. This includes a parameter
# that sets the density of the mesh, and also the radius of the secondary
# sphere, because we will need that in the computation of the analytical flux.
# The second sphere will have a radius half of the primary.    
sphere1 = phoebe.ParameterSet(frame='phoebe',context='star',add_constraints=True)
sphere1['radius'] = 1.,'Rsol'
sphere1['teff'] = 5777.,'K'
sphere1['atm'] = 'kurucz'
sphere1['ld_func'] = 'claret'
sphere1['ld_coeffs'] = 'kurucz'
sphere1['rotperiod'] = 5.,'d'
sphere1['shape'] = 'equipot'
sphere1['incl'] = 0.

mesh1 = phoebe.ParameterSet(frame='phoebe',context='mesh:marching',alg='c')
mesh1['delta'] = 0.05

# To compute light curves, the parameters for to compute the light curve have
# to be given. We compute the light curve analytically and numerically, so two
# sets are necessary.
lcdep1a = phoebe.ParameterSet(frame='phoebe',context='lcdep')
lcdep1a['ld_func'] = 'claret'
lcdep1a['ld_coeffs'] = 'kurucz'
lcdep1a['atm'] = 'kurucz'
lcdep1a['ref'] = 'SDSS.GP'
lcdep1a['passband'] = 'SDSS.GP'

lcdep1b = lcdep1a.copy()
lcdep1b['ref'] = 'SDSS.RP'
lcdep1b['passband'] = 'SDSS.RP'

# to be given. We compute the light curve analytically and numerically, so two
# sets are necessary.
lcdep1a = phoebe.ParameterSet(frame='phoebe',context='lcdep')
lcdep1a['ld_func'] = 'claret'
lcdep1a['ld_coeffs'] = 'kurucz'
lcdep1a['atm'] = 'kurucz'
lcdep1a['ref'] = 'SDSS.GP'
lcdep1a['passband'] = 'SDSS.GP'

lcdep1b = lcdep1a.copy()
lcdep1b['ref'] = 'SDSS.RP'
lcdep1b['passband'] = 'SDSS.RP'

# The secondary sphere is a dark sphere. We set its temperature equal to zero
# Kelvin.
sphere2 = phoebe.ParameterSet(frame='phoebe',context='star',add_constraints=True)
sphere2['radius'] = 1,'Rjup'
sphere2['teff'] = 0,'K'
sphere2['atm'] = 'blackbody'
sphere2['incl'] = 0.
sphere2['rotperiod'] = np.inf
sphere2['ld_func'] = 'uniform'
sphere2['ld_coeffs'] = [1]
sphere2.add_constraint('{radius2} = %.6g'%(sphere1.get_value('radius','m')))

mesh2 = phoebe.ParameterSet(frame='phoebe',context='mesh:marching')
mesh2['delta'] = 0.1

ps = phoebe.ParameterSet(frame='phoebe',context='orbit',add_constraints=True,c1label=sphere1['label'],c2label=sphere2['label'])
ps['sma'] = 5.

# Body setup