PHOEBE 2.1 Documentation

2.1 Docs

IPython Notebook | Python Script

Sun (single rotating star)


As always, let’s do imports and initialize a logger and a new bundle. See Building a System for more details.

%matplotlib inline
import phoebe
from phoebe import u # units
import numpy as np
import matplotlib.pyplot as plt

logger = phoebe.logger()

b = phoebe.default_star(starA='sun')

Setting Parameters

print b['sun']
ParameterSet: 22 parameters
            requiv@sun@component: 1.0 solRad
        requiv_max@sun@component: 1.0 solRad
              teff@sun@component: 6000.0 K
              abun@sun@component: 0.0
            period@sun@component: 1.0 d
*             freq@sun@component: 6.283185 rad / d
             pitch@sun@component: 0.0 deg
               yaw@sun@component: 0.0 deg
              incl@sun@component: 90.0 deg
           long_an@sun@component: 0.0 deg
         gravb_bol@sun@component: 0.32
  irrad_frac_refl_bol@sun@com...: 0.6
* irrad_frac_lost_bol@sun@com...: 0.4
       ld_func_bol@sun@component: logarithmic
     ld_coeffs_bol@sun@component: [0.5 0.5]
              mass@sun@component: 1.0 solMass
             freq@sun@constraint: 6.283185 / {period@sun@component}
  irrad_frac_lost_bol@sun@con...: 1.000000 - {irrad_frac_refl_bol@sun@component}
  mesh_method@sun@phoebe01@co...: marching
  ntriangles@sun@phoebe01@com...: 1500
  distortion_method@sun@phoeb...: rotstar
        atm@sun@phoebe01@compute: ck2004

Let’s set all the values of the sun based on the nominal solar values provided in the units package.

b.set_value('teff', 1.0*u.solTeff)
b.set_value('requiv', 1.0*u.solRad)
b.set_value('mass', 1.0*u.solMass)
b.set_value('period', 24.47*u.d)

And so that we can compare with measured/expected values, we’ll observe the sun from the earth - with an inclination of 23.5 degrees and at a distance of 1 AU.

b.set_value('incl', 23.5*u.deg)
b.set_value('distance', 1.0*u.AU)

Checking on the set values, we can see the values were converted correctly to PHOEBE’s internal units.

print b.get_quantity('teff')
print b.get_quantity('requiv')
print b.get_quantity('mass')
print b.get_quantity('period')
print b.get_quantity('incl')
print b.get_quantity('distance')
5772.0 K
1.0 solRad
1.0 solMass
24.47 d
23.5 deg
1.495978707e+11 m

Running Compute

Let’s add a light curve so that we can compute the flux at a single time and compare it to the expected value. We’ll set the passband luminosity to be the nominal value for the sun. We’ll also add a mesh dataset so that we can plot the temperature distributions and test the size of the sun verse known values.

b.add_dataset('lc', times=[0.], pblum=1*u.solLum)
b.add_dataset('mesh', times=[0.], columns=['teffs', 'loggs', 'rs'])
<ParameterSet: 4 parameters | contexts: compute, dataset>
b.run_compute(irrad_method='none', distortion_method='rotstar')
<ParameterSet: 8 parameters | kinds: mesh, lc>

Comparing to Expected Values

afig, mplfig = b['mesh'].plot(fc='teffs', x='xs', y='ys', show=True)
afig, mplfig = b['mesh'].plot(fc='teffs', x='us', y='vs', show=True)
print "teff: {} ({})".format(b.get_value('teffs').mean(),
                             b.get_value('teff', context='component'))
teff: 5771.99999826 (5772.0)

For a rotating sphere, the minimum radius should occur at the pole and the maximum should occur at the equator.

print "rmin (pole): {} ({})".format(b.get_value('rs').min(),
                             b.get_value('requiv', context='component'))
rmin (pole): 0.999999812391 (1.0)
print "rmax (equator): {} (>{})".format(b.get_value('rs').max(),
                              b.get_value('requiv', context='component'))
rmax (equator): 1.00000009416 (>1.0)
print "logg: {}".format(b.get_value('loggs').mean())
logg: 4.43806746133
print "flux: {}".format(b.get_quantity('fluxes@model')[0])
flux: 1358.70971154 W / m2
Prev: Intensity Weighting Next: Minimal Example to Produce a Synthetic Light Curve
Last update: 10/29/2018 9:20 a.m. (CET)