"""Private numerical source qualification, not a mission background prediction.

HEALPix centre conversion adapted from Healpix_cxx3.83 pix2loc/nest2xyf,
Copyright (C)2003-2024 Max-Planck-Society, author Martin Reinecke.
This script is GPL-2.0-or-later. No warranty. Reference URLs retained below.
"""
from pathlib import Path
from datetime import datetime,timezone
import numpy as np,json,gzip,hashlib
BASE=Path('D:/Photometer_research_20261022/research/backgrounds')
N=BASE/'round09_metadata';S=BASE/'round10_starlight'
H=6.62607015e-34;C=299792458.;RP=1e10/(4*np.pi)
def bandfold(x,y,centre,width,gaussian=False):
 lo=centre-(5*width if gaussian else width/2);hi=centre+(5*width if gaussian else width/2)
 grid=np.unique(np.r_[lo,x[(x>lo)&(x<hi)],np.linspace(lo,hi,int(np.ceil((hi-lo)/.005))+1) if gaussian else [],hi]);v=np.interp(grid,x,y)
 filt=np.exp(-4*np.log(2)*((grid-centre)/width)**2)if gaussian else np.ones_like(grid)
 return float(np.trapezoid(v*filt,grid)),float(np.trapezoid(filt,grid))
def nest_vectors(nside):
 order=int(np.log2(nside));assert 2**order==nside
 p=np.arange(12*nside*nside,dtype=np.int64);face=p>>(2*order);local=p&(nside*nside-1)
 ix=np.zeros_like(p);iy=np.zeros_like(p)
 for k in range(order):ix|=((local>>(2*k))&1)<<k;iy|=((local>>(2*k+1))&1)<<k
 jrll=np.array([2]*4+[3]*4+[4]*4);jpll=np.array([1,3,5,7,0,2,4,6,1,3,5,7])
 jr=jrll[face]*nside-ix-iy-1;nr=np.where(jr<nside,jr,np.where(jr>3*nside,4*nside-jr,nside))
 z=np.where(jr<nside,1-nr**2/(3.*nside**2),np.where(jr>3*nside,nr**2/(3.*nside**2)-1,(2*nside-jr)*2/(3.*nside)))
 t=jpll[face]*nr+ix-iy;t=np.where(t<0,t+8*nr,t);phi=np.pi*t/(4*nr)
 v=np.column_stack((np.sqrt(1-z*z)*np.cos(phi),np.sqrt(1-z*z)*np.sin(phi),z))
 assert np.max(abs(np.linalg.norm(v,axis=1)-1))<1e-14
 if nside==1:
  assert np.allclose(np.mod(phi,2*np.pi),np.pi/180*np.array([45,135,225,315,0,90,180,270,45,135,225,315]))
 return v

rows=np.genfromtxt(N/'ngc_contcomp_n4.dat',names=True,skip_header=20)
x=rows['lam'];cont=rows['ftot'];uv=(x>=325)&(x<=350)
np.savetxt(Path('D:\\Photometer_research_20261022\\revision_20261022_r4\\analysis\\analytic\\background_components\\Noll2024_325_350nm_mean_continuum.csv'),np.column_stack([x[uv],cont[uv]]),delimiter=',',header='wavelength_nm_pipeline_air_inferred,mean_continuum_R_per_nm',comments='')
noll=[]
for width in [1.,2.,4.,8.]:
 for gaussian in [False,True]:
  r,eq=bandfold(x,cont,337.1,width,gaussian)
  # Air wavelengths in pipeline inferred from official ESO; energy at nominal vacuum
  # is a <0.03% convention sensitivity here, not a high-resolution registration claim.
  e,_=bandfold(x,cont*RP*H*C/(x*1e-9),337.1,width,gaussian)
  noll.append({'centre_in_tabulated_nm':337.1,'width_nm':width,'filter':'unit_peak_Gaussian_FWHM'if gaussian else'unit_peak_tophat_full_width','equivalent_width_nm':eq,'integrated_R':r,'band_mean_R_per_nm':r/eq,'radiance_W_m2_sr':e,'mean_W_m2_sr_nm':e/eq,'photon_radiance_m2_sr_s':r*RP})

