"""Synthetic reflectance, areal mixing and normalized top-hat band averages."""
import csv,json,math
from pathlib import Path
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
P=Path(__file__).resolve().parent
BANDS=[(400,500),(500,600),(600,750),(750,900),(900,1000)]
def spectra(w):
 a=.55-.12*np.exp(-.5*((w-700)/70)**2);return a,a-.2*np.exp(-.5*((w-850)/8)**2)
def band(alpha,lo,hi,scale=1,tilt=0,step=.1):
 if not all(math.isfinite(v) for v in [alpha,lo,hi,scale,tilt,step]) or not 0<=alpha<=1 or not 400<=lo<hi<=1000 or scale<0 or abs(tilt)>.8 or step<=0:raise ValueError('Outside model domain')
 if scale==0:return None
 w=np.linspace(lo,hi,round((hi-lo)/step)+1);a,b=spectra(w);e=scale*(1+tilt*(w-700)/300)
 return float(np.trapezoid(e*(alpha*a+(1-alpha)*b),w)/np.trapezoid(e,w))
spec=[]
for w in range(400,1001):
 a,b=spectra(w);spec.append(dict(wavelength_nm=w,A_reflectance=float(a),B_reflectance=float(b),equal_areal_mix_reflectance=float((a+b)/2)))
rows=[]
for alpha in np.linspace(0,1,21):
 for tilt in [-.8,0,.8]:
  for lo,hi in BANDS:rows.append(dict(A_area_fraction=float(alpha),illumination_tilt=tilt,band_low_nm=lo,band_high_nm=hi,normalized_band_reflectance=band(float(alpha),lo,hi,tilt=tilt)))
widths=[dict(band_width_nm=w,A_band=band(1,850-w/2,850+w/2),B_band=band(0,850-w/2,850+w/2)) for w in range(10,201,2)]
for r in widths:r['A_minus_B']=r['A_band']-r['B_band']
for name,data in [('spectra.csv',spec),('bands.csv',rows),('width-sensitivity.csv',widths)]:
 with (P/name).open('w',newline='') as f:w=csv.DictWriter(f,fieldnames=data[0]);w.writeheader();w.writerows(data)
assert all(0<=r[k]<=1 for r in spec for k in ['A_reflectance','B_reflectance','equal_areal_mix_reflectance'])
assert all(0<=r['normalized_band_reflectance']<=1 for r in rows)
err=0;conv=0;analytic=0
for r in rows:
 a=r['A_area_fraction'];lo=r['band_low_nm'];hi=r['band_high_nm'];tilt=r['illumination_tilt'];v=r['normalized_band_reflectance']
 err=max(err,abs(v-(a*band(1,lo,hi,tilt=tilt)+(1-a)*band(0,lo,hi,tilt=tilt))))
 conv=max(conv,abs(v-band(a,lo,hi,tilt=tilt,step=.05)))
 assert abs(v-band(a,lo,hi,scale=7,tilt=tilt))<1e-14
 assert band(a,lo,hi,scale=0,tilt=tilt) is None
# Independent closed-form Gaussian integral for flat illumination.
def gaussian_integral(lo,hi,mu,sd):return sd*math.sqrt(math.pi/2)*(math.erf((hi-mu)/(math.sqrt(2)*sd))-math.erf((lo-mu)/(math.sqrt(2)*sd)))
for lo,hi in BANDS:
 a=.55-.12*gaussian_integral(lo,hi,700,70)/(hi-lo);b=a-.2*gaussian_integral(lo,hi,850,8)/(hi-lo)
 analytic=max(analytic,abs(a-band(1,lo,hi)),abs(b-band(0,lo,hi)))
