#!/usr/bin/env python3

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

from phoebe.parameters import hierarchy as _hierarchy

def test_spectroscopy6(verbose=True):

#    phoebe.conf.devel_on()

    b = phoebe.default_star()

    b.set_value('distance', context='system', value=100*u_.pc)
    b.set_value('requiv@starA@component', value=1.0*u_.solRad)
    b.set_value('ntriangles@starA', context='compute', value=2000)

    dir_ = os.path.dirname(os.path.realpath(__file__))

    times, wavelengths, fluxes, sigmas = np.loadtxt(os.path.join(dir_, "Spe.dat"), usecols=[0, 1, 2, 3], unpack=True)

    b.add_dataset('spe', times=times, wavelengths=wavelengths, fluxes=fluxes, sigmas=sigmas, passband='Johnson:R', use_instrumental=False, use_rotational=True)

    b.add_dataset('mesh', times=[0.25])
    b.add_dataset('lc', times=[0.25])

    b.set_value('columns', value=['abs_intensities@lc01'])

    f = open('twigs.txt', 'w')
    for twig in b.twigs:
      f.write("%s\n" % (twig))
    f.close()
 
    for val in [0.16, 0.17, 0.20, 0.30, 0.50, 99.00]:

        b.set_value('period@starA@component', value=val*u_.d)

        b.run_compute()

        fig, plt = b['mesh'].plot(show=True, fc='abs_intensities', ec=None, draw_sidebars=True, fcmap='gray', fclim=[0.0, 5.0e13])
        plt.savefig("mesh_%0.3f.png" % val)

        times = b['times@spe01@phoebe01@latest@spe@model'].value
        wavelengths = b['wavelengths@spe01@phoebe01@latest@spe@model'].value
        fluxes = b['fluxes@spe01@phoebe01@latest@spe@model'].value

        np.savetxt('test_spectroscopy6_%.3f.out' % val, np.c_[times, wavelengths, fluxes], header='times wavelenghts fluxes')


if __name__ == "__main__":

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

    test_spectroscopy6()