Open aprsa opened 7 years ago
Parameters:
b = phoebe.default_binary()
b.set_value('rpole', component='primary', value=1.8) b.set_value('rpole', component='secondary', value=0.96)
b.set_value('teff', component='primary', value=10000) # A0 b.set_value('teff', component='secondary', value=5200) # K0
b.set_value('frac_refl_bol', component='primary', value=1.0) b.set_value('frac_refl_bol', component='secondary', value=0.6)
b.set_value('q', component='binary', value=0.96/1.8) b.set_value('incl', component='binary', value=88)
b.add_dataset('lc', passband='Johnson:V', intens_weighting='energy', dataset='data')
This perhaps could be fixed a more proper averaging produce, where the mesh is truly randomized. Btw, on y-axis is the intensity and what is labeled on x-axis?
log Period. See the paper for the complete version of this figure.
Can you give me a complete script?
import phoebe
import numpy as np
import matplotlib.pyplot as plt
GMSun = 1.3271244e20
RSun = 6.957e8
day = 86400.
def sma(M1, q, P):
return (GMSun*day**2*M1*(1+q)*P**2/(4*np.pi**2))**(1./3)/RSun
naver = 10
b = phoebe.default_binary()
b.set_value('rpole', component='primary', value=1.8)
b.set_value('rpole', component='secondary', value=0.96)
b.set_value('teff', component='primary', value=10000) # A0
b.set_value('teff', component='secondary', value=5200) # K0
b.set_value('frac_refl_bol', component='primary', value=1.0)
b.set_value('frac_refl_bol', component='secondary', value=0.6)
b.set_value('q', component='binary', value=0.96/1.8)
b.set_value('incl', component='binary', value=88)
b.add_dataset('lc', passband='Johnson:V', intens_weighting='energy', dataset='data')
b.add_compute(compute='elv_only', distortion_method='roche', reflection_method='none', boosting_method='none')
logps = np.linspace(0, 2, 100)
elv_amp = []
elv_fluxes = np.empty(naver)
fout = open('effect_amplitudes.res', 'w')
for logp in logps:
period = 10**logp
b.set_value('sma', component='binary', value=sma(1.8, 0.96/1.8, period))
# set the period and timestamp that corresponds to quarter-phase
b.set_value('period', component='binary', value=period)
b.set_value('times', dataset='data', context='dataset', value=[0.25*period])
# ELVs need to be averaged into a smooth curve
for i in range(naver):
b.set_value('mesh_init_phi', compute='elv_only', value=2*np.pi*np.random.rand())
b.run_compute(compute='elv_only', model='elv_model')
elv_fluxes[i] = b.get_value('fluxes', model='elv_model')[0]
elv_only_flux = np.mean(elv_fluxes)
elv_amp.append(elv_only_flux)
fout.write('%f %f %f\n' % (logp, sma(1.8, 0.96/1.8, period), elv_amp[-1]))
print logp, sma(1.8, 0.96/1.8, period), elv_amp[-1]
fout.close()
plt.plot(logps, elv_amp, 'b-', lw=2, label='ELV')
plt.show()
The bottom fig shows elv (blue), reflection (green) and boosting (red); elv jitters because of mesh effects.