width_error=max(abs(r['A_minus_B']-.2*gaussian_integral(850-r['band_width_nm']/2,850+r['band_width_nm']/2,850,8)/r['band_width_nm']) for r in widths)
assert width_error<3e-6
assert err<1e-14 and conv<1e-7 and analytic<1e-7
assert all(a['A_minus_B']>b['A_minus_B'] for a,b in zip(widths,widths[1:]))
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=(11,4.7));x=[r['wavelength_nm'] for r in spec]
for k,label,color in [('A_reflectance','Synthetic A','#087e91'),('B_reflectance','Synthetic B','#a14786'),('equal_areal_mix_reflectance','Equal areal mixture','#777')]:ax[0].plot(x,[r[k] for r in spec],label=label,color=color)
for lo,hi in BANDS:ax[0].axvspan(lo,hi,alpha=.04 if lo%3 else .09,color='#664488');ax[0].axvline(lo,color='#bbb',lw=.5)
ax[0].set(xlabel='Wavelength (nm)',ylabel='Reflectance (dimensionless)',title='A  |  A narrow feature differs');ax[0].legend(fontsize=8)
for alpha,label,color in [(1,'A','#087e91'),(0,'B','#a14786'),(.5,'Equal mixture','#777')]:ax[1].plot(range(5),[band(alpha,lo,hi) for lo,hi in BANDS],'-o',label=label,color=color)
ax[1].set(xticks=range(5),xticklabels=[f'{lo}–{hi}' for lo,hi in BANDS],xlabel='Top-hat wavelength band (nm)',ylabel='Normalized band reflectance',title='B  |  Broad bands average the feature');ax[1].legend(fontsize=8)
for a in ax:a.grid(alpha=.2);a.set_ylim(.2,.6)
fig.suptitle('Synthetic spectra: a band is an average, not a full spectrum',fontsize=14);fig.text(.5,.015,'Analytic synthetic reflectances, flat illumination and ideal top-hat response. No named material, measured spectrum or classifier is represented.',ha='center',fontsize=8);fig.tight_layout(rect=[0,.055,1,.93])
for ext in ['png','svg','pdf']:fig.savefig(P/f'spectra-and-bands.{ext}',dpi=190)
fig,ax=plt.subplots(1,2,figsize=(10.5,4.6))
ax[0].plot([r['band_width_nm'] for r in widths],[r['A_minus_B'] for r in widths],color='#75469c');ax[0].set(xlabel='Band width centred at 850 nm (nm)',ylabel='A − B band-reflectance difference',title='A  |  Wider bands dilute a narrow difference')
for tilt,color in [(-.8,'#087e91'),(0,'#777'),(.8,'#a14786')]:
 rr=[r for r in rows if r['illumination_tilt']==tilt and r['band_low_nm']==750];ax[1].plot([r['A_area_fraction'] for r in rr],[r['normalized_band_reflectance'] for r in rr],color=color,label=f'Illumination tilt {tilt:+.1f}')
ax[1].set(xlabel='Synthetic A area fraction',ylabel='Normalized 750–900 nm band reflectance',title='B  |  Mixture and weighting both matter');ax[1].legend(fontsize=8)
for a in ax:a.grid(alpha=.2)
fig.suptitle('Information depends on the band and the viewing assumptions',fontsize=14);fig.text(.5,.015,'Common geometry and illumination; unresolved adjacent patches, not intimate powder mixing. Zero illumination gives undefined normalized reflectance.',ha='center',fontsize=8);fig.tight_layout(rect=[0,.055,1,.93])
for ext in ['png','svg','pdf']:fig.savefig(P/f'band-sensitivity.{ext}',dpi=190)
meta={'title':'Synthetic spectral bands and areal mixtures','status':'Analytical hypothetical spectra; no real material or classification result','rows':{'spectra':len(spec),'bands':len(rows),'width_sensitivity':len(widths)},'equations':{'A':'.55-.12*exp(-.5*((lambda_nm-700)/70)^2)','B':'A-.20*exp(-.5*((lambda_nm-850)/8)^2)','mixture':'alpha*A+(1-alpha)*B','illumination':'scale*(1+tilt*(lambda_nm-700)/300)','band':'integral(E*R*d_lambda)/integral(E*d_lambda), top-hat response'},'parameters':{'wavelength_nm':[400,1000],'bands_nm':BANDS,'A_fractions':[i/20 for i in range(21)],'tilts':[-.8,0,.8],'integration_step_nm':.1},'checks':{'all_pass':True,'mix_linearity_max_abs_error':err,'half_step_max_abs_difference':conv,'closed_form_flat_band_max_abs_error':analytic,'width_sensitivity_closed_form_max_abs_error':width_error,'bounds':True,'illumination_scale_cancellation':True,'dark_denominator_null':True,'width_dilution_monotone':True},'limitations':['Normalized band response presumes a same-illumination reference. It is not raw detector counts or a direct uncalibrated reflectance observation.','Zero illumination produces no recoverable normalized reflectance, represented by null, not a zero-reflectance material.','Linear areal mixing assumes unresolved adjacent patches under common viewing and lighting; it excludes intimate mixture scattering, shadows and nonlinear interactions.','No atmospheric absorption, emissivity, thermal emission, spatial misregistration, detector noise or calibration error is modeled.','No named mineral, classifier accuracy, abundance estimate or confidence probability follows.']}
(P/'study.json').write_text(json.dumps(meta,indent=2)+'\n');(P/'validation.json').write_text(json.dumps(meta['checks'],indent=2)+'\n');print(meta['checks'])