import math,json,csv
from pathlib import Path
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
out=Path(__file__).resolve().parent/'orbital-study';out.mkdir(exist_ok=True)
mu=398600.435507;R=6371;g0=9.80665
orbit=[{'altitude_km':h,'speed_km_s':math.sqrt(mu/(R+h)),'period_min':2*math.pi*math.sqrt((R+h)**3/mu)/60,'gravity_m_s2':1000*mu/(R+h)**2} for h in range(200,2001,10)]
ring=[]
for target in [.3,.5,1]:
 for rpm in [.5,1,2,3,4]:
  w=2*math.pi*rpm/60;r=target*g0/w**2;ring.append({'floor_g0':target,'rpm':rpm,'radius_m':r,'rim_speed_m_s':w*r,'head_foot_decrease_percent':200/r,'coriolis_at_radial_1m_s_m_s2':2*w})
for name,rows in [('orbits',orbit),('rotation',ring)]:
 with (out/f'{name}.csv').open('w',newline='') as f:writer=csv.DictWriter(f,fieldnames=rows[0]);writer.writeheader();writer.writerows(rows)
(out/'study.json').write_text(json.dumps({'model':'Ideal spherical two-body circular orbit and rigid steady rotation; independent models.','mu_km3_s2':mu,'earth_radius_km':R,'standard_gravity_m_s2':g0,'orbits':orbit,'rotation':ring,'limitations':['No drag or orbital lifetime prediction','No structural or biological validation','No uncertainty bands: deterministic formula evaluations']},indent=2))
plt.rcParams.update({'font.size':10,'axes.spines.top':False,'axes.spines.right':False,'svg.fonttype':'none','pdf.fonttype':42})
fig,ax=plt.subplots(1,2,figsize=(10,4),layout='constrained')
for a,key,label in zip(ax,['speed_km_s','period_min'],['Circular speed (km/s)','Orbital period (min)']):
 a.plot([r['altitude_km'] for r in orbit],[r[key] for r in orbit],color='#146b8c');a.set(xlabel='Altitude above model Earth (km)',ylabel=label);a.grid(alpha=.2)
 base=next(r for r in orbit if r['altitude_km']==400);a.plot(400,base[key],'o',color='#8943a2');a.annotate('400 km', (400,base[key]),xytext=(15,8),textcoords='offset points')
for ext in ['svg','png','pdf']:fig.savefig(out/f'orbits.{ext}',dpi=200,bbox_inches='tight')
plt.close(fig);fig,ax=plt.subplots(1,2,figsize=(10,4),layout='constrained')
for target,col in zip([.3,.5,1],['#146b8c','#8943a2','#c76d31']):
 x=[.5+i*.01 for i in range(351)];rad=[target*g0/(2*math.pi*rpm/60)**2 for rpm in x]
 for a,y in zip(ax,[rad,[200/r for r in rad]]):a.plot(x,y,color=col,label=f'{target:g} g₀')
 rr=[r for r in ring if r['floor_g0']==target]
 ax[0].scatter([r['rpm'] for r in rr],[r['radius_m'] for r in rr],color=col,s=18);ax[1].scatter([r['rpm'] for r in rr],[r['head_foot_decrease_percent'] for r in rr],color=col,s=18)
ax[0].set_yscale('log');ax[0].set_ylabel('Required radius (m; logarithmic scale)');ax[1].set_ylabel('2 m head-to-foot decrease (% of floor value)')
for a in ax:a.set_xlabel('Rotation rate (rpm)');a.grid(alpha=.2);a.legend(frameon=False,title='Floor acceleration')
for ext in ['svg','png','pdf']:fig.savefig(out/f'rotation.{ext}',dpi=200,bbox_inches='tight')
for row in orbit:
 r=R+row['altitude_km'];assert abs(row['speed_km_s']*row['period_min']*60-2*math.pi*r)<1e-8
assert orbit[0]['speed_km_s']>orbit[-1]['speed_km_s'] and orbit[0]['period_min']<orbit[-1]['period_min']
print('181 orbital rows, 15 rotation rows, six figures generated; circumference invariant passed.')