import math,json,csv,sys,hashlib
from pathlib import Path
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
P=Path(__file__).resolve().parent

def allocation(work,n,mode):
    if not 0<=work<=500 or n not in range(1,6) or mode not in ['off','idle']:raise ValueError('Invalid input')
    feasible=work<=100*n
    if not feasible:return dict(allocated_cpu_units=work,active_hosts=n,unused_mode=mode,feasible=False,allocated_cpu_load=None,power_W=None,facility_energy_J=None,energy_per_cpu_unit=None)
    u=work/(100*n);power=n*(70+180*u*u)+(5-n)*70*(mode=='idle');e=power*60*1.3
    return dict(allocated_cpu_units=work,active_hosts=n,unused_mode=mode,feasible=True,allocated_cpu_load=work,power_W=power,facility_energy_J=e,energy_per_cpu_unit=e/work if work else None)
energy=[allocation(w,n,mode)for w in range(0,501,5)for n in range(1,6)for mode in ['off','idle']]
thermal=[]
for n in range(1,6):
    power=70+180/n**2
    for t in range(0,1801,10):thermal.append(dict(active_hosts=n,time_s=t,active_host_power_W=power,temperature_C=25+.1*power*(-math.expm1(-t/300)),steady_temperature_C=25+.1*power))
MU=398600.435507;RE=6378.137

def orbit(h,w=550,d=.1):
    r=RE+h;rate=math.sqrt(MU/r**3);e=math.pi/6;psi=math.acos(RE/r*math.cos(e))-e;scale=1.22*w*1e-9*h*1000/d;duration=2*psi/rate
    return dict(altitude_km=h,circular_speed_km_s=math.sqrt(MU/r),period_min=2*math.pi/rate/60,near_nadir_Rayleigh_scale_m=scale,central_visibility_s=duration,feasible=scale<=5 and duration>=120)
orb=[orbit(h)for h in range(200,2001,100)]
sens=[]
for w in [450,550,650]:
    for d in [.05,.1,.2]:
        vals=[r['altitude_km']for r in [orbit(h,w,d)for h in range(200,2001,100)]if r['feasible']];sens.append(dict(wavelength_nm=w,aperture_m=d,feasible_altitudes_km=';'.join(map(str,vals)),feasible_count=len(vals)))
def frontier(rows):
    return [r for r in rows if not any(s[0]<=r[0]and s[1]>=r[1]and(s[0]<r[0]or s[1]>r[1])for s in rows)]
assert len(frontier([(r['near_nadir_Rayleigh_scale_m'],r['central_visibility_s'])for r in orb]))==19
assert frontier([(1,1),(1,1),(2,1),(1,0)])==[(1,1),(1,1)]
assert [r['altitude_km']for r in orb if r['feasible']]==[300,400,500,600,700]
maxgeo=0
for row in orb:
    r=RE+row['altitude_km'];e=math.pi/6
    # Ray-sphere intersection independently gives slant distance, then atan2 central angle.
    s=-RE*math.sin(e)+math.sqrt(r*r-RE*RE*math.cos(e)**2)
    psi=math.atan2(s*math.cos(e),RE+s*math.sin(e));duration=2*psi/math.sqrt(MU/r**3)
    maxgeo=max(maxgeo,abs(duration-row['central_visibility_s']))
assert maxgeo<1e-10
for row in energy:
    if row['feasible']:
        assert math.isclose(row['facility_energy_J']/60/1.3,row['power_W'],rel_tol=1e-14)
        if row['allocated_cpu_units']==0:assert row['energy_per_cpu_unit']is None
for n in range(1,6):
    seq=[r for r in thermal if r['active_hosts']==n];temp=25;p=seq[0]['active_host_power_W']
    for row in seq:
        if row['time_s']:temp=25+(temp-25)*math.exp(-10/300)+.1*p*(-math.expm1(-10/300))
        assert abs(temp-row['temperature_C'])<1e-12 and row['temperature_C']<=row['steady_temperature_C']
assert allocation(0,1,'idle')['power_W']==350 and allocation(0,1,'off')['power_W']==70
oldpath=P/'numeric-baseline.json';old=json.loads(oldpath.read_text())
keys={'circular_speed_km_s':'circular_speed_km_s','period_min':'period_min','near_nadir_Rayleigh_scale_m':'near_nadir_Rayleigh_scale_m','central_visibility_s':'central_visibility_s'}
parity=max(abs(a[k]-b[v])for a,b in zip(orb,old['rows'])for k,v in keys.items());assert parity<1e-10
for a,b in zip(sens,old['sensitivity']):assert a['feasible_count']==b['feasible_count'] and a['feasible_altitudes_km']==';'.join(map(str,b['feasible_altitudes_km']))
def save(name,rows):
    with(P/name).open('w',newline='',encoding='utf8')as f:w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)
