Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
sun_mesh = phoebe.ParameterSet(context='mesh:marching',alg='c')
sun_mesh['delta'] = 0.05
lcdep1 = phoebe.ParameterSet(frame='phoebe',context='lcdep')
lcdep1['ld_func'] = 'claret'
lcdep1['ld_coeffs'] = 'kurucz'
lcdep1['atm'] = 'kurucz'
lcdep1['passband'] = 'OPEN.BOL'
lcdep1['ref'] = 'Bolometric (numerical)'
lcdep2 = lcdep1.copy()
lcdep2['method'] = 'analytical'
lcdep2['ref'] = 'Bolometric (analytical)'
the_sun = phoebe.Star(sun,sun_mesh,pbdep=[lcdep1,lcdep2], position=globals)
the_sun.set_time(0)
the_sun.lc()
params = the_sun.get_parameters()
nflux = the_sun.params['syn']['lcsyn']['Bolometric (numerical)']['flux'][0]
aflux = the_sun.params['syn']['lcsyn']['Bolometric (analytical)']['flux'][0]
mupos = the_sun.mesh['mu']>0
num_error_area = np.abs(np.pi-((the_sun.mesh['size']*the_sun.mesh['mu'])[mupos]).sum())/np.pi*100
num_error_flux = np.abs(nflux-aflux)/aflux*100
print aflux, nflux, the_sun.projected_intensity()
real_error_flux = np.abs(1368.000-aflux)/aflux*100
# them differently, since the ifdeps have different references.
ifobs_simul2 = ifobs_simul1.copy()
ifobs_simul3 = ifobs_simul1.copy()
ifobs_simul4 = ifobs_simul1.copy()
ifobs_simul1['ref'] = ifdep1['ref']
ifobs_simul2['ref'] = ifdep2['ref']
ifobs_simul3['ref'] = ifdep3['ref']
ifobs_simul4['ref'] = ifdep4['ref']
# Body setup
# ----------
# Build the Star Bodies
star1 = phoebe.Star(vega, mesh, pbdep=ifdep1, obs=ifobs_simul1, position=globs)
star2 = phoebe.Star(vega2, mesh, pbdep=ifdep2, obs=ifobs_simul2, position=globs)
star3 = phoebe.Star(vega3, mesh, pbdep=ifdep3, obs=ifobs_simul3, position=globs)
star4 = phoebe.Star(vega4, mesh, pbdep=ifdep4, obs=ifobs_simul4, position=globs)
# For convenience, we put everything in a BodyBag. Then we can set the time of
# all bodies simultaneously and don't have to cycle through them.
system = phoebe.BodyBag([star1, star2, star3, star4])
# Computation of observables
# --------------------------
system.compute(eclipse_alg='only_horizon')
# Analysis of results
# -------------------
# Make some check images
phoebe.image(star1,savefig='vega_image1.png',ref=0,context='ifdep')
lcdep2['method'] = 'analytical'
lcdep2['ref'] = 'Bolometric (analytical)'
# Body setup
# ----------
# Next, a Body needs to be created in the Universe. The Body best representing
# the Sun is of course a single :py:class:`Star `.
# The parameters for the *pbdeps* need
# to be passed along. The parameters of the Body itself can always be accessed
# via ``the_sun.params['star']``, the parameters of the *pbdeps* can be
# accessed via ``the_sun.params['pbdep']``, which is an ordered dictionary.
# Thus, the light curves can be accessed via
# ``the_sun.params['pbdep']['lcdep'].values()[0]``, and the numerical light
# curve on its turn via ``the_sun.params['pbdep']['lcdep'].values()[1]``.
the_sun = phoebe.Star(sun, mesh=sun_mesh, pbdep=[lcdep1, lcdep2], position=globs)
# Computation of observables
# --------------------------
# Upon creation, the Body does not yet exist in the Universe. We need to set the
# time of the object via :py:func:`set_time() `, so that it knows where and in what orientiation it needs
# to put itself. For a single, non-rotating star this is easy; it will be put
# in the origin of the Universe but inclined with respect to the line of sight.
# The value of the time is in this case unimportant, since it is
# time-independent.
the_sun.set_time(0)
# Now compute the observables.
the_sun.lc()
"""
And finally make some check images using :py:class:`plot2D() `:
freq_pars2['m'] = 1
freq_pars2['ampl'] = 0.05
freq_pars3 = freq_pars2.copy()
freq_pars3['m'] = -1
# Create a ParameterSet with parameters for the light curve
lcdep1 = phoebe.ParameterSet(frame='phoebe',context='lcdep',ref='lc')
lcdep1['ld_func'] = 'claret'
lcdep1['ld_coeffs'] = 'kurucz_p00'
lcdep1['atm'] = 'kurucz_p00'
# Body setup
# ----------
star = phoebe.Star(lac,mesh,puls=[freq_pars1,freq_pars2,freq_pars3],pbdep=[lcdep1])
# Computation of observables
# --------------------------
n = 400
times = np.linspace(0,17,n)*lac['rotperiod']
rotperiods = np.ones(n)
rotperiods[:n/2] = np.inf
rotperiods[n/2:] = lac['rotperiod']
for i,itime in enumerate(times):
if i==0: continue
print(i)
star.reset()
star.params['star']['rotperiod'] = rotperiods[i]
star.set_time(itime,ref='lc')
plt.savefig('{}_teff_{:05d}.png'.format(name,i_time));plt.close()
xlim,ylim,p = phoebe.image(system,ref='light curve',select='proj2',vmax=1.2e5)
plt.xlim(-9,9);plt.ylim(-9,9)
plt.savefig('{}_proj2_{:05d}.png'.format(name,i_time));plt.close()
# Make sure to always have the same time period, so that we can see the frequency
# shifts
times = np.linspace(0,1./(freq_pars1['freq']-4./lac.get_value('rotperiod','d')),50)
# Create a DataSet with the template spectra and light curve
spobs1 = phoebe.SPDataSet(wavelength=np.linspace(4544, 4556, 200),
time=times, ref='line profile')
lcobs1 = phoebe.LCDataSet(time=times, ref='light curve')
star = phoebe.Star(lac,mesh,puls=[freq_pars1],
pbdep=[lcdep1,spdep1], obs=[lcobs1, spobs1])
phoebe.compute(star, subdiv_num=0, extra_func=[extra_func])
"""
For synthetic Gaussian line profile templates:
+--------------------------------------------------------+--------------------------------------------------------------+--------------------------------------------------------+--------------------------------------------------------------+
| l=4, m=4: line profile variations | l=4, m=4 effective temperature | l=4, m=4: radial velocity | l=4, m=4 intensity |
+--------------------------------------------------------+--------------------------------------------------------------+--------------------------------------------------------+--------------------------------------------------------------+
| .. image:: images_tut/pulsating_star_wrot_l4m4_.gif | .. image:: images_tut/pulsating_star_wrot_l4m4__teff.gif | .. image:: images_tut/pulsating_star_wrot_l4m4__rv.gif | .. image:: images_tut/pulsating_star_wrot_l4m4__proj2.gif |
| :scale: 100 % | :scale: 100 % | :scale: 100 % | :scale: 100 % |
| :width: 233px | :width: 233px | :width: 233px | :width: 233px |
| :height: 233px | :height: 233px | :height: 233px | :height: 233px |
+--------------------------------------------------------+--------------------------------------------------------------+--------------------------------------------------------+--------------------------------------------------------------+
# 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
# --------------------------
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.