from pathlib import Path
import csv,datetime,json,math,sys
sys.path.insert(0,r'D:/Photometer_research_20261022/research/cloud_transport/private_readers')
from scipy.integrate import quad
from scipy.optimize import brentq

p=Path(r'D:/Photometer_research_20261022/revision_20261022_r4/analysis/analytic/reflection')
p.mkdir(exist_ok=True)
c=299792458.;area=math.pi*(.025/2)**2
def integrate_x(f,lo,hi,H,h):
 points=sorted({x for x in [h/H*.1,h/H,h/H*10,h/H*100,.001,.01,.1] if lo<x<hi})
 return quad(f,lo,hi,points=points,epsabs=1e-11,epsrel=1e-10,limit=200)[0]
def setup(top_km,h_m,albedo,fov):
 H=(550-top_km)*1000; h=h_m
 if h==0:
  return H,1.,2*albedo,lambda x:0.,lambda x:2*albedo
 xmax_radius=H*math.tan(math.radians(fov/2))
 xmin=h/math.hypot(h,xmax_radius)
 def radius(x):return h*math.sqrt(max(0.,1/x**2-1))
 def delay(x):
  r=radius(x)
  # Stablepositiveexpression avoids subtractingtwo~500km lengths.
  return (h/x+h+r*r/(math.hypot(H,r)+H))/c
 def weight(x):
  r=radius(x)
  return 2*albedo*H*H*(H-h)**2/(H*H+r*r)**2
 B=integrate_x(weight,xmin,1,H,h)
 return H,xmin,B,delay,weight

def rise(top,h,a,fov,input_rise_us):
 H,xmin,B,delay,weight=setup(top,h,a,fov)
 ramp=input_rise_us*1e-6/.8
 if h==0 or a==0:return B,input_rise_us
 dmin=delay(1);dmax=delay(xmin)
 def boundary(t):
  if t<=dmin:return 1.
  if t>=dmax:return xmin
  return brentq(lambda x:delay(x)-t,xmin,1.,xtol=1e-14)
 def received(t):
  direct=min(1.,max(0.,t/ramp))
  lo=boundary(t);hi=boundary(t-ramp)
  middle=integrate_x(lambda x:weight(x)*(t-delay(x))/ramp,lo,hi,H,h) if hi>lo else 0.
  flat=integrate_x(weight,hi,1.,H,h) if hi<1 else 0.
  return (direct+middle+flat)/(1+B)
 t10=brentq(lambda t:received(t)-.1,0,ramp+dmax,xtol=1e-12)
 t90=brentq(lambda t:received(t)-.9,0,ramp+dmax,xtol=1e-12)
 assert abs(received(t10)-.1)<1e-6 and abs(received(t90)-.9)<1e-6
 return B,(t90-t10)*1e6

rows=[]
for top,h_values in [(15,[0,1,100,1000,3000]),(17,[0,1,100,1000])]:
 for h in h_values:
  for albedo in [.8,1.]:
   for fov in [60,90]:
    H,xmin,B,delay,weight=setup(top,h,albedo,fov)
    incident_fraction_in_view=.5*(1-xmin) if h else .5
    # Independent radialcoordinate quadrature checks exactweight/Jacobian.
    if h:
     R=H*math.tan(math.radians(fov/2))
     ratio_radial=quad(lambda v:2*albedo*H*H*(H-h)**2*math.sinh(v)/(math.cosh(v)**2*(H*H+h*h*math.sinh(v)**2)**2),0,math.asinh(R/h),epsabs=1e-11,epsrel=1e-10,limit=300)[0]
     if abs(ratio_radial-B)>=2e-7: raise RuntimeError(f'coordinatecheckfailedtop={top},h={h},rho={albedo},fov={fov},B={B},radial={ratio_radial}')
    else:ratio_radial=B
    for tr in [30.,80.,200.]:
     boost_part,observed=rise(top,h,albedo,fov,tr)
     assert 0<=B<=2*albedo+1e-9
     rows.append(dict(plane_altitude_km=top,source_above_plane_m=h,source_altitude_km=top+h/1000,
      satellite_altitude_km=550,full_circular_FOV_deg=fov,assumed_Lambert_albedo=albedo,
      nadir_fluence_ratio_to_vacuum=1+B,reflected_to_direct_ratio=B,independent_radial_integral=ratio_radial,
      downward_source_fraction_intercepted_within_view=incident_fraction_in_view,
      all_directions_total_upward_escape_fraction=(1+albedo)/2,
      ideal_25mm_aperture_fraction_of_emitted_photons=area/(4*math.pi*(H-h)**2)*(1+B),
      minimum_reflected_extra_delay_us=0 if h==0 else 2*h/c*1e6,
      maximum_reflected_extra_delay_in_FOV_us=0 if h==0 else delay(xmin)*1e6,
      assumed_intrinsic_10_90_linear_rise_us=tr,observed_10_90_of_final_plateau_rise_us=observed,
      field_observer='nadir, horizontalflataperture; directpointsourceontheaxis',
      status='OWNEXACTSURFACETOYNOTCLOUDRT'))
