"""Synthetic known-correspondence attitude and relative same-band brightness."""
import math,csv,json
from pathlib import Path
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
P=Path(__file__).resolve().parent
SEED=20260927
uv=np.array([[-.2,-.12],[-.12,.15],[0,-.18],[.08,.1],[.18,-.04],[.22,.17]])
b=np.column_stack([uv,np.ones(6)]);b/=np.linalg.norm(b,axis=1)[:,None]
x,y,z=np.deg2rad([3,-5,7]);Rx=np.array([[1,0,0],[0,np.cos(x),-np.sin(x)],[0,np.sin(x),np.cos(x)]]);Ry=np.array([[np.cos(y),0,np.sin(y)],[0,1,0],[-np.sin(y),0,np.cos(y)]]);Rz=np.array([[np.cos(z),-np.sin(z),0],[np.sin(z),np.cos(z),0],[0,0,1]]);R=Rz@Ry@Rx
s=b@R

def fit(src,dst):
 if len(src)<2 or np.linalg.matrix_rank(src,tol=1e-10)<2:return None
 U,_,Vt=np.linalg.svd(dst.T@src);return U@np.diag([1,1,np.linalg.det(U@Vt)])@Vt

def angle(r):
 # atan2 avoids loss of precision near identity; independent acos check below.
 return math.atan2(np.linalg.norm([r[2,1]-r[1,2],r[0,2]-r[2,0],r[1,0]-r[0,1]])/2,(np.trace(r)-1)/2)
rng=np.random.default_rng(SEED);trials=[]
for noise in [0,.1,1]:
 for trial in range(100):
  pixels=uv*800+[640,480]+rng.normal(0,noise,(6,2));obs=np.column_stack([(pixels-[640,480])/800,np.ones(6)]);obs/=np.linalg.norm(obs,axis=1)[:,None];est=fit(s,obs);err=angle(est@R.T)
  trials.append(dict(centroid_noise_sd_px=noise,trial=trial+1,rotation_error_rad=err,rotation_error_deg=math.degrees(err),fit_vector_RMS=float(np.sqrt(np.mean(np.sum((s@est.T-obs)**2,axis=1))))))
vec=[dict(direction_id=i+1,inertial_x=s[i,0],inertial_y=s[i,1],inertial_z=s[i,2],camera_x=b[i,0],camera_y=b[i,1],camera_z=b[i,2],pixel_x=uv[i,0]*800+640,pixel_y=uv[i,1]*800+480) for i in range(6)]
def magnitude(q):
 if not math.isfinite(q) or q<=0:return None
 return -2.5*math.log10(q)
assert all(magnitude(q) is None for q in [0,-1,math.nan,math.inf])
assert magnitude(1)==0 and magnitude(100)==-5 and magnitude(.01)==5
# Independent trace-based rotation angle at a finite test rotation.
assert abs(angle(R)-math.acos(float(np.clip((np.trace(R)-1)/2,-1,1))))<1e-14
mag=[dict(flux_ratio=q,relative_magnitude=-2.5*math.log10(q),extinction_mag=A,attenuation_factor=10**(-.4*A),observed_flux_ratio=q*10**(-.4*A),observed_relative_magnitude=-2.5*math.log10(q)+A) for q in np.logspace(-2,2,101) for A in [0,1,2]]
summ=[]
for noise in [0,.1,1]:
 e=[r['rotation_error_deg'] for r in trials if r['centroid_noise_sd_px']==noise];summ.append(dict(centroid_noise_sd_px=noise,trials=100,mean_rotation_error_deg=float(np.mean(e)),p95_rotation_error_deg=float(np.quantile(e,.95))))
for name,rows in [('directions.csv',vec),('attitude-trials.csv',trials),('attitude-summary.csv',summ),('magnitudes.csv',mag)]:
 with (P/name).open('w',newline='') as f:w=csv.DictWriter(f,fieldnames=rows[0]);w.writeheader();w.writerows(rows)
