Review index · Download raw source

reference/analytic/auroral_counts.py

1from pathlib import Path2import csv,datetime,hashlib,importlib.util,json,math,sys3sys.path.insert(0,r'D:/Photometer_research_20261022/research/cloud_transport/private_readers')4from scipy.stats import poisson5from scipy.optimize import brentq67out=Path(r'D:/Photometer_research_20261022/revision_20261022_r4/analysis/analytic/auroral_counts')8src=Path(__file__).with_name('count_requirements.py')9spec=importlib.util.spec_from_file_location('read_only_count_methods',src)10m=importlib.util.module_from_spec(spec);spec.loader.exec_module(m)11hc=6.62607015e-34*299792458.;wave=337.1e-9;energy=hc/wave12eta=.2;Rph=1e10/(4*math.pi);Rw=Rph*energy;dark=50.;tau=17e-913rows=[];max_tail_error=0.;max_power_error=0.14for diameter in [1.,5.,25.]:15 area=math.pi*(diameter*.001)**2/416 for fov in [60.,90.]:17  projected=math.pi*math.sin(math.radians(fov/2))**218  G=area*projected19  for radiance_R in [15.,100.,1000.,2270.]:20   rate=radiance_R*Rph*G*eta21   for dt in [30e-6,80e-6,.0002,.003]:22    bg=(rate+dark)*dt23    for criterion,alpha in [('single_window_5sigma_tail',m.ALPHA_SINGLE),('illustrative_0p1_false_per_observing_day',m.ALPHA_OPERATIONAL)]:24     threshold,signal=m.count_at_power(bg,alpha,.95)25     sf=poisson.sf(threshold-1,bg)26     power=poisson.sf(threshold-1,bg+signal)27     assert sf<=alpha*(1+1e-7) and poisson.sf(threshold-2,bg)>alpha28     assert abs(power-.95)<1e-729     alt_signal=brentq(lambda s:poisson.sf(threshold-1,bg+s)-.95,0,signal*2+100)30     assert abs(alt_signal-signal)<1e-531     max_power_error=max(max_power_error,abs(power-.95));max_tail_error=max(max_tail_error,abs(sf-m.poisson_sf(threshold,bg)))32     rows.append(dict(diameter_mm=diameter,full_circular_FOV_deg=fov,projected_etendue_m2_sr=G,eta=eta,33      assumed_received_band_equivalent_R=radiance_R,nominal_four_nm_mean_W_m2_sr_nm=radiance_R*Rw/4,34      incident_linear_photoelectron_rate_s=rate,unmultiplied_mean_photocurrent_A=rate*1.602176634e-19,35      observation_window_s=dt,known_background_mean_count=bg,Poisson_sigma_count=math.sqrt(bg),36      aperture_photon_equivalent_rms=math.sqrt(bg)/eta,criterion=criterion,per_test_alpha=alpha,37      critical_observed_count=threshold,signal_mean_for_95pct_power=signal,aperture_photons_for_95pct_power=signal/eta,38      aperture_fluence_ph_m2_for_95pct_power=signal/(eta*area),39      rectangular_source_irradiance_W_m2_for_95pct_power=signal*energy/(eta*area*dt),40      hypothetical_17ns_deadtime_occupancy=rate*tau,41      hypothetical_nonparalyzable_fraction_lost=rate*tau/(1+rate*tau),42      model_validity='IdealindependentlinearPoissoncountsonly; count-losscolumnsarediagnosticsnotacorrectiontothelikelihood'))43with (out/'conditional_auroral_count_scenarios.csv').open('w',newline='',encoding='utf-8') as f:44 w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)45review={'recorded_utc':datetime.datetime.now(datetime.timezone.utc).isoformat(),46 'status':'EXECUTED_CONDITIONAL_DETECTOR_COUNTS_NOT_A_POLAR_BACKGROUND_FORECAST',47 'constants':{'photon_energy_J':energy,'one_R_ph_m2_sr_s':Rph,'one_R_W_m2_sr':Rw},48 'formulae':['G=A*pi*sin(full_circular_FOV/2)^2','r=Rband*1e10/(4*pi)*G*eta','B=(r+rdark)*dt','FindsmallestintegernwithP(Poisson(B)>=n)<=alpha; solveP(Poisson(B+S)>=n)=.95'],49 'conditions':{'uniform_received_band_radiance_R':[15,100,1000,2270],'end_to_end_efficiency':eta,'dark_cps':dark,'nominal_rectangular_width_nm':4,'arrival_counts_before_deadtime':True},50 'relation_to_observations':'2270RusesSivjee1976selectedN2(0,0)bandintensityonlyasanamplitudescale. Complete4nmcapture, uniformreceivedradianceandsharedgeometryareassumptions, notobservedbythatpaper. No1971photometerbrightnessistreatedasactual4nmR.',51 'no_silent_transfers':['NoKpclassoroccurrenceassigned','Noextrapolationfrom10degreeaircraftFOVtofull60/90degreearc','Nocloudreflectionorlimbtonadirtransfercomputed','Noactualsourcephotonfluenceorpopulationinferred','Noactualfilter/PDE/temperature/radiationdrift','NoPMTmultiplicationexcessnoise,SiPMcorrelatednoise,particles,baselineuncertaintyordeadtimeinPoissonlikelihood','Single-window5sigmatailisnotglobalfalse-triggercontrol; operationalallocationishypothetical50kstarts/s*86400s*3families'],52 'validation':{'scipy_independent_tail_and_root_checks':len(rows),'max_absolute_tail_difference':max_tail_error,'max_absolute_power_error':max_power_error,'read_only_original_method_sha256':hashlib.sha256(src.read_bytes()).hexdigest()},53 'examples':[r for r in rows if r['diameter_mm']==25 and r['full_circular_FOV_deg']==60 and r['observation_window_s']==30e-6 and r['criterion']=='single_window_5sigma_tail']}54(out/'conditional_auroral_count_scenarios.json').write_text(json.dumps(review,indent=2),encoding='utf-8')55print(json.dumps({'rows':len(rows),'constants':review['constants'],'examples':review['examples'],'validation':review['validation']},indent=2))56