#!/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_star()

times, wavelengths, fluxes, sigmas = np.loadtxt("Sed2.dat", usecols=[0, 1, 2, 3], unpack=True)

b.add_dataset('sed', times=times, wavelengths=wavelengths, fluxes=fluxes, sigmas=sigmas)

print("b['sed@sed01@sed@dataset'] = ", b['sed@sed01@sed@dataset'])
print("b['times@sed01@sed@dataset'] = ", b['times@sed01@sed@dataset'])
print("b['wavelengths@sed01@sed@dataset'] = ", b['wavelengths@sed01@sed@dataset'])
print("b['fluxes@sed01@sed@dataset'] = ", b['fluxes@sed01@sed@dataset'])

b.set_value('distance', context='system', value=1.*units.au)
b.set_value('teff@starA', context='component', value=30000.)
b.set_value('ntriangles@starA', context='compute', value=100)

b.run_compute()

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

# Note: Model (synthetic) baselines and wavelengths must be copied ex-post!
#b.set_value('u1@latest@model', value=b.get_value('u1@dataset'), ignore_readonly=True)
#b.set_value('v1@latest@model', value=b.get_value('v1@dataset'), ignore_readonly=True)
#b.set_value('u2@latest@model', value=b.get_value('u2@dataset'), ignore_readonly=True)
#b.set_value('v2@latest@model', value=b.get_value('v2@dataset'), ignore_readonly=True)
#b.set_value('wavelengths@latest@model', value=b.get_value('wavelengths@dataset'), ignore_readonly=True)

print("")
print("b['sed@sed01@phoebe01@latest@sed@model'] = ", b['sed@sed01@phoebe01@latest@sed@model'])
#print("")
#print("b['times@sed01@phoebe01@latest@sed@model'] = ", b['times@sed01@phoebe01@latest@sed@model'])
#print("b['wavelengths@sed01@phoebe01@latest@sed@model'] = ", b['wavelengths@sed01@phoebe01@latest@sed@model'])
#print("b['fluxes@sed01@phoebe01@latest@sed@model'] = ", b['fluxes@sed01@phoebe01@latest@sed@model'])

times = b['times@sed01@phoebe01@latest@sed@model'].value
wavelengths = b['wavelengths@sed01@phoebe01@latest@sed@model'].value
fluxes = b['fluxes@sed01@phoebe01@latest@sed@model'].value

np.savetxt('test_phoebe13.out', np.c_[times, wavelengths, fluxes], header='times wavelenghts fluxes')

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