j=np.loadtxt(S/'j_lambda.dat');lam=j[:,0]*1000
with gzip.open(S/'i_lambda.dat.gz','rt')as f:sky=np.loadtxt(f)
assert sky.shape==(49152,99) and j.shape==(99,5) and np.isfinite(sky).all() and (sky>=0).all()
mean_check=np.max(abs(sky.mean(axis=0)-j[:,1])/j[:,1]);assert mean_check<5e-5
v=nest_vectors(64);domega=4*np.pi/len(v)
check={ 'sky_cube_mean_vs_tabulated_J_max_relative_difference':float(mean_check),'vector_mean':v.mean(axis=0).tolist(),'mean_z_squared':float(np.mean(v[:,2]**2)),'nside1_z':nest_vectors(1)[:,2].tolist()}
normals=np.vstack((np.eye(3),-np.eye(3)))
check['unit_radiance_hemisphere_integrals_sr']=[float(np.maximum(v@n,0).sum()*domega)for n in normals]
assert np.max(abs(v.mean(axis=0)))<1e-14 and abs(np.mean(v[:,2]**2)-1/3)<2e-5
assert max(abs(np.array(check['unit_radiance_hemisphere_integrals_sr'])-np.pi))<0.0003
star=[]
for centre in [337.1,777.4]:
 k=int(np.searchsorted(lam,centre));frac=(centre-lam[k-1])/(lam[k]-lam[k-1]);radiance=(sky[:,k-1]*(1-frac)+sky[:,k]*frac)*1e-9
 jval=np.array([np.interp(centre,lam,j[:,col])*1e-9 for col in [1,2,3,4]])
 hemi=[float(np.sum(radiance*np.maximum(v@n,0))*domega)for n in normals]
 star.append({'requested_nm':centre,'bracketing_tabulated_nm':lam[k-1:k+1].tolist(),'linear_interpolation_weight':float(frac),'J_total_W_m2_sr_nm':float(jval[0]),'tabulated_total_uncertainty_interpolated_W_m2_sr_nm':float(jval[1]),'J_direct_W_m2_sr_nm':float(jval[2]),'J_diffuse_W_m2_sr_nm':float(jval[0]-jval[2]),'J_loglog_interpolation_total_W_m2_sr_nm':float(np.exp(np.interp(np.log(centre),np.log(lam),np.log(j[:,1])))*1e-9),'hypothetical_cloud_normal_unitvectors_celestial_xyz':normals.tolist(),'no_atmosphere_downwelling_irradiance_W_m2_nm':hemi,'Lambert_albedo0p8_no_atmosphere_reflected_radiance_W_m2_sr_nm':[.8*y/np.pi for y in hemi]})
result={'status':'PRIVATE_QUALIFICATION_NOT_MISSION_PREDICTION','recorded_utc':datetime.now(timezone.utc).isoformat(),'noll_primary':'https://doi.org/10.5194/acp-24-1143-2024','noll_dataset_version':'https://doi.org/10.5281/zenodo.8335837','noll_limits':['Mean10633selected2009-2019CerroParanal spectra, >=10min exposure; Moonabovehorizonwithillumination>50%excluded. Not a polar or microsecond distribution.','Continuum extracted by median/percentile filters; strong resolved lines removed, unresolvedO2UVbands retained partly. Extraterrestrial components modeled/subtracted; extinction and vanRhijn90km corrected. Not total night sky or total777line background.','Header only says wavelengthnm. Official ESO says Xshooter pipeline wavelength in air; authors could have subsequently converted, not explicitly verified here. Fold centre337.1uses table convention; photonenergyuses nominal wavelength with <0.03%air/vac sensitivity.','4component NMF is statistical; labelledcomponents must not automatically be treated as puremolecular emissions.','0.5nmgrid/broad extraction smoothing cannot validate sub-nm filter choices.'], 'noll_filter_folds':noll,'starlight_primary':'https://doi.org/10.1051/0004-6361/202451086','starlight_catalogue':'https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/691/A330','starlight_limits':['Totalmeans direct+interstellar-dust-scattered starlight; no zodiacal, moon, airglow. Opticalphotometry+stellar-atmosphere/exctinction extrapolation, not a narrow337measurement.','99logbins span~0.1-4.3um: around337bracket~13nm wide. Cannot infer actual4nm spectralfeatures or own covariance from tabulated broad-model errors.','Hypotheticalsix cloudnormaldirections, noactualorbit/date/Sun/moon/cloudprofile. Not globalminimum/maximum or occurrenceweighteddistribution.','Earthfilled nadir does not see this external field directly; need hemisphericillumination and cloudBRDF/atmosphericscattering. Lambertalbedo.8is an explicit example only.','CubeisHEALPixNESTEDorder6in celestialRA/Dec despite paper figures plottedinGalacticcoordinates. Do not reinterpret cube inGalacticcoordinates.'], 'starlight_checks':check,'starlight_examples':star,'healpix_reference':['https://healpix.sourceforge.io/doc/html/Healpix_cxx/healpix__base_8cc_source.html#l00928','https://healpix.sourceforge.io/doc/html/Healpix_cxx/healpix__tables_8cc_source.html#l00068']}
Path('D:\\Photometer_research_20261022\\revision_20261022_r4\\analysis\\analytic\\background_components\\new_spectral_background_qualification.json').write_text(json.dumps(result,indent=2),encoding='utf8')
print(json.dumps({'noll4nm':[r for r in noll if r['width_nm']==4],'checks':check,'starlight':star},indent=2))