#!/usr/bin/env python3

import sys
import numpy as np
import phoebe

b = phoebe.default_binary()

times, epochs, sigmas = np.loadtxt("Etv.dat", usecols=[0, 1, 2], unpack=True)

t0 = 0.0
period = 1.0
#t0 = 0.5
#period = 1.5

print("times = ", times)
print("epochs = ", epochs)
print("sigmas = ", sigmas)

b.add_dataset('etv', times=times, epochs=epochs, sigmas=sigmas, t0=t0, period=period)

print("b['times@etv01@etv@dataset'] = ", b['times@etv01@etv@dataset'])
print("b['epochs@etv01@etv@dataset'] = ", b['epochs@etv01@etv@dataset'])
print("b['sigmas@etv01@etv@dataset'] = ", b['sigmas@etv01@etv@dataset'])
print("b['time_ecls@etv01@etv@dataset'] = ", b['time_ecls@etv01@etv@dataset'])
print("b['time_ephs@etv01@etv@dataset'] = ", b['time_ephs@etv01@etv@dataset'])
print("b['etvs@etv01@etv@dataset'] = ", b['etvs@etv01@etv@dataset'])

#print("b['t0@etv01@etv@dataset'] = ", b['t0@etv01@etv@dataset'])
#print("b['period@etv01@etv@dataset'] = ", b['period@etv01@etv@dataset'])
#print("b['dpdt@etv01@etv@dataset'] = ", b['dpdt@etv01@etv@dataset'])

#b.plot(x='times', y='sigmas', show=True)
#sys.exit(1)

b.run_compute(etv_tol=0.1)

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

print("b['times@etv01@phoebe01@latest@etv@model'] = ", b['times@etv01@phoebe01@latest@etv@model'])
print("b['time_ecls@etv01@phoebe01@latest@etv@model'] = ", b['time_ecls@etv01@phoebe01@latest@etv@model'])
print("b['time_ephs@etv01@phoebe01@latest@etv@model'] = ", b['time_ephs@etv01@phoebe01@latest@etv@model'])
print("b['etvs@etv01@phoebe01@latest@etv@model'] = ", b['etvs@etv01@phoebe01@latest@etv@model'])
print("")

times = b['times@etv01@phoebe01@latest@etv@model'].value
time_ecls = b['time_ecls@etv01@phoebe01@latest@etv@model'].value
time_ephs = b['time_ephs@etv01@phoebe01@latest@etv@model'].value
etvs = b['etvs@etv01@phoebe01@latest@etv@model'].value

np.savetxt('test_phoebe25.out', np.c_[times, time_ecls, time_ephs, etvs], header='times times_ecls times_ephs etvs')

b.plot(show=True)
#b.plot(x='times', y='fluxes', show=True)
#b.plot(yerror='sigmas', show=True)  # IT DOES NOT WORK, WHY?

sys.exit(1)