#!/usr/bin/env python3
"""
Reproduce the design-mode north-aligned phase-core toy simulation.
This script intentionally assigns proxy_axis_azimuth_deg = 0 degrees true north.
It is a hypothesis/proxy simulation, not a diagnostic result.
"""
from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import json

outdir = Path('outputs/design_mode_north_aligned')
outdir.mkdir(parents=True, exist_ok=True)
G = 6.67430e-11; M_earth = 5.9722e24; R_earth = 6.371e6
M_sun = 1.98847e30; AU = 1.495978707e11; R_sun = 6.957e8
m_p = 1.67262192369e-27; k_B = 1.380649e-23
height_m = 65.0; g_parent = G*M_earth/(R_earth+height_m)**2
Gamma_lab=1.10; Gamma_wall=0.95; Gamma_b=1.00; Gamma_total=Gamma_lab*Gamma_wall*Gamma_b
A_N=1.0; A_opp_max=1.25
M_core0=0.08*M_sun; M_env=0.92*M_sun; R0=120*AU; R_stable=0.0030*AU
T_ignition=4.0e6; rho_ref=2e-14; alpha_conf=75000.0; beta_accrete=5.5
n_steps=420; tau=np.linspace(0,1,n_steps)
C_internal=1/(1+np.exp(-12*(tau-0.45))); C_internal=(C_internal-C_internal.min())/(C_internal.max()-C_internal.min())
A_opp=np.zeros(n_steps); A_opp[0]=A_opp_max; R_core=np.zeros(n_steps); R_core[0]=R0
records=[]; ignition_step=None; stabilization_step=None; coherence_step=None; prev_fusion_to_cooling=1e-9
for i in range(n_steps):
    C=C_internal[i]
    if i>0:
        A_required=1/(Gamma_total*max(C,1e-6)*A_N)
        A_target=min(A_opp_max,max(0.20,A_required))
        response_rate=0.06 if ignition_step is None else 0.16
        A_opp[i]=A_opp[i-1]+response_rate*(A_target-A_opp[i-1])
    R_phi=Gamma_total*C*A_N*A_opp[i]
    gnet_over_gparent=1-R_phi; g_net_m_s2=g_parent*gnet_over_gparent
    stability_index=float(np.clip(np.exp(-(gnet_over_gparent/0.09)**2)*np.clip(C,0,1)**0.35,0,1))
    compression_index=float(np.clip(C*np.exp(-(gnet_over_gparent/0.16)**2),0,1.2))
    M_core=M_core0+M_env*(1-np.exp(-beta_accrete*C**2))
    R_pre=R0/(1+alpha_conf*compression_index**3)
    fusion_pressure=0.0 if ignition_step is None else 1/(1+np.exp(-(np.log10(max(prev_fusion_to_cooling,1e-30))-0.5)*4))
    R_target=max(R_pre,R_stable*(1+0.03*fusion_pressure))
    if i>0: R_core[i]=R_core[i-1]+0.08*(R_target-R_core[i-1])
    T_core=max(30.0, G*M_core*0.61*m_p/(5*k_B*R_core[i])*(0.90+0.10*stability_index))
    rho_core=M_core/((4/3)*np.pi*R_core[i]**3); rho_rel=rho_core/rho_ref
    temp_gate=1/(1+np.exp(-(np.log10(max(T_core,1))-6.2)*8))
    fusion_to_cooling=(T_core/T_ignition)**5*(rho_rel/1e19)**0.25*temp_gate*(max(stability_index,1e-9)/0.60)**2
    if stability_index<0.05: fusion_to_cooling*=stability_index/0.05
    prev_fusion_to_cooling=fusion_to_cooling
    field_coherent=(abs(gnet_over_gparent)<=0.05) and (stability_index>=0.80)
    ignited=(T_core>=T_ignition) and (fusion_to_cooling>=1) and (stability_index>=0.55)
    radius_change_fraction=0 if i==0 else abs((R_core[i]-R_core[i-1])/max(R_core[i],1e-30))
    stabilized=ignited and abs(gnet_over_gparent)<=0.02 and stability_index>=0.85 and fusion_to_cooling>=1 and radius_change_fraction<=0.02
    if coherence_step is None and field_coherent: coherence_step=i
    if ignition_step is None and ignited: ignition_step=i
    if stabilization_step is None and ignition_step is not None and i>ignition_step+10 and stabilized: stabilization_step=i
    state='parent_field_dominated'
    if field_coherent and not ignited: state='field_coherent_pre_ignition'
    if ignited and not stabilized: state='ignited_seeking_stabilization'
    if stabilized: state='ignited_stabilized'
    records.append(dict(step=i,tau=tau[i],C_phi_internal=C,A_opp=A_opp[i],R_phi_design=R_phi,gnet_over_gparent=gnet_over_gparent,g_net_m_s2=g_net_m_s2,field_coherence_index=stability_index,M_core_Msun=M_core/M_sun,R_core_AU=R_core[i]/AU,R_core_Rsun=R_core[i]/R_sun,T_core_K=T_core,fusion_to_cooling=fusion_to_cooling,field_coherent=field_coherent,ignited=ignited,stabilized=stabilized,state=state))
df=pd.DataFrame(records); df.to_csv(outdir/'design_mode_north_aligned_coherence_timeseries.csv',index=False)
print({'coherence_step':coherence_step,'ignition_step':ignition_step,'stabilization_step':stabilization_step})
for y,label in [('R_phi_design','R_phi'),('T_core_K','T_core'),('fusion_to_cooling','fusion_to_cooling')]:
    plt.figure(figsize=(8,5)); plt.plot(df['tau'],df[y]); plt.title(label); plt.xlabel('tau'); plt.tight_layout(); plt.savefig(outdir/f'{label}.png',dpi=160)
