#!/usr/bin/env python3

import sys
import numpy as np
import phoebe
from astropy import units

#logger = phoebe.logger(clevel='INFO')
#logger = phoebe.logger(clevel='DEBUG')

b = phoebe.default_binary()

times, u, v, wavelengths, vises, sigmas = np.loadtxt("Vis.dat", usecols=[0, 1, 2, 3, 5, 6], unpack=True)

b.add_dataset('vis', times=times, u=u, v=v, wavelengths=wavelengths, vises=vises, sigmas=sigmas, if_method='integrate')

#print("b['times@vis01@vis@dataset'] = ", b['times@vis01@vis@dataset'])
#print("b['u@vis01@vis@dataset'] = ", b['u@vis01@vis@dataset'])
#print("b['v@vis01@vis@dataset'] = ", b['v@vis01@vis@dataset'])
#print("b['wavelengths@vis01@vis@dataset'] = ", b['wavelengths@vis01@vis@dataset'])
#print("b['vises@vis01@vis@dataset'] = ", b['vises@vis01@vis@dataset'])
#print("b['sigmas@vis01@vis@dataset'] = ", b['sigmas@vis01@vis@dataset'])
#print("")
#sys.exit(1)

b.flip_constraint('mass@primary@constraint', solve_for='period@binary')

#print("b['constraint'] = ", b['constraint'])

b.set_value('distance', context = 'system', value=100*units.pc)
b.set_value('sma@binary@component', value=5.3*units.solRad)

#print("b['sma@binary@component'] =", b['sma@binary@component'])
#print("b['period@binary@component'] =", b['period@binary@component'])
#print("b['t0_perpass@binary@component'] =", b['t0_perpass@binary@component'])

#b.add_dataset('mesh', compute_phases=[0.25])
#b.run_compute()
#b.plot(kind='mesh', show=True)

print("b['fit_parameters'] = ", b['fit_parameters'])

b.add_solver('optimizer.nelder_mead', solver='nm_solver')

print("b['fit_parameters'] = ", b['fit_parameters'])

b.set_value('maxiter', 20)
#b.set_value('maxiter', 1)
b.set_value('fit_parameters', ['sma@binary@component'])
b.get_value('fit_parameters', expand=True)

b.run_solver('nm_solver', solution='sol')

print(b.filter(solution='sol'))

b.adopt_solution('sol')
b.run_compute()

f = open('twigs.txt', 'w')
for twig in b.twigs:
  f.write("%s\n" % (twig))
f.close()

times = b['times@vis01@phoebe01@latest@vis@model'].value
u = b['u@vis01@phoebe01@latest@vis@model'].value
v = b['v@vis01@phoebe01@latest@vis@model'].value
wavelengths = b['wavelengths@vis01@phoebe01@latest@vis@model'].value
vises = b['vises@vis01@phoebe01@latest@vis@model'].value

np.savetxt('test_phoebe6.out', np.c_[times, u, v, wavelengths, vises], header='times u v wavelenghts vises')

#b.plot(show=True)
b.plot(x='u', marker='.', linestyle='none', show=True)
#b.plot(x='v', marker='.', linestyle='none', show=True)