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

    phoebe.conf.devel_on()

    b = phoebe.default_binary()

    b.flip_constraint('period@primary@constraint', solve_for='syncpar@primary@component')
    b.flip_constraint('mass@primary@constraint', solve_for='sma@binary')
    b.flip_constraint('mass@secondary@constraint', solve_for='q@binary')

    b.set_value('distance', context='system', value=10*u_.pc)
    b.set_value('period@binary', value=1000.0*u_.d)
    b.set_value('mass@primary@component', value=1.0*u_.solMass)
    b.set_value('requiv@primary@component', value=1.0*u_.solRad)
    b.set_value('ntriangles@primary@compute', value=2000)

    # https://phoebe-project.org/docs/2.4/examples/distortion_method_none
    b.set_value('mass@secondary@component', value=0.00001*u_.solMass)
    b.set_value('requiv@secondary@component', value=0.1*u_.solRad)
    b.set_value('ntriangles@secondary@compute', value=100)
    b.set_value('distortion_method@secondary@compute', value='none')

    if verbose:
        print("m1 = ", b['mass@primary@component'].value)
        print("m2 = ", b['mass@secondary@component'].value)
        print("a1 = ", b['sma@binary@component'].value)

    f = open('twigs.txt', 'w')
    for twig in b.twigs:
      f.write("%s\n" % (twig))
    f.close()
 
    dir_ = os.path.dirname(os.path.realpath(__file__))

    times, u, v, wavelengths, vises, sigmas = np.loadtxt(os.path.join(dir_, "Vis.dat"), usecols=[0, 1, 2, 3, 5, 6], unpack=True)

    b.add_dataset('vis', times=times, u=u, v=v, wavelengths=wavelengths, vises=vises, sigmas=sigmas, passband='Johnson:V')

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

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

    for val in [0.158, 0.16, 0.20, 0.30, 0.50, 99.00]:

        b.set_value('period@primary@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@vis01@phoebe01@latest@vis@model'].value
        u = b['u@vis01@phoebe01@latest@vis@model'].value
        v = b['v@vis01@phoebe01@latest@vis@model'].value
        wavelengths = b['wavelengths@vis01@phoebe01@latest@vis@model'].value
        vises = b['vises@vis01@phoebe01@latest@vis@model'].value

        np.savetxt('test_interferometry6_%.3f.out' % val, np.c_[times, u, v, wavelengths, vises], header='times u v wavelenghts vises')


if __name__ == "__main__":

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

    test_interferometry6()