"""Independent reproducible calculations for the bounded Gaia/2MASS study.
Python3 standard library only. Run after acquire.py, or against packaged CSVs.
Uniform rectilinear geometric barycentric model; see methods.json.
"""
from pathlib import Path
import csv,json,math,hashlib,statistics
BASE=Path(__file__).resolve().parent
MASRAD=math.pi/(180*3_600_000)
PC_KM=149597870.7*648000/math.pi
YEAR_SECONDS=365.25*86400
CORR=[(0,1,'ra_dec_corr'),(0,2,'ra_parallax_corr'),(0,3,'ra_pmra_corr'),(0,4,'ra_pmdec_corr'),(1,2,'dec_parallax_corr'),(1,3,'dec_pmra_corr'),(1,4,'dec_pmdec_corr'),(2,3,'parallax_pmra_corr'),(2,4,'parallax_pmdec_corr'),(3,4,'pmra_pmdec_corr')]
def num(s):return None if s is None or s=='' or s.lower() in ('null','nan') else float(s)
def wrap(x):return (x+180)%360-180
def position(g,year,delta=None):
 if not 1990<=year<=2040:raise ValueError('Epoch outside bounded1990–2040 model')
 d=delta or [0]*6
 a=math.radians(g['ra'])+d[0]*MASRAD/math.cos(math.radians(g['dec']));b=math.radians(g['dec'])+d[1]*MASRAD
 plx=g['parallax']+d[2]
 if plx<=0:raise ValueError('Positive parallax required')
 dist=1000/plx;rv=(g['radial_velocity'] if g['radial_velocity'] is not None else 0)+d[5]
 ca,sa,cb,sb=math.cos(a),math.sin(a),math.cos(b),math.sin(b)
 n=[cb*ca,cb*sa,sb];p=[-sa,ca,0];q=[-sb*ca,-sb*sa,cb]
 vr=rv*YEAR_SECONDS/PC_KM
 v=[dist*MASRAD*((g['pmra']+d[3])*p[i]+(g['pmdec']+d[4])*q[i])+vr*n[i] for i in range(3)]
 vec=[dist*n[i]+v[i]*(year-g['ref_epoch']) for i in range(3)];norm=math.hypot(*vec)
 return {'raDeg':math.degrees(math.atan2(vec[1],vec[0]))%360,'decDeg':math.degrees(math.atan2(vec[2],math.hypot(vec[0],vec[1]))),'distancePc':norm}
def covariance(g,year):
 base=position(g,year);errors=[g[x] for x in ['ra_error','dec_error','parallax_error','pmra_error','pmdec_error']]+[g['radial_velocity_error'] or 0]
 c=[[errors[i]**2 if i==j else 0 for j in range(6)] for i in range(6)]
 for i,j,key in CORR:c[i][j]=c[j][i]=g[key]*errors[i]*errors[j]
 steps=[.1,.1,max(.001,g['parallax']*1e-5),.01,.01,.05];jac=[[],[]]
 for k,h in enumerate(steps):
  delta=[0]*6;delta[k]=h;plus=position(g,year,delta);delta[k]=-h;minus=position(g,year,delta)
  jac[0].append(wrap(plus['raDeg']-minus['raDeg'])*math.cos(math.radians(base['decDeg']))*3_600_000/(2*h));jac[1].append((plus['decDeg']-minus['decDeg'])*3_600_000/(2*h))
 out=[[sum(jac[i][k]*c[k][l]*jac[j][l] for k in range(6) for l in range(6)) for j in range(2)] for i in range(2)]
 xx,xy,yy=out[0][0],(out[0][1]+out[1][0])/2,out[1][1];disc=math.hypot(xx-yy,2*xy)
 return {'xxMas2':xx,'xyMas2':xy,'yyMas2':yy,'majorMas':math.sqrt(max(0,(xx+yy+disc)/2)),'minorMas':math.sqrt(max(0,(xx+yy-disc)/2)),'angleRad':.5*math.atan2(2*xy,xx-yy)}