for name,rows in [('energy.csv',energy),('temperature.csv',thermal),('orbital-trade.csv',orb),('optical-sensitivity.csv',sens)]:save(name,rows)
plt.rcParams.update({'font.size':10,'axes.spines.top':False,'axes.spines.right':False,'savefig.dpi':180})
fig,ax=plt.subplots(1,2,figsize=(11,4.5),layout='constrained')
for mode in ['off','idle']:
    rr=[r for r in energy if r['allocated_cpu_units']==100 and r['unused_mode']==mode];ax[0].plot([r['active_hosts']for r in rr],[r['facility_energy_J']/1000 for r in rr],marker='o',label='Unused hosts '+mode)
ax[0].set(xlabel='Hosts sharing an assigned 100 CPU-unit load',ylabel='Facility energy over 60 s (kJ)',xticks=range(1,6),title='A  Power assumption changes the comparison');ax[0].legend(fontsize=8)
for n in range(1,6):
    rr=[r for r in thermal if r['active_hosts']==n];ax[1].plot([r['time_s']for r in rr],[r['temperature_C']for r in rr],label=f'{n} active')
ax[1].set(xlabel='Time under constant load (s)',ylabel='One active host temperature (°C)',title='B  Exact hypothetical RC response');ax[1].legend(fontsize=8);fig.suptitle('Synthetic identical servers: same assigned CPU load; no completion or latency model')
for ext in ['png','svg','pdf']:fig.savefig(P/f'energy-and-temperature.{ext}')
plt.close(fig)
fig,ax=plt.subplots(figsize=(8,5),layout='constrained')
ax.plot([r['near_nadir_Rayleigh_scale_m']for r in orb],[r['central_visibility_s']for r in orb],color='gray')
for r in orb:
    ax.scatter(r['near_nadir_Rayleigh_scale_m'],r['central_visibility_s'],color='#278c89'if r['feasible']else'#7957a8');ax.annotate(str(r['altitude_km']),(r['near_nadir_Rayleigh_scale_m'],r['central_visibility_s']),xytext=(4,3),textcoords='offset points',fontsize=7)
ax.axvline(5,ls='--',color='gray');ax.axhline(120,ls='--',color='gray');ax.set(xlabel='Near-nadir Rayleigh separation scale (m; smaller preferred)',ylabel='Ideal central-pass visibility (s; larger preferred)',title='Illustrative orbit trade: all 19 points are non-dominated\nLabels: altitude in km; green: chosen 5 m / 120 s requirements')
for ext in ['png','svg','pdf']:fig.savefig(P/f'orbital-trade.{ext}')
plt.close(fig)
(P/'validation.json').write_text(json.dumps(dict(all_pass=True,rows=dict(energy=1010,thermal=905,orbital=19,sensitivity=9),feasible_energy_rows=sum(r['feasible']for r in energy),orbit_legacy_max_difference=parity,independent_triangle_max_difference_s=maxgeo,exact_rc_semigroup=True,strict_pareto_and_equal_points=True,zero_work_undefined=True),indent=2)+'\n')
(P/'study.json').write_text(json.dumps(dict(runtime=dict(python=sys.version,numpy=np.__version__,matplotlib=matplotlib.__version__),synthetic=True,server=dict(hosts=5,capacity_cpu_units=100,idle_W=70,dynamic_W=180,PUE=1.3,energy_duration_s=60,thermal_ambient_C=25,thermal_initial_C=25,Rth_K_W=.1,tau_s=300),orbit=dict(mu_km3_s2=MU,radius_km=RE,wavelength_nm=550,aperture_m=.1,elevation_deg=30),baseline_sha256=dict(orbital_results=hashlib.sha256(oldpath.read_bytes()).hexdigest(),generator=hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),original_results=old['source_sha256']),load_interpretation='Assigned CPU load is maintained over 60 s; no executed work, completion or latency model',scope='Server energy/RC is a new corrected model, not legacy scheduler reproduction; orbit numeric comparator reproduced with corrected optical terminology'),indent=2)+'\n')
print('PASS1943rows',parity,maxgeo)