Review index · Download raw source
reference/analytic/reflection.py
1from pathlib import Path2import csv,datetime,json,math,sys3sys.path.insert(0,r'D:/Photometer_research_20261022/research/cloud_transport/private_readers')4from scipy.integrate import quad5from scipy.optimize import brentq67p=Path(r'D:/Photometer_research_20261022/revision_20261022_r4/analysis/analytic/reflection')8p.mkdir(exist_ok=True)9c=299792458.;area=math.pi*(.025/2)**210def integrate_x(f,lo,hi,H,h):11 points=sorted({x for x in [h/H*.1,h/H,h/H*10,h/H*100,.001,.01,.1] if lo<x<hi})12 return quad(f,lo,hi,points=points,epsabs=1e-11,epsrel=1e-10,limit=200)[0]13def setup(top_km,h_m,albedo,fov):14 H=(550-top_km)*1000; h=h_m15 if h==0:16 return H,1.,2*albedo,lambda x:0.,lambda x:2*albedo17 xmax_radius=H*math.tan(math.radians(fov/2))18 xmin=h/math.hypot(h,xmax_radius)19 def radius(x):return h*math.sqrt(max(0.,1/x**2-1))20 def delay(x):21 r=radius(x)22 # Stablepositiveexpression avoids subtractingtwo~500km lengths.23 return (h/x+h+r*r/(math.hypot(H,r)+H))/c24 def weight(x):25 r=radius(x)26 return 2*albedo*H*H*(H-h)**2/(H*H+r*r)**227 B=integrate_x(weight,xmin,1,H,h)28 return H,xmin,B,delay,weight2930def rise(top,h,a,fov,input_rise_us):31 H,xmin,B,delay,weight=setup(top,h,a,fov)32 ramp=input_rise_us*1e-6/.833 if h==0 or a==0:return B,input_rise_us34 dmin=delay(1);dmax=delay(xmin)35 def boundary(t):36 if t<=dmin:return 1.37 if t>=dmax:return xmin38 return brentq(lambda x:delay(x)-t,xmin,1.,xtol=1e-14)39 def received(t):40 direct=min(1.,max(0.,t/ramp))41 lo=boundary(t);hi=boundary(t-ramp)42 middle=integrate_x(lambda x:weight(x)*(t-delay(x))/ramp,lo,hi,H,h) if hi>lo else 0.43 flat=integrate_x(weight,hi,1.,H,h) if hi<1 else 0.44 return (direct+middle+flat)/(1+B)45 t10=brentq(lambda t:received(t)-.1,0,ramp+dmax,xtol=1e-12)46 t90=brentq(lambda t:received(t)-.9,0,ramp+dmax,xtol=1e-12)47 assert abs(received(t10)-.1)<1e-6 and abs(received(t90)-.9)<1e-648 return B,(t90-t10)*1e64950rows=[]51for top,h_values in [(15,[0,1,100,1000,3000]),(17,[0,1,100,1000])]:52 for h in h_values:53 for albedo in [.8,1.]:54 for fov in [60,90]:55 H,xmin,B,delay,weight=setup(top,h,albedo,fov)56 incident_fraction_in_view=.5*(1-xmin) if h else .557 # Independent radialcoordinate quadrature checks exactweight/Jacobian.58 if h:59 R=H*math.tan(math.radians(fov/2))60 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]61 if abs(ratio_radial-B)>=2e-7: raise RuntimeError(f'coordinatecheckfailedtop={top},h={h},rho={albedo},fov={fov},B={B},radial={ratio_radial}')62 else:ratio_radial=B63 for tr in [30.,80.,200.]:64 boost_part,observed=rise(top,h,albedo,fov,tr)65 assert 0<=B<=2*albedo+1e-966 rows.append(dict(plane_altitude_km=top,source_above_plane_m=h,source_altitude_km=top+h/1000,67 satellite_altitude_km=550,full_circular_FOV_deg=fov,assumed_Lambert_albedo=albedo,68 nadir_fluence_ratio_to_vacuum=1+B,reflected_to_direct_ratio=B,independent_radial_integral=ratio_radial,69 downward_source_fraction_intercepted_within_view=incident_fraction_in_view,70 all_directions_total_upward_escape_fraction=(1+albedo)/2,71 ideal_25mm_aperture_fraction_of_emitted_photons=area/(4*math.pi*(H-h)**2)*(1+B),72 minimum_reflected_extra_delay_us=0 if h==0 else 2*h/c*1e6,73 maximum_reflected_extra_delay_in_FOV_us=0 if h==0 else delay(xmin)*1e6,74 assumed_intrinsic_10_90_linear_rise_us=tr,observed_10_90_of_final_plateau_rise_us=observed,75 field_observer='nadir, horizontalflataperture; directpointsourceontheaxis',76 status='OWNEXACTSURFACETOYNOTCLOUDRT'))77with (p/'conditional_reflection_and_rise.csv').open('w',newline='',encoding='utf-8') as f:78 w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)79derivation={80 'utc':datetime.datetime.now(datetime.timezone.utc).isoformat(),81 'status':'PRIVATE_ANALYTIC_DIAGNOSTIC_NOT_A_REPLACEMENT_CLOUDSCAT_RUN',82 'question_tags':['A0','A2','G22'],83 'purpose':'Showwhydirectionalfluencecanexceedvacuumwithoutcreatingphotons, andwhetherreflectedgeometrycanbroadenanevenabove-topsourceunderoneexplicitidealmodel.',84 'primary_BRDF_reference':'https://www.pbr-book.org/4ed/Reflection_Models/Diffuse_Reflection',85 'derived_own_equations':{86 'incident_fluence_plane_ph_m2':'N*h/[4*pi*(r^2+h^2)^(3/2)]',87 'Lambert_reflected_radiance_fluence_ph_m2_sr':'rho/pi timesincidentfluence',88 '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]',89 'direct_aperture_fraction':'A/[4*pi*(H-h)^2]',90 '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]',91 'substitution':'x=h/sqrt(r^2+h^2); integratex_min...1 withweight2*rho*H^2*(H-h)^2/(r(x)^2+H^2)^2',92 'delay_seconds':'[sqrt(h^2+r^2)+sqrt(H^2+r^2)-(H-h)]/c',93 'surface_limit_h_to_zero':'Nadirfluencefactor1+2*rho; hemisphereescape(1+rho)/2',94 'convolution':'y(t)=[s(t)+integralweight(x)*s(t-delay(x))dx]/(1+B); s(t)clippedlinearrampwithindefiniteplateau',95 'rise_definition':'t_at_90percent_finalplateau minus t_at_10percent_finalplateau. NOTarrival-timequantiles, FWHM oronsettopeak.'},96 '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'],97 'interpretation_limits':['Morelowerreflectancemayaddphotons; thisdoesnotmeanmoreOVERHEADextinctionmustincreasefluence','Factor3forrho1isthistoynadirlimit, notuniversalMie/icecloudbound','Photonhemisphereescapefraction<=1evenwhendirectionalratio>1','Geometry-onlyrisechangesarenotclaimedrealdepthcalibrationorvalidationof30–80versus200usclassifier','Raisingtelescopealtitudemainlychangesdilution; pulsebroadeningdependsonpathdifferences, notcommonflighttime'],98 'validation':{'rows':len(rows),'independent_radial_coordinate_check':True,'cumulative_escape_bound_passed':True,'threshold_root_checks_passed':True},99 '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]}100(p/'derivation_and_qualification.json').write_text(json.dumps(derivation,indent=2),encoding='utf-8')101print(json.dumps({'rows':len(rows),'examples':derivation['selected_examples']},indent=2))102