#SimpleBinary: fitting
import phoebe
import numpy as np
import matplotlib.pyplot as plt

logger = phoebe.get_basic_logger()
star1 = phoebe.ParameterSet('star', mass=1.2, radius=1.1, teff=6321.,
                            atm='kurucz', ld_coeffs='kurucz', ld_func='claret')
star2 = phoebe.ParameterSet('star', mass=0.8, radius=0.9, teff=4123.,
                            atm='kurucz', ld_coeffs='kurucz', ld_func='claret')
comp1, comp2, orbit = phoebe.create.binary_from_stars(star1, star2, period=(10.,'d'))
comp1['syncpar'] = 1.0
comp2['syncpar'] = 1.0
orbit['ecc'] = 0.1
mesh = phoebe.ParameterSet('mesh:marching', delta=0.07)
lcdep = phoebe.ParameterSet('lcdep', passband='KEPLER.V', atm='kurucz',
                            ld_func='claret', ld_coeffs='kurucz', ref='my kepler lc')
rvdep = phoebe.ParameterSet('rvdep', passband='JOHNSON.V', atm='kurucz',
                            ld_func='claret', ld_coeffs='kurucz')
#loading observations
input_lc, lcdep_=phoebe.parse_lc('lc.obs')
input_lc['ref']='my kepler lc'
input_lc.get_parameter('pblum').set_adjust(True)
input_lc.get_parameter('l3').set_adjust(True)
phoebe.parameters.tools.add_esinw_ecosw(orbit)
orbit['incl'] = 88.
incl = orbit.get_parameter('incl')
incl.set_adjust(True)
incl.set_prior(distribution='uniform', lower=0.0, upper=180.)
orbit['ecosw'] = 0.05
ecosw = orbit.get_parameter('ecosw')
ecosw.set_adjust(True)
ecosw.set_prior(distribution='uniform', lower=-1, upper=1)
orbit['esinw'] = 0.05
esinw = orbit.get_parameter('esinw')
esinw.set_adjust(True)
esinw.set_prior(distribution='uniform', lower=-1, upper=1)
fitparams = phoebe.ParameterSet('fitting:lmfit', method='nelder')
star1 = phoebe.BinaryRocheStar(comp1, mesh=mesh, pbdep=[lcdep,rvdep])#, orbit=orbit)
star2 = phoebe.BinaryRocheStar(comp2, mesh=mesh, pbdep=[lcdep,rvdep])#, orbit=orbit)
system = phoebe.BinaryBag([star1, star2], orbit=orbit, obs=[input_lc])
params = phoebe.ParameterSet('compute', eclipse_alg='binary')
feedback = phoebe.fitting.run(system, params, fitparams=fitparams, mpi=None, accept=True)
system.save('mymodel.phoebe')
feedback.save('myfeedback.phoebe')
system = phoebe.load_body('mymodel.phoebe')
feedback = phoebe.load_ps('myfeedback.phoebe')
syn = system.get_synthetic(category='lc', ref='my kepler lc').asarray()
plt.figure()
plt.plot(syn['time'], syn['flux'], 'ko-')
syn1 = system[0].get_synthetic(category='rv').asarray()
syn2 = system[1].get_synthetic(category='rv').asarray()
plt.figure()
plt.plot(syn1['time'], syn1['rv'], 'ko-')
plt.plot(syn2['time'], syn2['rv'], 'ro-')
plt.show()