#!/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, u1, v1, u2, v2, wavelengths, clos, sigmas = np.loadtxt("Clo.dat", usecols=[0, 1, 2, 3, 4, 5, 7, 8], unpack=True)

b.add_dataset('clo', times=times, u1=u1, v1=v1, u2=u2, v2=v2, wavelengths=wavelengths, clos=clos, sigmas=sigmas, if_method='integrate', passband='Johnson:V')

#b.set_value_all('atm', 'blackbody')
#b.set_value_all('ld_mode', 'manual')
#b.set_value_all('ld_func', 'linear')
#b.set_value_all('ld_coeffs', [0.0])

print("b['times@clo01@clo@dataset'] = ", b['times@clo01@clo@dataset'])
print("b['u1@clo01@clo@dataset'] = ", b['u1@clo01@clo@dataset'])
print("b['v1@clo01@clo@dataset'] = ", b['v1@clo01@clo@dataset'])
print("b['u2@clo01@clo@dataset'] = ", b['u2@clo01@clo@dataset'])
print("b['v2@clo01@clo@dataset'] = ", b['v2@clo01@clo@dataset'])
print("b['wavelengths@clo01@clo@dataset'] = ", b['wavelengths@clo01@clo@dataset'])
print("b['clos@clo01@clo@dataset'] = ", b['clos@clo01@clo@dataset'])

b.set_value('distance', context = 'system', value=100*units.pc)
b.set_value('incl@binary@orbit', value=80*units.deg)
b.set_value('teff@secondary@component', value=5000.0)

b.run_compute()

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

print("")
print("b['times@clo01@phoebe01@latest@clo@model'] = ", b['times@clo01@phoebe01@latest@clo@model'])
print("b['u1@clo01@phoebe01@latest@clo@model'] = ", b['u1@clo01@phoebe01@latest@clo@model'])
print("b['v1@clo01@phoebe01@latest@clo@model'] = ", b['v1@clo01@phoebe01@latest@clo@model'])
print("b['u2@clo01@phoebe01@latest@clo@model'] = ", b['u2@clo01@phoebe01@latest@clo@model'])
print("b['v2@clo01@phoebe01@latest@clo@model'] = ", b['v2@clo01@phoebe01@latest@clo@model'])
print("b['wavelengths@clo01@phoebe01@latest@clo@model'] = ", b['wavelengths@clo01@phoebe01@latest@clo@model'])
print("b['clos@clo01@phoebe01@latest@clo@model'] = ", b['clos@clo01@phoebe01@latest@clo@model'])

times = b['times@clo01@phoebe01@latest@clo@model'].value
u1 = b['u1@clo01@phoebe01@latest@clo@model'].value
v1 = b['v1@clo01@phoebe01@latest@clo@model'].value
u2 = b['u2@clo01@phoebe01@latest@clo@model'].value
v2 = b['v2@clo01@phoebe01@latest@clo@model'].value
wavelengths = b['wavelengths@clo01@phoebe01@latest@clo@model'].value
clos = b['clos@clo01@phoebe01@latest@clo@model'].value

np.savetxt('test_phoebe5.out', np.c_[times, u1, v1, u2, v2, wavelengths, clos], header='times u1 v1 u2 v2 wavelenghts clos')

b.plot(show=True)
b.plot(x='u1', show=True)
b.plot(x='u2', show=True)