Review index · Download raw source
reference/cloud/engine_snapshot.py
1"""Independent scalar finite-cylinder Monte Carlo. NumPy only.2Original implementation for this study; not CloudScat source.3See README.md and QA.json for validation scope and numerical limits.4"""5from dataclasses import dataclass,asdict6from pathlib import Path7import argparse,hashlib,json,math,time8import numpy as np9C=299792458.1011def sha(path):12 return hashlib.sha256(Path(path).read_bytes()).hexdigest()1314class Phase:15 def __init__(self,kind,mu=None,f=None,g=None):16 self.kind=kind17 if kind=="table":18 self.mu=np.asarray(mu,float); self.f=np.asarray(f,float)19 self.raw_integral=float(np.trapezoid(self.f,self.mu))20 self.f=self.f/self.raw_integral21 self.dx=np.diff(self.mu); self.slope=np.diff(self.f)/self.dx22 self.mass=.5*(self.f[:-1]+self.f[1:])*self.dx23 self.cdf=np.r_[0.,np.cumsum(self.mass)];self.cdf[-1]=1.24 # Exact first moment of each piecewise-linear segment.25 x=self.mu[:-1];d=self.dx;y=self.f[:-1];s=self.slope26 self.g=float(np.sum(x*y*d+(x*s+y)*d*d/2+s*d**3/3))27 else:28 self.g=float(g)29 def pdf_mu(self,x):30 x=np.clip(x,-1.,1.)31 if self.kind=="table":return np.interp(x,self.mu,self.f)32 g=self.g33 return (1-g*g)/(2*(1+g*g-2*g*x)**1.5)34 def pdf_sr(self,x):return self.pdf_mu(x)/(2*math.pi)35 def sample_mu(self,rng,n):36 return self.inverse(rng.random(n))37 def inverse(self,u):38 u=np.asarray(u)39 if self.kind!="table":40 g=self.g41 if abs(g)<1e-12:return 2*u-142 return np.clip((1+g*g-((1-g*g)/(1-g+2*g*u))**2)/(2*g),-1,1)43 k=np.minimum(np.searchsorted(self.cdf,u,side="right")-1,len(self.dx)-1)44 k=np.maximum(k,0);v=u-self.cdf[k];y=self.f[k];s=self.slope[k]45 disc=np.maximum(0,y*y+2*s*v)46 den=y+np.sqrt(disc)47 d=np.divide(2*v,den,out=np.zeros_like(v,dtype=float),where=den>0)48 return np.clip(self.mu[k]+np.clip(d,0,self.dx[k]),-1,1)4950def load_phases(lut,reff=10.,wavelength=337.1):51 z=np.load(lut)52 i=int(np.argmin(abs(z["reff_um"]-reff)))53 if abs(z["reff_um"][i]-reff)>1e-9:raise ValueError("Exact requested reff unavailable")54 a,b=z["wavelength_nm"];v=(wavelength-a)/(b-a)55 if not 0<=v<=1:raise ValueError("Not a bracketed interpolation")56 mu=np.cos(np.deg2rad(z["theta_deg"]))[::-1].copy()57 f=((1-v)*z["phase"][i,0,0]+v*z["phase"][i,1,0])[::-1].astype(float)58 tab=Phase("table",mu,f)59 meta={"input_file_name":Path(lut).name,"input_sha256":sha(lut),60 "original_full_nc_sha256":"2b37e27cffbbf0dd7e048df16c5fb612cd7a83e3c657c851bb50e758e11e870e",61 "reff_um":reff,"wavelength_nm":wavelength,"bracket_nm":[float(a),float(b)],62 "linear_weight_upper":float(v),"native_angle_nodes":len(mu),63 "raw_P11_integral_dmu":tab.raw_integral,64 "normalization":"f(mu)=P11/integral(P11 dmu); p_per_sr=f/(2pi)",65 "representation":"Linear interpolation of P11 in wavelength, then piecewise-linear probability density in mu; each inverse CDF segment solved quadratically.",66 "derived_g":tab.g,"scalar_SSA_assumed":1.,67 "input_source_url":"https://docs.hygeos.com/public.php/webdav/ic_baum_ghm_sol.nc",68 "registry_url":"https://raw.githubusercontent.com/hygeos/smartg/master/smartg/auxdata.py",69 "raw_LUT_redistributed":False,70 "license_status":"Original LUT is not included in release; permission for redistribution not established.",71 "physical_limits":"Randomly oriented severely rough bulk ice mixture; processed-LUT lineage and fresh-anvil PSD applicability remain uncertain. No vector polarization propagation or measured337nm absorption."}72 return {"ghm":tab,"hg":Phase("hg",g=tab.g)},meta7374@dataclass75class Geometry:76 base_m:float=14000.77 top_m:float=17000.78 radius_m:float=30000.79 observer_z_m:float=550000.80 full_fov_deg:float=60.81 @property82 def thickness(self):return self.top_m-self.base_m8384def intervals(pos,direction,geo):85 """Forward ray interval inside convex cylinder; no physical epsilon."""86 p=np.asarray(pos);u=np.asarray(direction);n=len(p)87 uz=u[:,2];z=p[:,2]88 moving=abs(uz)>1e-1589 az=np.full(n,-np.inf);bz=np.full(n,np.inf)90 t0=np.divide(geo.base_m-z,uz,out=np.zeros(n),where=moving)91 t1=np.divide(geo.top_m-z,uz,out=np.zeros(n),where=moving)92 az[moving]=np.minimum(t0,t1)[moving];bz[moving]=np.maximum(t0,t1)[moving]93 outz=(~moving)&((z<geo.base_m)|(z>geo.top_m));az[outz]=np.inf;bz[outz]=-np.inf94 aa=u[:,0]**2+u[:,1]**295 bb=p[:,0]*u[:,0]+p[:,1]*u[:,1]96 cc=p[:,0]**2+p[:,1]**2-geo.radius_m**297 disc=bb*bb-aa*cc98 radial=aa>1e-2899 ar=np.full(n,-np.inf);br=np.full(n,np.inf)100 good=radial&(disc>=0)101 sq=np.sqrt(np.maximum(disc,0))102 ar[good]=(-bb[good]-sq[good])/aa[good]103 br[good]=(-bb[good]+sq[good])/aa[good]104 bad=(radial&(disc<0))|((~radial)&(cc>0))105 ar[bad]=np.inf;br[bad]=-np.inf106 lo=np.maximum(0,np.maximum(az,ar));hi=np.minimum(bz,br)107 hit=hi>lo108 return lo,hi,hit109110def segment_inside(pos,direction,length,geo):111 lo,hi,hit=intervals(pos,direction,geo)112 return np.where(hit,np.maximum(0,np.minimum(hi,length)-lo),0)113114def rotate(u,mu,phi):115 n=len(u)116 a=np.sqrt(u[:,0]**2+u[:,1]**2)117 e=np.zeros((n,3));nz=a>1e-12118 e[nz,0]=-u[nz,1]/a[nz];e[nz,1]=u[nz,0]/a[nz];e[~nz,0]=1.119 f=np.cross(u,e)120 transverse=np.sqrt(np.maximum(0,1-mu*mu))121 out=mu[:,None]*u+transverse[:,None]*(np.cos(phi)[:,None]*e+np.sin(phi)[:,None]*f)122 # Roundoff only; physical directions should remain unit vectors.123 return out/np.linalg.norm(out,axis=1)[:,None]124125def launch(rng,n,phase,importance):126 u=rng.uniform(-1,1,n);phi=rng.uniform(0,2*math.pi,n)127 if importance:128 choose=rng.random(n)<importance129 u[choose]=phase.sample_mu(rng,int(choose.sum()))130 weight=1/(1-importance+importance*4*math.pi*phase.pdf_sr(u))131 else:weight=np.ones(n)132 st=np.sqrt(np.maximum(0,1-u*u))133 return np.column_stack((st*np.cos(phi),st*np.sin(phi),u)),weight134135def simulate(phase,source_m,total_tau,n,seed,geo=None,dt_s=2e-6,max_delay_s=.012,136 max_scatter=20000,importance=.5):137 """Per-emitted-photon plane-aperture fluence; vacuum exterior, SSA1.138 Initial physical source remains isotropic. Optional source-only proposal139 q=(1-a)/(4pi)+a*p_phase(mu_to_nadir) uses exact likelihood weights.140 Collision directions ALWAYS use the physical phase, identically evaluated141 in the next-event estimator. No biased walk or implicit absorption.142 """143 geo=geo or Geometry();start=time.perf_counter();rng=np.random.default_rng(seed)144 beta=total_tau/geo.thickness145 r0=geo.observer_z_m-source_m146 vacuum=1/(4*math.pi*r0*r0)147 src=np.array([[0.,0.,source_m]]);vertical=np.array([[0.,0.,1.]])148 direct_length=float(segment_inside(src,vertical,np.array([r0]),geo)[0])149 direct=math.exp(-beta*direct_length)150 nbins=int(round(max_delay_s/dt_s));hist=np.zeros(nbins+1)151 imp=importance if geo.base_m<source_m<geo.top_m else 0.152 dirs,weights=launch(rng,n,phase,imp)153 launch_sum=float(weights.sum())154 positions=np.tile(src,(n,1));lo,hi,hit=intervals(positions,dirs,geo)155 if beta==0:hit[:]=False156 escaped_initial=float(weights[~hit].sum())157 escaped_weight=escaped_initial;escaped_n=int((~hit).sum())158 escaped_up=float(weights[(~hit)&(dirs[:,2]>0)].sum())159 boundary_weights={"top":0.,"base":0.,"side":0.,"miss_or_vacuum":escaped_initial}160 ids=np.flatnonzero(hit)161 pos=positions[hit]+lo[hit,None]*dirs[hit]162 u=dirs[hit];w=weights[hit];path=lo[hit].copy()163 nee=np.zeros(n);collisions=0;iteration=0;min_delay=np.inf;negative=0164 first_hist=np.zeros(nbins+1)165 cone_weights={3.:0.,5.:0.,10.:0.}166 def cone_tally(positions,directions,ww):167 up=directions[:,2]>0168 if not np.any(up):return169 pp=positions[up];dd=directions[up];wt=ww[up]170 travel=(geo.observer_z_m-pp[:,2])/dd[:,2]171 xx=pp[:,0]+travel*dd[:,0];yy=pp[:,1]+travel*dd[:,1]172 rr=xx*xx+yy*yy173 for deg in cone_weights:174 radius=r0*math.tan(math.radians(deg))175 cone_weights[deg]+=float(wt[rr<=radius*radius].sum())176 max_dir_error=0.;weighted_path_exit=0.;exit_max_path=0.177 while len(ids) and iteration<max_scatter:178 iteration+=1179 _,exit_distance,ok=intervals(pos,u,geo)180 if not np.all(ok):181 raise RuntimeError("Interior ray lost cylinder intersection")182 free=rng.exponential(1/beta,len(ids))183 collide=free<exit_distance184 leaving=~collide185 if np.any(leaving):186 ww=w[leaving];endpos=pos[leaving]+exit_distance[leaving,None]*u[leaving]187 if iteration>1:cone_tally(endpos,u[leaving],ww)188 plen=path[leaving]+exit_distance[leaving]189 escaped_weight+=float(ww.sum());escaped_n+=int(leaving.sum())190 escaped_up+=float(ww[u[leaving,2]>0].sum())191 weighted_path_exit+=float(np.sum(ww*plen));exit_max_path=max(exit_max_path,float(plen.max()))192 top=abs(endpos[:,2]-geo.top_m)<1e-5193 base=(~top)&(abs(endpos[:,2]-geo.base_m)<1e-5)194 side=~(top|base)195 for label,mask in [("top",top),("base",base),("side",side)]:196 boundary_weights[label]+=float(ww[mask].sum())197 if not np.any(collide):198 ids=ids[:0];break199 ids=ids[collide];pos=pos[collide]+free[collide,None]*u[collide]200 path=path[collide]+free[collide];u=u[collide];w=w[collide]201 collisions+=len(ids)202 toward=-pos.copy();toward[:,2]+=geo.observer_z_m203 distance=np.linalg.norm(toward,axis=1);direction=toward/distance[:,None]204 cosine=np.einsum("ij,ij->i",u,direction)205 optical_length=segment_inside(pos,direction,distance,geo)206 aperture_cos=direction[:,2]207 accepted=aperture_cos>=math.cos(math.radians(geo.full_fov_deg/2))208 # Dimensionless relative to isotropic vacuum aperture fluence.209 contribution=w*4*math.pi*r0*r0*phase.pdf_sr(cosine)*aperture_cos/(distance*distance)*np.exp(-beta*optical_length)*accepted210 delay=(path+distance-r0)/C211 min_delay=min(min_delay,float(delay.min()))212 negative+=int(np.sum(delay < -1e-12))213 if negative:raise RuntimeError("Noncausal next-event contribution beyond roundoff")214 delay=np.maximum(delay,0)215 bins=np.minimum((delay/dt_s).astype(np.int64),nbins)216 part=np.bincount(bins,weights=contribution,minlength=nbins+1)217 hist+=part218 if iteration==1:first_hist+=part219 nee[ids]+=contribution220 mu=phase.sample_mu(rng,len(ids));phi=rng.uniform(0,2*math.pi,len(ids))221 u=rotate(u,mu,phi)222 max_dir_error=max(max_dir_error,float(np.max(abs(np.linalg.norm(u,axis=1)-1))))223 truncated_n=len(ids);truncated_weight=float(w.sum()) if truncated_n else 0.224 residual=launch_sum-escaped_weight-truncated_weight225 info={"seed":int(seed),"n_packets":n,"total_tau":float(total_tau),226 "overlying_tau":beta*direct_length,"source_altitude_m":source_m,227 "phase":phase.kind,"g":phase.g,"initial_importance_fraction":imp,228 "vacuum_transfer_m_inverse2":vacuum,229 "direct_relative_to_vacuum":direct,230 "scattered_relative_to_vacuum":float(hist.sum()/n),231 "total_relative_to_vacuum":direct+float(hist.sum()/n),232 "weighted_launch_mean":launch_sum/n,233 "weighted_escape_fraction":escaped_weight/n,"weighted_upward_escape_fraction":escaped_up/n,234 "weighted_truncation_fraction":truncated_weight/n,235 "energy_balance_residual_per_packet":residual/n,236 "escaped_packet_count":escaped_n,"truncated_packet_count":truncated_n,237 "boundary_escape_weights_per_packet":{k:v/n for k,v in boundary_weights.items()},238 "mean_collisions_per_packet":collisions/n,"max_scatter_iterations":iteration,239 "weighted_insideplus_entry_path_exit_mean_m":weighted_path_exit/n,240 "maximum_entryplus_cloud_path_m":exit_max_path,241 "min_scattered_delay_s":None if not np.isfinite(min_delay) else min_delay,242 "negative_delay_beyond_tolerance_count":negative,243 "max_direction_norm_error":max_dir_error,244 "response_overflow_relative_to_vacuum":float(hist[-1]/n),245 "first_scatter_relative_to_vacuum":float(first_hist.sum()/n),246 "finite_receiver_scattered_gain":{str(k):4*v/n/math.tan(math.radians(k))**2 for k,v in cone_weights.items()},247 "single_packet_standard_error_relative":float(np.std(nee,ddof=1)/math.sqrt(n)),248 "elapsed_s":time.perf_counter()-start}249 return info,hist/n250251def phase_qa(phase,seed=821013,n=300000):252 rng=np.random.default_rng(seed);mu=phase.sample_mu(rng,n)253 # Independent Gauss integration of piecewise-linear input.254 if phase.kind=="table":255 x,w=np.polynomial.legendre.leggauss(8)256 xx=(phase.mu[1:]+phase.mu[:-1])[:,None]/2+np.diff(phase.mu)[:,None]*x/2257 ww=np.diff(phase.mu)[:,None]*w/2258 f=phase.pdf_mu(xx)259 integrals=[float(np.sum(ww*f)),float(np.sum(ww*f*xx)),float(np.sum(ww*f*(3*xx*xx-1)/2))]260 else:261 integrals=[1.,phase.g,phase.g**2]262 sample=[float(mu.mean()),float(np.mean((3*mu*mu-1)/2))]263 se=[float(mu.std(ddof=1)/math.sqrt(n)),float(np.std((3*mu*mu-1)/2,ddof=1)/math.sqrt(n))]264 u=np.linspace(.000001,.999999,10000);inv=phase.inverse(u)265 if phase.kind=="table":266 k=np.minimum(np.searchsorted(phase.mu,inv,side="right")-1,len(phase.dx)-1);k=np.maximum(0,k)267 dx=inv-phase.mu[k]268 back=phase.cdf[k]+phase.f[k]*dx+.5*phase.slope[k]*dx*dx269 inverse_error=float(np.max(abs(back-u)))270 else:271 g=phase.g272 back=(1-g*g)/(2*g)*(1/np.sqrt(1+g*g-2*g*inv)-1/(1+g))273 inverse_error=float(np.max(abs(back-u)))274 return {"normalization":integrals[0],"g_integral":integrals[1],"g_parameter":phase.g,275 "second_Legendre_moment":integrals[2],"sample_size":n,"sample_moments":sample,276 "sample_standard_errors":se,"sample_z_scores":[(sample[i]-integrals[i+1])/se[i]for i in range(2)],277 "inverse_CDF_max_abs_error":inverse_error,278 "sample_in_bounds":bool(np.all((mu>=-1)&(mu<=1))),279 "pass":bool(abs(integrals[0]-1)<1e-10 and inverse_error<1e-9 and all(abs((sample[i]-integrals[i+1])/se[i])<5 for i in range(2)))}280281if __name__=="__main__":282 ap=argparse.ArgumentParser();ap.add_argument("--lut",required=True);ap.add_argument("--out",required=True)283 ap.add_argument("--pilot",action="store_true");args=ap.parse_args()284 out=Path(args.out);out.mkdir(exist_ok=True)285 phases,meta=load_phases(args.lut)286 q={"phase_input":meta,"phases":{k:phase_qa(v)for k,v in phases.items()}}287 (out/"phase_QA.json").write_text(json.dumps(q,indent=2),encoding="utf-8")288 cases=[]289 for name in ["ghm","hg"]:290 for altitude,tau in [(15000.,1.),(15000.,45.),(17000.,10.),(18000.,10.)]:291 info,h=simulate(phases[name],altitude,tau,5000,810022+len(cases))292 cases.append(info);print(json.dumps(info),flush=True)293 (out/"pilot_results.json").write_text(json.dumps(cases,indent=2),encoding="utf-8")294