def separation(a,b):
 ar,ad,br,bd=map(math.radians,[a['raDeg'],a['decDeg'],b['raDeg'],b['decDeg']]);da=br-ar
 y=math.hypot(math.cos(bd)*math.sin(da),math.cos(ad)*math.sin(bd)-math.sin(ad)*math.cos(bd)*math.cos(da));x=math.sin(ad)*math.sin(bd)+math.cos(ad)*math.cos(bd)*math.cos(da)
 return math.degrees(math.atan2(y,x))*3600
def analyse_record(record):
 g=record['gaia'];before=[];after=[]
 for t in record['tmass']:
  observed={'raDeg':t['raDeg'],'decDeg':t['decDeg']};old=separation({'raDeg':g['ra'],'decDeg':g['dec']},observed);new=separation(position(g,t['epochJulianYear']),observed)
  before.append(old);after.append(new)
 return {'targetId':record['id'],'gaiaSourceId':g['source_id'],'candidates':len(before),'unpropagatedNearestArcsec':min(before) if before else None,'epochPropagatedNearestArcsec':min(after) if after else None,'withinOneArcsecBefore':sum(x<=1 for x in before),'withinOneArcsecAfter':sum(x<=1 for x in after)}
def build():
 gaia=list(csv.DictReader((BASE/'gaia-dr3-source.csv').open()));tmass=list(csv.DictReader((BASE/'tmass-neighbours.csv').open()));records=[]
 for row in gaia:
  g={k:(v if k=='source_id' else num(v)) for k,v in row.items()};assert len(g['source_id'])>=15 and g['ref_epoch']==2016
  assert all(g[k] is not None for k in ['ra','dec','parallax','pmra','pmdec','ra_error','dec_error','parallax_error','pmra_error','pmdec_error']+[x[2] for x in CORR])
  ts=[]
  for r in tmass:
   if r['source_id']!=g['source_id']:continue
   t={'designation':r['designation'],'raDeg':num(r['tm_ra']),'decDeg':num(r['tm_dec']),'julianDate':num(r['j_date']),'epochJulianYear':2000+(num(r['j_date'])-2451545)/365.25,'errorMajorArcsec':num(r['err_maj']),'errorMinorArcsec':num(r['err_min']),'errorAngleDeg':num(r['err_ang']),'jMag':num(r['j_m']),'hMag':num(r['h_m']),'ksMag':num(r['ks_m']),'photometryQuality':r['ph_qual'],'officialCandidateAngularDistanceArcsec':num(r['angular_distance']),'officialCandidateScore':num(r['score']),'officialCandidateFlag':int(r['xm_flag'])}
   if t['designation'] not in {v['designation'] for v in ts}:ts.append(t)
  records.append({'id':'gaia-'+g['source_id'],'name':'Gaia DR3 '+g['source_id'],'gaia':g,'tmass':ts})
 assert len(records)==20
 sources=[json.loads((BASE/f).read_text()) for f in ['gaia-source-provenance.json','tmass-provenance.json']]
 for source in sources:assert source['sha256']==hashlib.sha256((BASE/source['file']).read_bytes()).hexdigest()
 methods={'id':'uniform-rectilinear-barycentric-v1','epochRange':[1990,2040],'epochUnit':'Julian year; Gaia reference2016.0(TCB)','positionFrame':'ICRS barycentric geometric direction','radialVelocityPolicy':'Use reported Gaia radial velocity when available; otherwise conditional zero radial velocity. Missing RV uncertainty is unquantified, never measured zero.','covariance':'Published Gaia5-parameter covariance with finite-difference2x6Jacobian; RV variance independent when reported. Principal-axis1sigma lengths are not a68percent2Dcoverage region.','excluded':['Annual observer parallax','Solar-system aberration','Light-time/relativistic corrections','Galactic gravitational acceleration','Binary orbital motion','Extinction and Gaia parallax zero-point correction'],'tmassEpoch':'Approximate Julian year from reported2MASS JD; UTC/TCB seconds and JD±30s uncertainty are neglected at this educational precision.','crossmatch':'Only the packaged published Gaia2MASS good-neighbourhood candidate list, not a complete cone search. Radius selection is not proof of physical identity.'}
 packet={'schema':'org.spacetrc.stellar-archive','schemaVersion':1,'edition':'9.0.0','sources':sources,'selectionProvenance':json.loads((BASE/'selection-provenance.json').read_text()),'methods':methods,'credits':{'Gaia':'ESA/Gaia/DPAC. Gaia Collaboration etal2023; originalDR3credit instructions linked inpage.','2MASS':'2MASS:UniversityofMassachusetts/IPAC-Caltech; NASA/NSF.Skrutskie etal2006.'},'records':records}
 (BASE/'data.json').write_text(json.dumps(packet,ensure_ascii=False,separators=(',',':'))+'\n')
 targets={'schemaVersion':1,'targets':[{'id':r['id'],'name':r['name'],'raDeg':r['gaia']['ra'],'decDeg':r['gaia']['dec'],'kind':'star','epoch':2016.0,'gaiaSourceId':r['gaia']['source_id'],'sourceUrl':sources[0]['url']} for r in records]};(BASE/'targets.json').write_text(json.dumps(targets,ensure_ascii=False,indent=2)+'\n')
 results=[analyse_record(r) for r in records];usable=[r for r in results if r['candidates']]
 summary={'sources':20,'sourcesWithPublishedCandidates':len(usable),'candidateRows':sum(r['candidates'] for r in results),'sourcesWithRV':sum(r['gaia']['radial_velocity'] is not None for r in records),'medianUnpropagatedArcsec':statistics.median(r['unpropagatedNearestArcsec'] for r in usable),'medianPropagatedArcsec':statistics.median(r['epochPropagatedNearestArcsec'] for r in usable),'sourcesWithCandidateWithinOneArcsecBefore':sum(r['withinOneArcsecBefore']>0 for r in results),'sourcesWithCandidateWithinOneArcsecAfter':sum(r['withinOneArcsecAfter']>0 for r in results)}
 study={'schema':'org.spacetrc.stellar-repeat-study','schemaVersion':1,'edition':'9.0.0','status':'archived-data-worked-repeat; no independent specialist review or newdiscovery','question':{'tr':'Aynı yıldız adayını farklı kataloglarda ararken gözlem dönemi neden önemlidir?','en':'Why does observation epoch matter when comparing candidate records across catalogues?'},'selection':'20publishedGaiaDR3sources:parallax_over_error>20,ruwe<1.4,Gbetween8and13,properMotion>300mas/year,declinationbetween−60and60; first20byfixedrandom_index. Existingpublished2MASSneighbours arepreselected; this is not an unbiasedcompleteness/precisionbenchmark.','methods':methods,'summary':summary,'results':results,'dataSha256':hashlib.sha256((BASE/'data.json').read_bytes()).hexdigest()}
 (BASE/'study.json').write_text(json.dumps(study,ensure_ascii=False,indent=2)+'\n')
 with (BASE/'study.csv').open('w',newline='') as f:
  w=csv.DictWriter(f,fieldnames=list(results[0]));w.writeheader();w.writerows(results)
 fixtures=[]
 for r in records:
  for year in [1999.0,2016.0,2026.0,2040.0]:fixtures.append({'targetId':r['id'],'epoch':year,'position':position(r['gaia'],year),'covariance':covariance(r['gaia'],year)})
 (BASE/'reference-results.json').write_text(json.dumps(fixtures,separators=(',',':'))+'\n')
 print(json.dumps(summary,indent=2))
if __name__=='__main__':build()