Fluxes jitter at ~1e-4 level

[aprsa]

The bottom fig shows elv (blue), reflection (green) and boosting (red); elv jitters because of mesh effects.

[aprsa]

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

[horvatm]

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?

[aprsa]

log Period. See the paper for the complete version of this figure.

[horvatm]

Can you give me a complete script?

[aprsa]

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