with (p/'conditional_reflection_and_rise.csv').open('w',newline='',encoding='utf-8') as f:
 w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)
derivation={
 'utc':datetime.datetime.now(datetime.timezone.utc).isoformat(),
 'status':'PRIVATE_ANALYTIC_DIAGNOSTIC_NOT_A_REPLACEMENT_CLOUDSCAT_RUN',
 'question_tags':['A0','A2','G22'],
 'purpose':'Showwhydirectionalfluencecanexceedvacuumwithoutcreatingphotons, andwhetherreflectedgeometrycanbroadenanevenabove-topsourceunderoneexplicitidealmodel.',
 'primary_BRDF_reference':'https://www.pbr-book.org/4ed/Reflection_Models/Diffuse_Reflection',
 'derived_own_equations':{
  'incident_fluence_plane_ph_m2':'N*h/[4*pi*(r^2+h^2)^(3/2)]',
  'Lambert_reflected_radiance_fluence_ph_m2_sr':'rho/pi timesincidentfluence',
  'reflected_detector_fraction_dN_over_N':'rho*A*h*H^2*r*dr/[2*pi*(r^2+h^2)^(3/2)*(r^2+H^2)^2]',
  'direct_aperture_fraction':'A/[4*pi*(H-h)^2]',
  'reflected_to_direct_integrand_dr':'2*rho*h*H^2*(H-h)^2*r/[(r^2+h^2)^(3/2)*(r^2+H^2)^2]',
  'substitution':'x=h/sqrt(r^2+h^2); integratex_min...1 withweight2*rho*H^2*(H-h)^2/(r(x)^2+H^2)^2',
  'delay_seconds':'[sqrt(h^2+r^2)+sqrt(H^2+r^2)-(H-h)]/c',
  'surface_limit_h_to_zero':'Nadirfluencefactor1+2*rho; hemisphereescape(1+rho)/2',
  'convolution':'y(t)=[s(t)+integralweight(x)*s(t-delay(x))dx]/(1+B); s(t)clippedlinearrampwithindefiniteplateau',
  'rise_definition':'t_at_90percent_finalplateau minus t_at_10percent_finalplateau. NOTarrival-timequantiles, FWHM oronsettopeak.'},
 'assumptions':['Isotropicpointsourceaboveplane; finiteaperturesmallcomparedwithrange','InstantaneouslocalLambertreflectionwithassumedrho, noresidencetime/lateraltransportinsidecloud','Infiniteflatplaneandunlimitedradiancefield; detectorcircularFOVtruncatesreceivedreflection','Vacuumaboveplane; noRayleigh,aerosol,ozone,Earthcurvature,sourceextent,fieldresponse/electronics/sampling','Infiniteplateauafterintrinsiclinear-rise; nofinite3mspulseendordroop','No mappingofrhotooverheadcloudtau. h=0 meanslimitingjust-abovesurfaceposition, notanembeddedcloud-topscatteringboundary'],
 'interpretation_limits':['Morelowerreflectancemayaddphotons; thisdoesnotmeanmoreOVERHEADextinctionmustincreasefluence','Factor3forrho1isthistoynadirlimit, notuniversalMie/icecloudbound','Photonhemisphereescapefraction<=1evenwhendirectionalratio>1','Geometry-onlyrisechangesarenotclaimedrealdepthcalibrationorvalidationof30–80versus200usclassifier','Raisingtelescopealtitudemainlychangesdilution; pulsebroadeningdependsonpathdifferences, notcommonflighttime'],
 'validation':{'rows':len(rows),'independent_radial_coordinate_check':True,'cumulative_escape_bound_passed':True,'threshold_root_checks_passed':True},
 'selected_examples':[r for r in rows if r['assumed_Lambert_albedo']==1 and r['full_circular_FOV_deg']==60 and r['assumed_intrinsic_10_90_linear_rise_us']==30]}
(p/'derivation_and_qualification.json').write_text(json.dumps(derivation,indent=2),encoding='utf-8')
print(json.dumps({'rows':len(rows),'examples':derivation['selected_examples']},indent=2))