#!/usr/bin/env python3

"""
Problem: The c. o. m. is NOT in the centre of (u, v) coordinates!

Problem: Body 3 is computed in a Roche potential of the c. o. m. (1+2),
not of individual bodies! For larger distances, it should not matter.

Problem: Bodies 1, 2 are computed in a Roche potential of bodies 1, 2;
unaffected by body 3! For larger distances, it should not matter.

"""

import sys
import numpy as np
import phoebe
from astropy import units as u

from phoebe.parameters import hierarchy as _hierarchy

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

phoebe.conf.devel_on()

b = phoebe.Bundle()

b.add_star(component='starA', color='red')
b.add_star(component='starB', color='green')
b.add_star(component='starC', color='blue')

b.add_orbit(component='orbit1', period=1.0)
b.add_orbit(component='orbit2', period=1.0)  # cf. non-physical value for testing!!!

hier1 = _hierarchy.binaryorbit(b['orbit1'], b['starA'], b['starB'])
hier2 = _hierarchy.binaryorbit(b['orbit2'], hier1, b['starC'])

b.set_hierarchy(hier2)

print("hier = ", b.get_hierarchy())

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

print("DONE!")

b.add_constraint(phoebe.parameters.constraint.keplers_third_law_hierarchical, 'orbit2', 'orbit1')
b.add_compute(dynamics_method='keplerian')

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

b.set_value('requiv@starA', value=2.0*u.solRad)
b.set_value('requiv@starB', value=2.0*u.solRad)
b.set_value('requiv@starC', value=2.0*u.solRad)
b.set_value('ntriangles@starA', context='compute', value=500)
b.set_value('ntriangles@starB', context='compute', value=500)
b.set_value('ntriangles@starC', context='compute', value=500)

b.run_compute()

#system = b['system']
#meshes = system.meshes
#mesh.get_columns_flat('')

fig, plt = b.plot(show=True)
plt.savefig('test_phoebe18_compact.png')

#print("b['00.000000@xyz_elements@starA@mesh01@phoebe01@latest@mesh@model'] = ", b['00.000000@xyz_elements@starA@mesh01@phoebe01@latest@mesh@model'])
#print("b['00.000000@xyz_elements@starB@mesh01@phoebe01@latest@mesh@model'] = ", b['00.000000@xyz_elements@starB@mesh01@phoebe01@latest@mesh@model'])
#print("b['00.000000@xyz_elements@starC@mesh01@phoebe01@latest@mesh@model'] = ", b['00.000000@xyz_elements@starC@mesh01@phoebe01@latest@mesh@model'])

sys.exit(1)