assert np.max(abs(s@s.T-b@b.T))<1e-14
assert np.linalg.norm(fit(s,b)-R)<1e-14
assert np.linalg.norm(fit(s,s)-np.eye(3))<1e-14
perm=[4,1,5,0,3,2];assert np.linalg.norm(fit(s[perm],b[perm])-R)<1e-14
assert fit(s[:1],b[:1]) is None and fit(np.repeat(s[:1],2,axis=0),np.repeat(b[:1],2,axis=0)) is None
assert angle(fit(s,b[perm])@R.T)>.01
assert all(abs(10**(-.4*r['relative_magnitude'])/r['flux_ratio']-1)<1e-14 for r in mag)
assert all(abs(-2.5*math.log10(r['observed_flux_ratio'])-r['observed_relative_magnitude'])<1e-14 for r in mag)
assert all(0<=v['pixel_x']<1280 and 0<=v['pixel_y']<960 for v in vec)
assert max(r['rotation_error_deg'] for r in trials if r['centroid_noise_sd_px']==0)<1e-10
plt.rcParams.update({'font.family':'DejaVu Sans','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.5,4.7));ax[0].scatter([v['pixel_x'] for v in vec],[v['pixel_y'] for v in vec],color='#087f92',s=45)
for v in vec:ax[0].annotate(str(v['direction_id']),(v['pixel_x']+8,v['pixel_y']+8),fontsize=9)
ax[0].set(xlabel='Image x (pixels)',ylabel='Image y (pixels, upward convention)',xlim=(0,1280),ylim=(0,960),title='A  |  Six synthetic, known directions');ax[0].set_aspect('equal')
for i,(noise,color) in enumerate([(.1,'#087f92'),(1,'#a44786')]):
 vals=sorted(r['rotation_error_deg'] for r in trials if r['centroid_noise_sd_px']==noise);ax[1].plot(vals,np.arange(1,101)/100,label=f'{noise:g} pixel SD',color=color)
ax[1].set(xlabel='Rotation error (degrees)',ylabel='Empirical cumulative fraction of 100 trials',title='B  |  Conditional on correct correspondences');ax[1].legend(fontsize=8)
for a in ax:a.grid(alpha=.2)
fig.suptitle('Orientation from directions: a synthetic calibration exercise',fontsize=14);fig.text(.5,.04,'800-pixel focal length; fixed six-direction geometry; independent Gaussian centroid errors. Zero-noise recovery is at numerical precision.',ha='center',fontsize=8);fig.text(.5,.012,'No real stars, catalogue matching, false detections or hardware accuracy. The fitted rotation maps inertial directions into camera coordinates.',ha='center',fontsize=8);fig.tight_layout(rect=[0,.09,1,.93])
for ext in ['png','svg','pdf']:fig.savefig(P/f'attitude-study.{ext}',dpi=190)
fig,ax=plt.subplots(figsize=(8,4.5))
for A,color in [(0,'#087f92'),(1,'#a44786'),(2,'#7550a0')]:
 rr=[r for r in mag if r['extinction_mag']==A];ax.semilogx([r['flux_ratio'] for r in rr],[r['observed_relative_magnitude'] for r in rr],label=f'Assigned extinction {A} mag',color=color)
ax.set(xlabel='Unattenuated flux / reference flux (same band)',ylabel='Observed magnitude − reference magnitude',title='Greater magnitude means less detected flux');ax.legend(fontsize=9);ax.grid(alpha=.2);fig.text(.5,.025,'Relative same-band scale only: a factor 100 brighter corresponds to −5 magnitudes. Extinction adds to magnitude.',ha='center',fontsize=8);fig.tight_layout(rect=[0,.065,1,1])
for ext in ['png','svg','pdf']:fig.savefig(P/f'relative-magnitude.{ext}',dpi=190)
meta={'title':'Synthetic known-correspondence attitude and relative magnitude','status':'Seeded numerical illustration; no star recognition or hardware test','seed':SEED,'rows':{'directions':6,'trials':300,'summary':3,'magnitudes':303},'rotation_inertial_to_camera':R.tolist(),'parameters':{'focal_length_px':800,'image_size_px':[1280,960],'principal_point_px':[640,480],'pixel_y_convention':'upward','centroid_noise_sd_px':[0,.1,1],'trials_each':100,'extinction_mag':[0,1,2]},'fit':'Equal-weight Wahba orthogonal Procrustes: H=sum(b*s^T), U,S,Vt=svd(H), R=U diag(1,1,det(UVt)) Vt','checks':{'all_pass':True,'unit_vector_pair_dot_invariance':True,'zero_noise_rotation_recovery':True,'identity_rotation':True,'joint_permutation_invariance':True,'wrong_correspondence_detectably_wrong':True,'one_vector_and_collinear_no_attitude':True,'magnitude_inverse_and_extinction':True},'limits':['Six synthetic directions have known correct labels. No catalogue search, source detection, recognition success rate or geographic position is estimated.','Gaussian pixel centroid errors are independent; focal length, distortion and principal point are exact. Conditional error distributions are not camera certification.','One direction or collinear repeated directions cannot determine full attitude; near-parallel geometry can also be ill-conditioned.','Magnitude uses positive flux in the same photometric band/system; zero or negative flux has no ordinary logarithmic magnitude.','Extinction is assigned per band. No distance, intrinsic luminosity, universal photometric zero point or signal-to-noise magnitude uncertainty is inferred.']}
(P/'study.json').write_text(json.dumps(meta,indent=2)+'\n');(P/'validation.json').write_text(json.dumps(meta['checks'],indent=2)+'\n');print(summ)