#!/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_spectroscopy14(verbose=True):

#    phoebe.conf.devel_on()

    b = phoebe.default_star()

    # Note: Using a tiny test grid.
    from phoebe.backend import spectroscopy
    from phoebe.backend import pyterpolmu

    spectroscopy.sg = pyterpolmu.SyntheticGrid(gridlist='gridlist_ostar', wmin=6500., wmax=6600.)
    #spectroscopy.sg = pyterpolmu.SyntheticGrid(gridlist='gridlist_ostar')  # NOT WORKING! cf. non-equidistant wavelengts

    b.flip_constraint('logg@starA', solve_for='requiv@starA')

    b.set_value('distance@system', value=100*u_.pc)
    b.set_value('mass@starA@component', value=10.0*u_.solMass)
    b.set_value('teff@starA@component', value=30000*u_.K)
    b.set_value('period@starA@component', value=1000.0*u_.d)
    b.set_value('ntriangles@starA@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', 'intensities@lc01'])

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

    f = open('twigs.txt', 'w')
    for twig in b.twigs:
      f.write("%s\n" % (twig))
    f.close()
 
    for val in [3.5, 3.0, 2.9, 2.8, 2.7]:

        b.set_value('logg@starA@component', value=val)

        print(b['logg@starA@component'])
        print(b['requiv@starA@component'])

        b.run_compute()

        fig, plt = b['mesh'].plot(show=True, fc='intensities', ec=None, draw_sidebars=True, fcmap='gray', fclim=[0.0, 1.0e-13])
        plt.savefig("mesh_%0.2f.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_spectroscopy14_%.2f.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_spectroscopy14()