Review index · Download raw source
reference/analytic/diffusion.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 brentq6p=Path(r'D:/Photometer_research_20261022/revision_20261022_r4/analysis/analytic/diffusion');p.mkdir(exist_ok=True)7c=299792458.8def normalized_at_mode_units(x):return math.exp(-1.5*math.log(x)-1.5*(1/x-1)) if x>0 else 0.9x10=brentq(lambda x:normalized_at_mode_units(x)-.1,1e-7,1)10x90=brentq(lambda x:normalized_at_mode_units(x)-.9,1e-7,1)11coef=(x90-x10)/612rows=[]13for D in [3e9,3.7e9,7.9e10]:14 for rise_us in [5.5,30.,80.,200.]:15 depth=math.sqrt(D*rise_us*1e-6/coef);mode=depth*depth/(6*D)16 for gamma in [0.,.1,10.,50.]:17 # L is inferred under gamma0, then absorption is varied conditionally.18 tpeak=mode if gamma==0 else 2*(depth*depth/(4*D))/(1.5+math.sqrt(2.25+gamma*depth*depth/D))19 def logh(t):return math.log(depth)-.5*math.log(4*math.pi*D)-1.5*math.log(t)-depth*depth/(4*D*t)-gamma*t20 lp=logh(tpeak)21 t10=brentq(lambda t:logh(t)-lp-math.log(.1),tpeak*1e-5,tpeak,xtol=1e-16)22 t90=brentq(lambda t:logh(t)-lp-math.log(.9),tpeak*1e-5,tpeak,xtol=1e-16)23 absorption_scale=gamma*depth*depth/(4*D)24 density=lambda u:2/math.sqrt(math.pi)*math.exp(-u*u-absorption_scale/(u*u)) if u>0 else 0.25 total=math.exp(-depth*math.sqrt(gamma/D))26 numeric=quad(density,0,math.inf,epsabs=1e-11,epsrel=1e-10)[0]27 assert abs(numeric-total)<1e-928 low=depth/math.sqrt(4*D*.003)29 captured=quad(density,low,math.inf,epsabs=1e-11,epsrel=1e-10)[0]30 if gamma==0:assert abs(captured-math.erfc(low))<1e-1031 rows.append(dict(assumed_D_m2_s=D,input_nonabsorbing_10_90_rise_us=rise_us,inferred_depth_m=depth,32 assumed_uniform_absorption_frequency_s=gamma,mode_onset_to_peak_us=tpeak*1e6,33 leading_10_90_us=(t90-t10)*1e6,earliest_physical_straight_path_us=depth/c*1e6,34 diffusion_t10_us=t10*1e6,diffusion_t10_precedes_light_flight=t10<depth/c,35 equivalent_transport_mean_free_path_m=3*D/c,equivalent_transport_depth=depth/(3*D/c),36 all_time_diffusion_exit_probability=total,within3ms_exit_probability=captured,37 fraction_of_eventually_exiting_weight_within3ms=captured/total))38with (p/'conditional_diffusion_rise_absorption.csv').open('w',newline='',encoding='utf-8') as f:39 w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)40times=[1e-6*(1.12**i) for i in range(95)]41family=[]42for depth,D in [(1000.,3e9),(2000.,1.2e10),(3000.,2.7e10)]:43 h=[depth/math.sqrt(4*math.pi*D)*t**-1.5*math.exp(-depth*depth/(4*D*t)) for t in times]44 family.append({'depth_m':depth,'D_m2_s':D,'if_cloud_top18km_source_altitude_km':18-depth/1000,'L2_over_D_s':depth*depth/D,'sampled_kernel_s_inverse':h})45for entry in family[1:]:assert max(abs(x-y)/max(1e-100,abs(y)) for x,y in zip(entry['sampled_kernel_s_inverse'],family[0]['sampled_kernel_s_inverse']))<1e-1046report={'recorded_utc':datetime.datetime.now(datetime.timezone.utc).isoformat(),47 'status':'EXACT_MATH_FOR_IDEAL_DIFFUSION_MODEL_NOT_REAL_CLOUD_VALIDATION',48 'questions':['A0','A2','F18','G21','G22'],49 'existing_primary_recheck':{'doi':'10.1029/2022JD037460','url':'https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2022JD037460','locator':'Bai2023sections6–7','access':'PrimaryfullHTMLrechecked; Eq6bitmapstandardlinkcachemiss, noimageclaimedinspected. Sourcealreadyreviewedearlier; thisisnotnewpaper.','reported_relation':'Tr10–90=.0742*L²/D,mode=L²/(6D)','D_range_context':'3.7e9–7.9e10modelparametersfromPSDdiscussion; notmeasuredprobabilitydistribution.'},50 'own_first_passage_derivation':{'PDE':'d_t u=D*d_zz u-gamma*u onz>0; absorbingplanez0; initialdeltaatzL; methodofimages','plane_first_exit_density':'h(t)=L/sqrt(4*pi*D)*t^(-3/2)*exp[-L²/(4Dt)-gamma*t]','nonabsorbing_mode':'L²/(6D)','leading_shape_at_t_over_mode_x':'h/hmode=x^(-3/2)*exp[-1.5*(1/x-1)]','x10':x10,'x90':x90,'computed_10_90_coefficient_times_L2_over_D':coef,'rise_over_onset_peak':x90-x10,'whole_kernel_identifiability':'SameL²/Dandgamma produceidenticalwholeh(t), notjustsamemode. LcannotbedeterminedwithoutD/source/boundaryinformation.','nonabsorbing_cumulative_exit':'erfc[L/sqrt(4Dt)]','all_time_exit_with_uniform_absorption':'exp[-L*sqrt(gamma/D)]'},51 'counterexample_family':family,'counterexample_times_s':times,52 'scope_limits':['Parabolicdiffusionhasacausalnonzerot<t_lightprecursor; report flagswhen10%thresholditselfisunphysical. Notanearly-timeRTEbenchmark.','Semi-infinitehomogeneousmediumwithabsorbingplaneandinstantaneousdeltaemission; nofinitecloudbase,phasefunction,sourceextent,upwardballistictermorinternalcloudtopboundary.','A>3msobservedpulseisnotassumed3msintrinsicemissionordwelltime; realemissionrequiresconvolution','Gammauniformeverywhere, notspatialozone/Rayleigh/iceprofile; example0.1/10/50notempiricalconfidenceintervals.','All-timeescapeand3mscapturebelongtothismodel; noactualNBEenergycorrection','ModelDvaluesandtop18kmincounterexampleareconditional, notindependentheightobservations','Published5.5usmodel-basedsubsamplefitdoesnotprove general5.5usmeasurementprecisionat10usbinning'],53 'validation':{'rows':len(rows),'absorption_exit_integral_vs_closed_formula':len(rows),'nonabsorbing_CDF_checks':12,'matched_family_sampled_kernel_checks':190,'known_published_rounded_coefficient_difference':coef-.0742},54 'nonabsorbing_3ms_capture_examples':[r for r in rows if r['assumed_D_m2_s']==3e9 and r['assumed_uniform_absorption_frequency_s']==0]}55(p/'derivation_and_qualification.json').write_text(json.dumps(report,indent=2),encoding='utf-8')56print(json.dumps({'coefficient':coef,'rise_over_peak':x90-x10,'rows':len(rows),'3ms_examples':report['nonabsorbing_3ms_capture_examples']},indent=2))57