Review index · Download raw source
reference/cloud/cylinder_mc.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,cone_deg=None):126 u=rng.uniform(-1,1,n);phi=rng.uniform(0,2*math.pi,n)127 if importance:128 choose=rng.random(n)<importance129 if cone_deg is None:130 u[choose]=phase.sample_mu(rng,int(choose.sum()))131 weight=1/(1-importance+importance*4*math.pi*phase.pdf_sr(u))132 else:133 if not 0<cone_deg<90:raise ValueError("Importance cone must be between0 and90 degrees")134 cut=math.cos(math.radians(cone_deg))135 u[choose]=rng.uniform(cut,1,int(choose.sum()))136 weight=1/(1-importance+importance*2/(1-cut)*(u>=cut))137 else:weight=np.ones(n)138 st=np.sqrt(np.maximum(0,1-u*u))139 return np.column_stack((st*np.cos(phi),st*np.sin(phi),u)),weight140141def simulate(phase,source_m,total_tau,n,seed,geo=None,dt_s=2e-6,max_delay_s=.012,142 max_scatter=20000,importance=.5,importance_cone_deg=None):143 """Per-emitted-photon plane-aperture fluence; vacuum exterior, SSA1.144 Initial physical source remains isotropic. Optional source-only proposal145 q=(1-a)/(4pi)+a*p_phase(mu_to_zenith) uses exact likelihood weights.146 Collision directions ALWAYS use the physical phase, identically evaluated147 in the next-event estimator. No biased walk or implicit absorption.148 """149 geo=geo or Geometry();start=time.perf_counter();rng=np.random.default_rng(seed)150 beta=total_tau/geo.thickness151 r0=geo.observer_z_m-source_m152 vacuum=1/(4*math.pi*r0*r0)153 src=np.array([[0.,0.,source_m]]);vertical=np.array([[0.,0.,1.]])154 direct_length=float(segment_inside(src,vertical,np.array([r0]),geo)[0])155 direct=math.exp(-beta*direct_length)156 nbins=int(round(max_delay_s/dt_s));hist=np.zeros(nbins+1)157 imp=importance if beta>0 and geo.base_m<source_m<geo.top_m else 0.158 dirs,weights=launch(rng,n,phase,imp,importance_cone_deg)159 launch_sum=float(weights.sum())160 positions=np.tile(src,(n,1));lo,hi,hit=intervals(positions,dirs,geo)161 if beta==0:hit[:]=False162 escaped_initial=float(weights[~hit].sum())163 escaped_weight=escaped_initial;escaped_n=int((~hit).sum())164 escaped_up=float(weights[(~hit)&(dirs[:,2]>0)].sum())165 boundary_weights={"top":0.,"base":0.,"side":0.,"miss_or_vacuum":escaped_initial}166 ids=np.flatnonzero(hit)167 pos=positions[hit]+lo[hit,None]*dirs[hit]168 u=dirs[hit];w=weights[hit];path=lo[hit].copy()169 nee=np.zeros(n);collisions=0;iteration=0;min_delay=np.inf;negative=0170 first_hist=np.zeros(nbins+1)171 cone_weights={3.:0.,5.:0.,10.:0.}172 def cone_tally(positions,directions,ww):173 up=directions[:,2]>0174 if not np.any(up):return175 pp=positions[up];dd=directions[up];wt=ww[up]176 travel=(geo.observer_z_m-pp[:,2])/dd[:,2]177 xx=pp[:,0]+travel*dd[:,0];yy=pp[:,1]+travel*dd[:,1]178 rr=xx*xx+yy*yy179 for deg in cone_weights:180 radius=r0*math.tan(math.radians(deg))181 cone_weights[deg]+=float(wt[rr<=radius*radius].sum())182 max_dir_error=0.;weighted_path_exit=0.;exit_max_path=0.183 while len(ids) and iteration<max_scatter:184 iteration+=1185 _,exit_distance,ok=intervals(pos,u,geo)186 if not np.all(ok):187 raise RuntimeError("Interior ray lost cylinder intersection")188 free=rng.exponential(1/beta,len(ids))189 collide=free<exit_distance190 leaving=~collide191 if np.any(leaving):192 ww=w[leaving];endpos=pos[leaving]+exit_distance[leaving,None]*u[leaving]193 if iteration>1:cone_tally(endpos,u[leaving],ww)194 plen=path[leaving]+exit_distance[leaving]195 escaped_weight+=float(ww.sum());escaped_n+=int(leaving.sum())196 escaped_up+=float(ww[u[leaving,2]>0].sum())197 weighted_path_exit+=float(np.sum(ww*plen));exit_max_path=max(exit_max_path,float(plen.max()))198 top=abs(endpos[:,2]-geo.top_m)<1e-5199 base=(~top)&(abs(endpos[:,2]-geo.base_m)<1e-5)200 side=~(top|base)201 for label,mask in [("top",top),("base",base),("side",side)]:202 boundary_weights[label]+=float(ww[mask].sum())203 if not np.any(collide):204 ids=ids[:0];break205 ids=ids[collide];pos=pos[collide]+free[collide,None]*u[collide]206 path=path[collide]+free[collide];u=u[collide];w=w[collide]207 collisions+=len(ids)208 toward=-pos.copy();toward[:,2]+=geo.observer_z_m209 distance=np.linalg.norm(toward,axis=1);direction=toward/distance[:,None]210 cosine=np.einsum("ij,ij->i",u,direction)211 optical_length=segment_inside(pos,direction,distance,geo)212 aperture_cos=direction[:,2]213 accepted=aperture_cos>=math.cos(math.radians(geo.full_fov_deg/2))214 # Dimensionless relative to isotropic vacuum aperture fluence.215 contribution=w*4*math.pi*r0*r0*phase.pdf_sr(cosine)*aperture_cos/(distance*distance)*np.exp(-beta*optical_length)*accepted216 delay=(path+distance-r0)/C217 min_delay=min(min_delay,float(delay.min()))218 negative+=int(np.sum(delay < -1e-12))219 if negative:raise RuntimeError("Noncausal next-event contribution beyond roundoff")220 delay=np.maximum(delay,0)221 bins=np.minimum((delay/dt_s).astype(np.int64),nbins)222 part=np.bincount(bins,weights=contribution,minlength=nbins+1)223 hist+=part224 if iteration==1:first_hist+=part225 nee[ids]+=contribution226 mu=phase.sample_mu(rng,len(ids));phi=rng.uniform(0,2*math.pi,len(ids))227 u=rotate(u,mu,phi)228 max_dir_error=max(max_dir_error,float(np.max(abs(np.linalg.norm(u,axis=1)-1))))229 truncated_n=len(ids);truncated_weight=float(w.sum()) if truncated_n else 0.230 residual=launch_sum-escaped_weight-truncated_weight231 info={"seed":int(seed),"n_packets":n,"total_tau":float(total_tau),232 "overlying_tau":beta*direct_length,"source_altitude_m":source_m,233 "phase":phase.kind,"g":phase.g,"initial_importance_fraction":imp,"initial_importance_cone_deg":importance_cone_deg,234 "vacuum_transfer_m_inverse2":vacuum,235 "direct_relative_to_vacuum":direct,236 "scattered_relative_to_vacuum":float(hist.sum()/n),237 "total_relative_to_vacuum":direct+float(hist.sum()/n),238 "weighted_launch_mean":launch_sum/n,239 "weighted_escape_fraction":escaped_weight/n,"weighted_upward_escape_fraction":escaped_up/n,240 "weighted_truncation_fraction":truncated_weight/n,241 "energy_balance_residual_per_packet":residual/n,242 "escaped_packet_count":escaped_n,"truncated_packet_count":truncated_n,243 "boundary_escape_weights_per_packet":{k:v/n for k,v in boundary_weights.items()},244 "mean_collisions_per_packet":collisions/n,"max_scatter_iterations":iteration,245 "weighted_insideplus_entry_path_exit_mean_m":weighted_path_exit/n,246 "maximum_entryplus_cloud_path_m":exit_max_path,247 "min_scattered_delay_s":None if not np.isfinite(min_delay) else min_delay,248 "negative_delay_beyond_tolerance_count":negative,249 "max_direction_norm_error":max_dir_error,250 "response_overflow_relative_to_vacuum":float(hist[-1]/n),251 "first_scatter_relative_to_vacuum":float(first_hist.sum()/n),252 "finite_receiver_scattered_gain":{str(k):4*v/n/math.tan(math.radians(k))**2 for k,v in cone_weights.items()},253 "max_single_packet_scattered_gain":float(nee.max()),254 "fraction_scattered_gain_from_top10_packets":float(np.sort(nee)[-10:].sum()/nee.sum()) if nee.sum()>0 else 0.,255 "single_packet_standard_error_relative":float(np.std(nee,ddof=1)/math.sqrt(n)) if n>1 else 0.,256 "elapsed_s":time.perf_counter()-start}257 return info,hist/n258259def phase_qa(phase,seed=821013,n=300000):260 rng=np.random.default_rng(seed);mu=phase.sample_mu(rng,n)261 # Independent Gauss integration of piecewise-linear input.262 if phase.kind=="table":263 x,w=np.polynomial.legendre.leggauss(8)264 xx=(phase.mu[1:]+phase.mu[:-1])[:,None]/2+np.diff(phase.mu)[:,None]*x/2265 ww=np.diff(phase.mu)[:,None]*w/2266 f=phase.pdf_mu(xx)267 integrals=[float(np.sum(ww*f)),float(np.sum(ww*f*xx)),float(np.sum(ww*f*(3*xx*xx-1)/2))]268 else:269 integrals=[1.,phase.g,phase.g**2]270 sample=[float(mu.mean()),float(np.mean((3*mu*mu-1)/2))]271 se=[float(mu.std(ddof=1)/math.sqrt(n)),float(np.std((3*mu*mu-1)/2,ddof=1)/math.sqrt(n))]272 u=np.linspace(.000001,.999999,10000);inv=phase.inverse(u)273 if phase.kind=="table":274 k=np.minimum(np.searchsorted(phase.mu,inv,side="right")-1,len(phase.dx)-1);k=np.maximum(0,k)275 dx=inv-phase.mu[k]276 back=phase.cdf[k]+phase.f[k]*dx+.5*phase.slope[k]*dx*dx277 inverse_error=float(np.max(abs(back-u)))278 else:279 g=phase.g280 back=(1-g*g)/(2*g)*(1/np.sqrt(1+g*g-2*g*inv)-1/(1+g))281 inverse_error=float(np.max(abs(back-u)))282 return {"normalization":integrals[0],"g_integral":integrals[1],"g_parameter":phase.g,283 "second_Legendre_moment":integrals[2],"sample_size":n,"sample_moments":sample,284 "sample_standard_errors":se,"sample_z_scores":[(sample[i]-integrals[i+1])/se[i]for i in range(2)],285 "inverse_CDF_max_abs_error":inverse_error,286 "sample_in_bounds":bool(np.all((mu>=-1)&(mu<=1))),287 "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)))}288289if __name__=="__main__":290 ap=argparse.ArgumentParser();ap.add_argument("--lut",required=True);ap.add_argument("--out",required=True)291 ap.add_argument("--pilot",action="store_true");args=ap.parse_args()292 out=Path(args.out);out.mkdir(exist_ok=True)293 phases,meta=load_phases(args.lut)294 q={"phase_input":meta,"phases":{k:phase_qa(v)for k,v in phases.items()}}295 (out/"phase_QA.json").write_text(json.dumps(q,indent=2),encoding="utf-8")296 cases=[]297 for name in ["ghm","hg"]:298 for altitude,tau in [(15000.,1.),(15000.,45.),(17000.,10.),(18000.,10.)]:299 info,h=simulate(phases[name],altitude,tau,5000,810022+len(cases))300 cases.append(info);print(json.dumps(info),flush=True)301 (out/"pilot_results.json").write_text(json.dumps(cases,indent=2),encoding="utf-8")302