#!/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, t3s, sigmas = np.loadtxt("T3.dat", usecols=[0, 1, 2, 3, 4, 5, 7, 8], unpack=True)

b.add_dataset('t3', times=times, u1=u1, v1=v1, u2=u2, v2=v2, wavelengths=wavelengths, t3s=t3s, sigmas=sigmas, if_method='simple', passband='Johnson:V')

print("b['times@t301@t3@dataset'] = ", b['times@t301@t3@dataset'])
print("b['u1@t301@t3@dataset'] = ", b['u1@t301@t3@dataset'])
print("b['v1@t301@t3@dataset'] = ", b['v1@t301@t3@dataset'])
print("b['u2@t301@t3@dataset'] = ", b['u2@t301@t3@dataset'])
print("b['v2@t301@t3@dataset'] = ", b['v2@t301@t3@dataset'])
print("b['wavelengths@t301@t3@dataset'] = ", b['wavelengths@t301@t3@dataset'])
print("b['t3s@t301@t3@dataset'] = ", b['t3s@t301@t3@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@t301@phoebe01@latest@t3@model'] = ", b['times@t301@phoebe01@latest@t3@model'])
print("b['u1@t301@phoebe01@latest@t3@model'] = ", b['u1@t301@phoebe01@latest@t3@model'])
print("b['v1@t301@phoebe01@latest@t3@model'] = ", b['v1@t301@phoebe01@latest@t3@model'])
print("b['u2@t301@phoebe01@latest@t3@model'] = ", b['u2@t301@phoebe01@latest@t3@model'])
print("b['v2@t301@phoebe01@latest@t3@model'] = ", b['v2@t301@phoebe01@latest@t3@model'])
print("b['wavelengths@t301@phoebe01@latest@t3@model'] = ", b['wavelengths@t301@phoebe01@latest@t3@model'])
print("b['t3s@t301@phoebe01@latest@t3@model'] = ", b['t3s@t301@phoebe01@latest@t3@model'])

times = b['times@t301@phoebe01@latest@t3@model'].value
u1 = b['u1@t301@phoebe01@latest@t3@model'].value
v1 = b['v1@t301@phoebe01@latest@t3@model'].value
u2 = b['u2@t301@phoebe01@latest@t3@model'].value
v2 = b['v2@t301@phoebe01@latest@t3@model'].value
wavelengths = b['wavelengths@t301@phoebe01@latest@t3@model'].value
t3s = b['t3s@t301@phoebe01@latest@t3@model'].value

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

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