"""Own scientific SVGs from computed Roche volumes and an attributed native aggregate."""
from pathlib import Path
import hashlib
import json
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from roche_reference import l1, omega, volume_radius
from envelope_reference import ce_final_separation

import argparse
parser=argparse.ArgumentParser()
parser.add_argument('--results',type=Path,default=Path('results'))
arguments=parser.parse_args()
OUTPUT = arguments.results/'figures'
OUTPUT.mkdir(parents=True, exist_ok=True)
plt.rcParams.update({'font.size': 11, 'svg.fonttype': 'path', 'svg.hashsalt': 'ac29-owned-figures', 'axes.spines.top': False, 'axes.spines.right': False})


def save(figure, name):
    figure.savefig(OUTPUT/name, bbox_inches='tight', metadata={'Date': None})
    figure.savefig(OUTPUT/(Path(name).stem+'.png'), bbox_inches='tight', dpi=140)
    plt.close(figure)


def roche_figure():
    angle = np.linspace(0, 2*np.pi, 721)
    mu = .5; upper = np.full_like(angle, l1(1)+mu); lower = np.zeros_like(angle)
    for _ in range(54):
        middle = (lower+upper)/2
        interior = omega(mu, -mu+middle*np.cos(angle), middle*np.sin(angle)) > omega(mu, 0, 0)
        lower = np.where(interior, middle, lower); upper = np.where(interior, upper, middle)
    radius = (lower+upper)/2; x = -.5+radius*np.cos(angle); y = radius*np.sin(angle)
    volume = volume_radius(1, 256)
    figure, axis = plt.subplots(figsize=(7.2, 4.3))
    axis.fill(x,y,color='#d3a060',alpha=.25); axis.plot(x,y,color='#ae6b19',lw=2)
    axis.fill(-x,y,color='#78b4d4',alpha=.22); axis.plot(-x,y,color='#286d99',lw=2)
    axis.plot(-.5+volume*np.cos(angle),volume*np.sin(angle),'--',color='#667689',lw=1.2)
    axis.scatter([-.5,.5,0],[0,0,0],c=['#ae6b19','#286d99','#4b5c69'],s=[32,32,18])
    axis.text(-.57,-.07,r'$M_d$'); axis.text(.46,-.07,r'$M_a$'); axis.text(.02,.035,r'$L_1$')
    axis.annotate('',xy=(0,.47),xytext=(-.5,.47),arrowprops={'arrowstyle':'<->','color':'#4b5c69'})
    axis.text(-.36,.50,r'$d_{L_1}=0.500a$'); axis.text(-.91,-.47,rf'$R_L={volume:.4f}a$')
    axis.set(xlabel=r'$x/a$',ylabel=r'$y/a$',xlim=(-1.03,1.03),ylim=(-.55,.58)); axis.set_aspect('equal'); axis.grid(alpha=.15)
    save(figure,'roche-volume.svg')


def orbit_figure():
    removed = np.linspace(0,.08,201)
    figure, axis = plt.subplots(figsize=(7.2,4.3))
    for q,color in [(.4,'#ae6b19'),(1.5,'#286d99'),(3,'#147866')]:
        separation=(q/((q-q*removed)*(1+q*removed)))**2
        axis.plot(removed,separation,color=color,lw=2,label=rf'$q_0={q:g}$')
    axis.set(xlabel=r'$x/M_{d0}$',ylabel=r'$a/a_0$'); axis.grid(alpha=.18); axis.legend(frameon=False)
    save(figure,'conservative-orbits.svg')


def envelope_figure():
    profile=json.loads((arguments.results/'structure-aggregate.json').read_text())[0]
    md,mc,ma=profile['massSolar'],profile['coreMassSolar'],1.4
    ai=profile['radiusSolar']/volume_radius(md/ma,256)
    alpha=np.linspace(.1,1,201)
    af=np.array([ce_final_separation(md,mc,ma,ai,profile['gravitationalCostUnits'],value) for value in alpha])
    contact=profile['coreCutRadiusSolar']+2e-5
    core_fit=profile['coreCutRadiusSolar']/volume_radius(mc/ma,256)
    figure,axis=plt.subplots(figsize=(7.2,4.3))
    axis.plot(alpha,af,color='#ae6b19',lw=2,label=r'$a_f$')
    axis.axhline(contact,color='#a34b49',linestyle=':',label=r'$R_c+R_a$')
    axis.axhline(core_fit,color='#667689',linestyle='--',label=r'$R_c/f(q_c)$')
    axis.fill_between(alpha,0,contact,color='#e9c1bf',alpha=.3)
    axis.set(xlabel=r'$\alpha_{\rm CE}$',ylabel=r'$a_f/R_\odot$',ylim=(0,.48)); axis.grid(alpha=.18); axis.legend(frameon=False)
    save(figure,'conditional-envelope.svg')


roche_figure(); orbit_figure(); envelope_figure()
rows=[]
for path in sorted(OUTPUT.glob('*.svg')):
    body=path.read_bytes(); rows.append({'path':str(path),'bytes':len(body),'sha256':hashlib.sha256(body).hexdigest()})
(arguments.results/'reader-figures.json').write_text(json.dumps({'kind':'Own computed scientific figures; no paper imagery copied', 'figures':rows,
    'nativeAttribution':'Envelope figure is our gravity/core-cut aggregate derived from Thoul & Chidester2022 Zenodo7082536, CC BY4.0; original SHA6200a8540d6d1da9e4b90c2bccb9c99c0ee609d5e3b352c8d6c3c5bf84ce204f. Not a new MESA run.'},indent=2)+'\n')
print('PASS',len(rows),'own scientific SVG figures')