From d2a901cab11c88a43068746afabba690e7c8a7ff Mon Sep 17 00:00:00 2001 From: Yingjie Wang Date: Sat, 12 Sep 2026 15:29:33 -0400 Subject: [PATCH] Feat: Add CSV CIE 1931 blackbody LUT backend Replace the in-tree Wyman analytic CIE fit with a repository CIE 1931 2-degree LUT derived from the 360-830 nm, 1 nm-linear CSV reference. The GRBBLUT3 table stores 1024 uniform log(T) XYZ nodes over [670.146556, 101408.88] K with four-point cubic Lagrange interpolation, and the same CSV reference supplies a three-term Rayleigh-Jeans form above the table and an inverse-temperature/log-XYZ crossover with channel-specific endpoint forms below it. The loader validates the header and payload checksum and rejects legacy formats and runtime generation. Initialize the backend in the renderer, frame test, and PSF capture tool so the new default is active wherever colors are produced. --- assets/CIE_xyz_1931_2deg.csv | 471 ++++++++++++++++++ assets/blackbody/README.md | 40 ++ .../blackbody/cie1931_2deg_xyz_1024.grbblut | Bin 0 -> 24736 bytes src/blackbody_common.c | 128 +++++ src/blackbody_internal.h | 12 + src/blackbody_lut.c | 175 +++++++ src/main.c | 23 + src/optics.c | 62 --- src/optics.h | 10 + tests/capture_psf.c | 6 +- tests/test_frame.c | 4 +- 11 files changed, 866 insertions(+), 65 deletions(-) create mode 100644 assets/CIE_xyz_1931_2deg.csv create mode 100644 assets/blackbody/README.md create mode 100644 assets/blackbody/cie1931_2deg_xyz_1024.grbblut create mode 100644 src/blackbody_common.c create mode 100644 src/blackbody_internal.h create mode 100644 src/blackbody_lut.c diff --git a/assets/CIE_xyz_1931_2deg.csv b/assets/CIE_xyz_1931_2deg.csv new file mode 100644 index 0000000..9b4e3f7 --- /dev/null +++ b/assets/CIE_xyz_1931_2deg.csv @@ -0,0 +1,471 @@ +360,0.000129900000,0.0000039170000,0.000606100000 +361,0.000145847000,0.0000043935810,0.000680879200 +362,0.000163802100,0.0000049296040,0.000765145600 +363,0.000184003700,0.0000055321360,0.000860012400 +364,0.000206690200,0.0000062082450,0.000966592800 +365,0.000232100000,0.0000069650000,0.001086000000 +366,0.000260728000,0.0000078132190,0.001220586000 +367,0.000293075000,0.0000087673360,0.001372729000 +368,0.000329388000,0.0000098398440,0.001543579000 +369,0.000369914000,0.0000110432300,0.001734286000 +370,0.000414900000,0.0000123900000,0.001946000000 +371,0.000464158700,0.0000138864100,0.002177777000 +372,0.000518986000,0.0000155572800,0.002435809000 +373,0.000581854000,0.0000174429600,0.002731953000 +374,0.000655234700,0.0000195837500,0.003078064000 +375,0.000741600000,0.0000220200000,0.003486000000 +376,0.000845029600,0.0000248396500,0.003975227000 +377,0.000964526800,0.0000280412600,0.004540880000 +378,0.001094949000,0.0000315310400,0.005158320000 +379,0.001231154000,0.0000352152100,0.005802907000 +380,0.001368000000,0.0000390000000,0.006450001000 +381,0.001502050000,0.0000428264000,0.007083216000 +382,0.001642328000,0.0000469146000,0.007745488000 +383,0.001802382000,0.0000515896000,0.008501152000 +384,0.001995757000,0.0000571764000,0.009414544000 +385,0.002236000000,0.0000640000000,0.010549990000 +386,0.002535385000,0.0000723442100,0.011965800000 +387,0.002892603000,0.0000822122400,0.013655870000 +388,0.003300829000,0.0000935081600,0.015588050000 +389,0.003753236000,0.0001061361000,0.017730150000 +390,0.004243000000,0.0001200000000,0.020050010000 +391,0.004762389000,0.0001349840000,0.022511360000 +392,0.005330048000,0.0001514920000,0.025202880000 +393,0.005978712000,0.0001702080000,0.028279720000 +394,0.006741117000,0.0001918160000,0.031897040000 +395,0.007650000000,0.0002170000000,0.036210000000 +396,0.008751373000,0.0002469067000,0.041437710000 +397,0.010028880000,0.0002812400000,0.047503720000 +398,0.011421700000,0.0003185200000,0.054119880000 +399,0.012869010000,0.0003572667000,0.060998030000 +400,0.014310000000,0.0003960000000,0.067850010000 +401,0.015704430000,0.0004337147000,0.074486320000 +402,0.017147440000,0.0004730240000,0.081361560000 +403,0.018781220000,0.0005178760000,0.089153640000 +404,0.020748010000,0.0005722187000,0.098540480000 +405,0.023190000000,0.0006400000000,0.110200000000 +406,0.026207360000,0.0007245600000,0.124613300000 +407,0.029782480000,0.0008255000000,0.141701700000 +408,0.033880920000,0.0009411600000,0.161303500000 +409,0.038468240000,0.0010698800000,0.183256800000 +410,0.043510000000,0.0012100000000,0.207400000000 +411,0.048995600000,0.0013620910000,0.233692100000 +412,0.055022600000,0.0015307520000,0.262611400000 +413,0.061718800000,0.0017203680000,0.294774600000 +414,0.069212000000,0.0019353230000,0.330798500000 +415,0.077630000000,0.0021800000000,0.371300000000 +416,0.086958110000,0.0024548000000,0.416209100000 +417,0.097176720000,0.0027640000000,0.465464200000 +418,0.108406300000,0.0031178000000,0.519694800000 +419,0.120767200000,0.0035264000000,0.579530300000 +420,0.134380000000,0.0040000000000,0.645600000000 +421,0.149358200000,0.0045462400000,0.718483800000 +422,0.165395700000,0.0051593200000,0.796713300000 +423,0.181983100000,0.0058292800000,0.877845900000 +424,0.198611000000,0.0065461600000,0.959439000000 +425,0.214770000000,0.0073000000000,1.039050100000 +426,0.230186800000,0.0080865070000,1.115367300000 +427,0.244879700000,0.0089087200000,1.188497100000 +428,0.258777300000,0.0097676800000,1.258123300000 +429,0.271807900000,0.0106644300000,1.323929600000 +430,0.283900000000,0.0116000000000,1.385600000000 +431,0.294943800000,0.0125731700000,1.442635200000 +432,0.304896500000,0.0135827200000,1.494803500000 +433,0.313787300000,0.0146296800000,1.542190300000 +434,0.321645400000,0.0157150900000,1.584880700000 +435,0.328500000000,0.0168400000000,1.622960000000 +436,0.334351300000,0.0180073600000,1.656404800000 +437,0.339210100000,0.0192144800000,1.685295900000 +438,0.343121300000,0.0204539200000,1.709874500000 +439,0.346129600000,0.0217182400000,1.730382100000 +440,0.348280000000,0.0230000000000,1.747060000000 +441,0.349599900000,0.0242946100000,1.760044600000 +442,0.350147400000,0.0256102400000,1.769623300000 +443,0.350013000000,0.0269585700000,1.776263700000 +444,0.349287000000,0.0283512500000,1.780433400000 +445,0.348060000000,0.0298000000000,1.782600000000 +446,0.346373300000,0.0313108300000,1.782968200000 +447,0.344262400000,0.0328836800000,1.781699800000 +448,0.341808800000,0.0345211200000,1.779198200000 +449,0.339094100000,0.0362257100000,1.775867100000 +450,0.336200000000,0.0380000000000,1.772110000000 +451,0.333197700000,0.0398466700000,1.768258900000 +452,0.330041100000,0.0417680000000,1.764039000000 +453,0.326635700000,0.0437660000000,1.758943800000 +454,0.322886800000,0.0458426700000,1.752466300000 +455,0.318700000000,0.0480000000000,1.744100000000 +456,0.314025100000,0.0502436800000,1.733559500000 +457,0.308884000000,0.0525730400000,1.720858100000 +458,0.303290400000,0.0549805600000,1.705936900000 +459,0.297257900000,0.0574587200000,1.688737200000 +460,0.290800000000,0.0600000000000,1.669200000000 +461,0.283970100000,0.0626019700000,1.647528700000 +462,0.276721400000,0.0652775200000,1.623412700000 +463,0.268917800000,0.0680420800000,1.596022300000 +464,0.260422700000,0.0709110900000,1.564528000000 +465,0.251100000000,0.0739000000000,1.528100000000 +466,0.240847500000,0.0770160000000,1.486111400000 +467,0.229851200000,0.0802664000000,1.439521500000 +468,0.218407200000,0.0836668000000,1.389879900000 +469,0.206811500000,0.0872328000000,1.338736200000 +470,0.195360000000,0.0909800000000,1.287640000000 +471,0.184213600000,0.0949175500000,1.237422300000 +472,0.173327300000,0.0990458400000,1.187824300000 +473,0.162688100000,0.1033674000000,1.138761100000 +474,0.152283300000,0.1078846000000,1.090148000000 +475,0.142100000000,0.1126000000000,1.041900000000 +476,0.132178600000,0.1175320000000,0.994197600000 +477,0.122569600000,0.1226744000000,0.947347300000 +478,0.113275200000,0.1279928000000,0.901453100000 +479,0.104297900000,0.1334528000000,0.856619300000 +480,0.095640000000,0.1390200000000,0.812950100000 +481,0.087299550000,0.1446764000000,0.770517300000 +482,0.079308040000,0.1504693000000,0.729444800000 +483,0.071717760000,0.1564619000000,0.689913600000 +484,0.064580990000,0.1627177000000,0.652104900000 +485,0.057950010000,0.1693000000000,0.616200000000 +486,0.051862110000,0.1762431000000,0.582328600000 +487,0.046281520000,0.1835581000000,0.550416200000 +488,0.041150880000,0.1912735000000,0.520337600000 +489,0.036412830000,0.1994180000000,0.491967300000 +490,0.032010000000,0.2080200000000,0.465180000000 +491,0.027917200000,0.2171199000000,0.439924600000 +492,0.024144400000,0.2267345000000,0.416183600000 +493,0.020687000000,0.2368571000000,0.393882200000 +494,0.017540400000,0.2474812000000,0.372945900000 +495,0.014700000000,0.2586000000000,0.353300000000 +496,0.012161790000,0.2701849000000,0.334857800000 +497,0.009919960000,0.2822939000000,0.317552100000 +498,0.007967240000,0.2950505000000,0.301337500000 +499,0.006296346000,0.3085780000000,0.286168600000 +500,0.004900000000,0.3230000000000,0.272000000000 +501,0.003777173000,0.3384021000000,0.258817100000 +502,0.002945320000,0.3546858000000,0.246483800000 +503,0.002424880000,0.3716986000000,0.234771800000 +504,0.002236293000,0.3892875000000,0.223453300000 +505,0.002400000000,0.4073000000000,0.212300000000 +506,0.002925520000,0.4256299000000,0.201169200000 +507,0.003836560000,0.4443096000000,0.190119600000 +508,0.005174840000,0.4633944000000,0.179225400000 +509,0.006982080000,0.4829395000000,0.168560800000 +510,0.009300000000,0.5030000000000,0.158200000000 +511,0.012149490000,0.5235693000000,0.148138300000 +512,0.015535880000,0.5445120000000,0.138375800000 +513,0.019477520000,0.5656900000000,0.128994200000 +514,0.023992770000,0.5869653000000,0.120075100000 +515,0.029100000000,0.6082000000000,0.111700000000 +516,0.034814850000,0.6293456000000,0.103904800000 +517,0.041120160000,0.6503068000000,0.096667480000 +518,0.047985040000,0.6708752000000,0.089982720000 +519,0.055378610000,0.6908424000000,0.083845310000 +520,0.063270000000,0.7100000000000,0.078249990000 +521,0.071635010000,0.7281852000000,0.073208990000 +522,0.080462240000,0.7454636000000,0.068678160000 +523,0.089739960000,0.7619694000000,0.064567840000 +524,0.099456450000,0.7778368000000,0.060788350000 +525,0.109600000000,0.7932000000000,0.057250010000 +526,0.120167400000,0.8081104000000,0.053904350000 +527,0.131114500000,0.8224962000000,0.050746640000 +528,0.142367900000,0.8363068000000,0.047752760000 +529,0.153854200000,0.8494916000000,0.044898590000 +530,0.165500000000,0.8620000000000,0.042160000000 +531,0.177257100000,0.8738108000000,0.039507280000 +532,0.189140000000,0.8849624000000,0.036935640000 +533,0.201169400000,0.8954936000000,0.034458360000 +534,0.213365800000,0.9054432000000,0.032088720000 +535,0.225749900000,0.9148501000000,0.029840000000 +536,0.238320900000,0.9237348000000,0.027711810000 +537,0.251066800000,0.9320924000000,0.025694440000 +538,0.263992200000,0.9399226000000,0.023787160000 +539,0.277101700000,0.9472252000000,0.021989250000 +540,0.290400000000,0.9540000000000,0.020300000000 +541,0.303891200000,0.9602561000000,0.018718050000 +542,0.317572600000,0.9660074000000,0.017240360000 +543,0.331438400000,0.9712606000000,0.015863640000 +544,0.345482800000,0.9760225000000,0.014584610000 +545,0.359700000000,0.9803000000000,0.013400000000 +546,0.374083900000,0.9840924000000,0.012307230000 +547,0.388639600000,0.9874182000000,0.011301880000 +548,0.403378400000,0.9903128000000,0.010377920000 +549,0.418311500000,0.9928116000000,0.009529306000 +550,0.433449900000,0.9949501000000,0.008749999000 +551,0.448795300000,0.9967108000000,0.008035200000 +552,0.464336000000,0.9980983000000,0.007381600000 +553,0.480064000000,0.9991120000000,0.006785400000 +554,0.495971300000,0.9997482000000,0.006242800000 +555,0.512050100000,1.0000000000000,0.005749999000 +556,0.528295900000,0.9998567000000,0.005303600000 +557,0.544691600000,0.9993046000000,0.004899800000 +558,0.561209400000,0.9983255000000,0.004534200000 +559,0.577821500000,0.9968987000000,0.004202400000 +560,0.594500000000,0.9950000000000,0.003900000000 +561,0.611220900000,0.9926005000000,0.003623200000 +562,0.627975800000,0.9897426000000,0.003370600000 +563,0.644760200000,0.9864444000000,0.003141400000 +564,0.661569700000,0.9827241000000,0.002934800000 +565,0.678400000000,0.9786000000000,0.002749999000 +566,0.695239200000,0.9740837000000,0.002585200000 +567,0.712058600000,0.9691712000000,0.002438600000 +568,0.728828400000,0.9638568000000,0.002309400000 +569,0.745518800000,0.9581349000000,0.002196800000 +570,0.762100000000,0.9520000000000,0.002100000000 +571,0.778543200000,0.9454504000000,0.002017733000 +572,0.794825600000,0.9384992000000,0.001948200000 +573,0.810926400000,0.9311628000000,0.001889800000 +574,0.826824800000,0.9234576000000,0.001840933000 +575,0.842500000000,0.9154000000000,0.001800000000 +576,0.857932500000,0.9070064000000,0.001766267000 +577,0.873081600000,0.8982772000000,0.001737800000 +578,0.887894400000,0.8892048000000,0.001711200000 +579,0.902318100000,0.8797816000000,0.001683067000 +580,0.916300000000,0.8700000000000,0.001650001000 +581,0.929799500000,0.8598613000000,0.001610133000 +582,0.942798400000,0.8493920000000,0.001564400000 +583,0.955277600000,0.8386220000000,0.001513600000 +584,0.967217900000,0.8275813000000,0.001458533000 +585,0.978600000000,0.8163000000000,0.001400000000 +586,0.989385600000,0.8047947000000,0.001336667000 +587,0.999548800000,0.7930820000000,0.001270000000 +588,1.009089200000,0.7811920000000,0.001205000000 +589,1.018006400000,0.7691547000000,0.001146667000 +590,1.026300000000,0.7570000000000,0.001100000000 +591,1.033982700000,0.7447541000000,0.001068800000 +592,1.040986000000,0.7324224000000,0.001049400000 +593,1.047188000000,0.7200036000000,0.001035600000 +594,1.052466700000,0.7074965000000,0.001021200000 +595,1.056700000000,0.6949000000000,0.001000000000 +596,1.059794400000,0.6822192000000,0.000968640000 +597,1.061799200000,0.6694716000000,0.000929920000 +598,1.062806800000,0.6566744000000,0.000886880000 +599,1.062909600000,0.6438448000000,0.000842560000 +600,1.062200000000,0.6310000000000,0.000800000000 +601,1.060735200000,0.6181555000000,0.000760960000 +602,1.058443600000,0.6053144000000,0.000723680000 +603,1.055224400000,0.5924756000000,0.000685920000 +604,1.050976800000,0.5796379000000,0.000645440000 +605,1.045600000000,0.5668000000000,0.000600000000 +606,1.039036900000,0.5539611000000,0.000547866700 +607,1.031360800000,0.5411372000000,0.000491600000 +608,1.022666200000,0.5283528000000,0.000435400000 +609,1.013047700000,0.5156323000000,0.000383466700 +610,1.002600000000,0.5030000000000,0.000340000000 +611,0.991367500000,0.4904688000000,0.000307253300 +612,0.979331400000,0.4780304000000,0.000283160000 +613,0.966491600000,0.4656776000000,0.000265440000 +614,0.952847900000,0.4534032000000,0.000251813300 +615,0.938400000000,0.4412000000000,0.000240000000 +616,0.923194000000,0.4290800000000,0.000229546700 +617,0.907244000000,0.4170360000000,0.000220640000 +618,0.890502000000,0.4050320000000,0.000211960000 +619,0.872920000000,0.3930320000000,0.000202186700 +620,0.854449900000,0.3810000000000,0.000190000000 +621,0.835084000000,0.3689184000000,0.000174213300 +622,0.814946000000,0.3568272000000,0.000155640000 +623,0.794186000000,0.3447768000000,0.000135960000 +624,0.772954000000,0.3328176000000,0.000116853300 +625,0.751400000000,0.3210000000000,0.000100000000 +626,0.729583600000,0.3093381000000,0.000086133330 +627,0.707588800000,0.2978504000000,0.000074600000 +628,0.685602200000,0.2865936000000,0.000065000000 +629,0.663810400000,0.2756245000000,0.000056933330 +630,0.642400000000,0.2650000000000,0.000049999990 +631,0.621514900000,0.2547632000000,0.000044160000 +632,0.601113800000,0.2448896000000,0.000039480000 +633,0.581105200000,0.2353344000000,0.000035720000 +634,0.561397700000,0.2260528000000,0.000032640000 +635,0.541900000000,0.2170000000000,0.000030000000 +636,0.522599500000,0.2081616000000,0.000027653330 +637,0.503546400000,0.1995488000000,0.000025560000 +638,0.484743600000,0.1911552000000,0.000023640000 +639,0.466193900000,0.1829744000000,0.000021813330 +640,0.447900000000,0.1750000000000,0.000020000000 +641,0.429861300000,0.1672235000000,0.000018133330 +642,0.412098000000,0.1596464000000,0.000016200000 +643,0.394644000000,0.1522776000000,0.000014200000 +644,0.377533300000,0.1451259000000,0.000012133330 +645,0.360800000000,0.1382000000000,0.000010000000 +646,0.344456300000,0.1315003000000,0.000007733333 +647,0.328516800000,0.1250248000000,0.000005400000 +648,0.313019200000,0.1187792000000,0.000003200000 +649,0.298001100000,0.1127691000000,0.000001333333 +650,0.283500000000,0.1070000000000,0.000000000000 +651,0.269544800000,0.1014762000000,0.000000000000 +652,0.256118400000,0.0961886400000,0.000000000000 +653,0.243189600000,0.0911229600000,0.000000000000 +654,0.230727200000,0.0862648500000,0.000000000000 +655,0.218700000000,0.0816000000000,0.000000000000 +656,0.207097100000,0.0771206400000,0.000000000000 +657,0.195923200000,0.0728255200000,0.000000000000 +658,0.185170800000,0.0687100800000,0.000000000000 +659,0.174832300000,0.0647697600000,0.000000000000 +660,0.164900000000,0.0610000000000,0.000000000000 +661,0.155366700000,0.0573962100000,0.000000000000 +662,0.146230000000,0.0539550400000,0.000000000000 +663,0.137490000000,0.0506737600000,0.000000000000 +664,0.129146700000,0.0475496500000,0.000000000000 +665,0.121200000000,0.0445800000000,0.000000000000 +666,0.113639700000,0.0417587200000,0.000000000000 +667,0.106465000000,0.0390849600000,0.000000000000 +668,0.099690440000,0.0365638400000,0.000000000000 +669,0.093330610000,0.0342004800000,0.000000000000 +670,0.087400000000,0.0320000000000,0.000000000000 +671,0.081900960000,0.0299626100000,0.000000000000 +672,0.076804280000,0.0280766400000,0.000000000000 +673,0.072077120000,0.0263293600000,0.000000000000 +674,0.067686640000,0.0247080500000,0.000000000000 +675,0.063600000000,0.0232000000000,0.000000000000 +676,0.059806850000,0.0218007700000,0.000000000000 +677,0.056282160000,0.0205011200000,0.000000000000 +678,0.052971040000,0.0192810800000,0.000000000000 +679,0.049818610000,0.0181206900000,0.000000000000 +680,0.046770000000,0.0170000000000,0.000000000000 +681,0.043784050000,0.0159037900000,0.000000000000 +682,0.040875360000,0.0148371800000,0.000000000000 +683,0.038072640000,0.0138106800000,0.000000000000 +684,0.035404610000,0.0128347800000,0.000000000000 +685,0.032900000000,0.0119200000000,0.000000000000 +686,0.030564190000,0.0110683100000,0.000000000000 +687,0.028380560000,0.0102733900000,0.000000000000 +688,0.026344840000,0.0095333110000,0.000000000000 +689,0.024452750000,0.0088461570000,0.000000000000 +690,0.022700000000,0.0082100000000,0.000000000000 +691,0.021084290000,0.0076237810000,0.000000000000 +692,0.019599880000,0.0070854240000,0.000000000000 +693,0.018237320000,0.0065914760000,0.000000000000 +694,0.016987170000,0.0061384850000,0.000000000000 +695,0.015840000000,0.0057230000000,0.000000000000 +696,0.014790640000,0.0053430590000,0.000000000000 +697,0.013831320000,0.0049957960000,0.000000000000 +698,0.012948680000,0.0046764040000,0.000000000000 +699,0.012129200000,0.0043800750000,0.000000000000 +700,0.011359160000,0.0041020000000,0.000000000000 +701,0.010629350000,0.0038384530000,0.000000000000 +702,0.009938846000,0.0035890990000,0.000000000000 +703,0.009288422000,0.0033542190000,0.000000000000 +704,0.008678854000,0.0031340930000,0.000000000000 +705,0.008110916000,0.0029290000000,0.000000000000 +706,0.007582388000,0.0027381390000,0.000000000000 +707,0.007088746000,0.0025598760000,0.000000000000 +708,0.006627313000,0.0023932440000,0.000000000000 +709,0.006195408000,0.0022372750000,0.000000000000 +710,0.005790346000,0.0020910000000,0.000000000000 +711,0.005409826000,0.0019535870000,0.000000000000 +712,0.005052583000,0.0018245800000,0.000000000000 +713,0.004717512000,0.0017035800000,0.000000000000 +714,0.004403507000,0.0015901870000,0.000000000000 +715,0.004109457000,0.0014840000000,0.000000000000 +716,0.003833913000,0.0013844960000,0.000000000000 +717,0.003575748000,0.0012912680000,0.000000000000 +718,0.003334342000,0.0012040920000,0.000000000000 +719,0.003109075000,0.0011227440000,0.000000000000 +720,0.002899327000,0.0010470000000,0.000000000000 +721,0.002704348000,0.0009765896000,0.000000000000 +722,0.002523020000,0.0009111088000,0.000000000000 +723,0.002354168000,0.0008501332000,0.000000000000 +724,0.002196616000,0.0007932384000,0.000000000000 +725,0.002049190000,0.0007400000000,0.000000000000 +726,0.001910960000,0.0006900827000,0.000000000000 +727,0.001781438000,0.0006433100000,0.000000000000 +728,0.001660110000,0.0005994960000,0.000000000000 +729,0.001546459000,0.0005584547000,0.000000000000 +730,0.001439971000,0.0005200000000,0.000000000000 +731,0.001340042000,0.0004839136000,0.000000000000 +732,0.001246275000,0.0004500528000,0.000000000000 +733,0.001158471000,0.0004183452000,0.000000000000 +734,0.001076430000,0.0003887184000,0.000000000000 +735,0.000999949300,0.0003611000000,0.000000000000 +736,0.000928735800,0.0003353835000,0.000000000000 +737,0.000862433200,0.0003114404000,0.000000000000 +738,0.000800750300,0.0002891656000,0.000000000000 +739,0.000743396000,0.0002684539000,0.000000000000 +740,0.000690078600,0.0002492000000,0.000000000000 +741,0.000640515600,0.0002313019000,0.000000000000 +742,0.000594502100,0.0002146856000,0.000000000000 +743,0.000551864600,0.0001992884000,0.000000000000 +744,0.000512429000,0.0001850475000,0.000000000000 +745,0.000476021300,0.0001719000000,0.000000000000 +746,0.000442453600,0.0001597781000,0.000000000000 +747,0.000411511700,0.0001486044000,0.000000000000 +748,0.000382981400,0.0001383016000,0.000000000000 +749,0.000356649100,0.0001287925000,0.000000000000 +750,0.000332301100,0.0001200000000,0.000000000000 +751,0.000309758600,0.0001118595000,0.000000000000 +752,0.000288887100,0.0001043224000,0.000000000000 +753,0.000269539400,0.0000973356000,0.000000000000 +754,0.000251568200,0.0000908458700,0.000000000000 +755,0.000234826100,0.0000848000000,0.000000000000 +756,0.000219171000,0.0000791466700,0.000000000000 +757,0.000204525800,0.0000738580000,0.000000000000 +758,0.000190840500,0.0000689160000,0.000000000000 +759,0.000178065400,0.0000643026700,0.000000000000 +760,0.000166150500,0.0000600000000,0.000000000000 +761,0.000155023600,0.0000559818700,0.000000000000 +762,0.000144621900,0.0000522256000,0.000000000000 +763,0.000134909800,0.0000487184000,0.000000000000 +764,0.000125852000,0.0000454474700,0.000000000000 +765,0.000117413000,0.0000424000000,0.000000000000 +766,0.000109551500,0.0000395610400,0.000000000000 +767,0.000102224500,0.0000369151200,0.000000000000 +768,0.000095394450,0.0000344486800,0.000000000000 +769,0.000089023900,0.0000321481600,0.000000000000 +770,0.000083075270,0.0000300000000,0.000000000000 +771,0.000077512690,0.0000279912500,0.000000000000 +772,0.000072313040,0.0000261135600,0.000000000000 +773,0.000067457780,0.0000243602400,0.000000000000 +774,0.000062928440,0.0000227246100,0.000000000000 +775,0.000058706520,0.0000212000000,0.000000000000 +776,0.000054770280,0.0000197785500,0.000000000000 +777,0.000051099180,0.0000184528500,0.000000000000 +778,0.000047676540,0.0000172168700,0.000000000000 +779,0.000044485670,0.0000160645900,0.000000000000 +780,0.000041509940,0.0000149900000,0.000000000000 +781,0.000038733240,0.0000139872800,0.000000000000 +782,0.000036142030,0.0000130515500,0.000000000000 +783,0.000033723520,0.0000121781800,0.000000000000 +784,0.000031464870,0.0000113625400,0.000000000000 +785,0.000029353260,0.0000106000000,0.000000000000 +786,0.000027375730,0.0000098858770,0.000000000000 +787,0.000025524330,0.0000092173040,0.000000000000 +788,0.000023793760,0.0000085923620,0.000000000000 +789,0.000022178700,0.0000080091330,0.000000000000 +790,0.000020673830,0.0000074657000,0.000000000000 +791,0.000019272260,0.0000069595670,0.000000000000 +792,0.000017966400,0.0000064879950,0.000000000000 +793,0.000016749910,0.0000060486990,0.000000000000 +794,0.000015616480,0.0000056393960,0.000000000000 +795,0.000014559770,0.0000052578000,0.000000000000 +796,0.000013573870,0.0000049017710,0.000000000000 +797,0.000012654360,0.0000045697200,0.000000000000 +798,0.000011797230,0.0000042601940,0.000000000000 +799,0.000010998440,0.0000039717390,0.000000000000 +800,0.000010253980,0.0000037029000,0.000000000000 +801,0.000009559646,0.0000034521630,0.000000000000 +802,0.000008912044,0.0000032183020,0.000000000000 +803,0.000008308358,0.0000030003000,0.000000000000 +804,0.000007745769,0.0000027971390,0.000000000000 +805,0.000007221456,0.0000026078000,0.000000000000 +806,0.000006732475,0.0000024312200,0.000000000000 +807,0.000006276423,0.0000022665310,0.000000000000 +808,0.000005851304,0.0000021130130,0.000000000000 +809,0.000005455118,0.0000019699430,0.000000000000 +810,0.000005085868,0.0000018366000,0.000000000000 +811,0.000004741466,0.0000017122300,0.000000000000 +812,0.000004420236,0.0000015962280,0.000000000000 +813,0.000004120783,0.0000014880900,0.000000000000 +814,0.000003841716,0.0000013873140,0.000000000000 +815,0.000003581652,0.0000012934000,0.000000000000 +816,0.000003339127,0.0000012058200,0.000000000000 +817,0.000003112949,0.0000011241430,0.000000000000 +818,0.000002902121,0.0000010480090,0.000000000000 +819,0.000002705645,0.0000009770578,0.000000000000 +820,0.000002522525,0.0000009109300,0.000000000000 +821,0.000002351726,0.0000008492513,0.000000000000 +822,0.000002192415,0.0000007917212,0.000000000000 +823,0.000002043902,0.0000007380904,0.000000000000 +824,0.000001905497,0.0000006881098,0.000000000000 +825,0.000001776509,0.0000006415300,0.000000000000 +826,0.000001656215,0.0000005980895,0.000000000000 +827,0.000001544022,0.0000005575746,0.000000000000 +828,0.000001439440,0.0000005198080,0.000000000000 +829,0.000001341977,0.0000004846123,0.000000000000 +830,0.000001251141,0.0000004518100,0.000000000000 diff --git a/assets/blackbody/README.md b/assets/blackbody/README.md new file mode 100644 index 0000000..531a28b --- /dev/null +++ b/assets/blackbody/README.md @@ -0,0 +1,40 @@ +# CIE 1931 blackbody LUT + +`cie1931_2deg_xyz_1024.grbblut` is the sole repository LUT. It contains 1,024 +XYZ nodes uniformly spaced in `log(T)` over `[670.146556, 101408.88] K`. +Wolfram integrates the Planck spectrum against the piecewise-linear 1 nm data +in `assets/CIE_xyz_1931_2deg.csv` at 40-digit working precision and rounds each +stored value once to little-endian IEEE binary64. + +The 160-byte `GRBBLUT3` header records the CIE 1931 2-degree, 360--830 nm, +1 nm-linear reference definition, node count, temperature limits, physical +constants, generator description, and payload FNV-1a checksum. The C loader +validates the header, exact payload size, and checksum, and rejects trailing +bytes. Runtime table generation and the former `GRBBLUT1`/`GRBBLUT2` formats are +not supported. The evaluator uses four-point cubic Lagrange interpolation in +`log(T)`. + +- file size: 24,736 bytes +- payload FNV-1a: `0e39867b0e70a809` +- file SHA-256: + `850446663212fe488c88fca3ac2f54d6f8119a04801dda84fc1a62acdec37d63` +- all 1,023 interval midpoints versus independent Gauss-16 integration: + max XYZ-channel relative `5.4633e-5`, no negative XYZ + +The LUT interpolation acceptance limit is `1e-4`. A separate `5e-4` overall +relative error budget accounts for the estimated error of linearly +interpolating the 1 nm CIE samples; the LUT's much smaller internal error must +not be mistaken for physical-reference accuracy. + +Outside the stored range, the same CSV reference supplies both models. Above +101408.88 K the evaluator uses three integrated Rayleigh--Jeans terms +`A T+B+C/T` (maximum validated XYZ-channel relative error `1.6635e-5` +through `10*Tmax`). Below 670.146556 K it uses a six-segment, degree-eight +inverse-temperature/log-XYZ Chebyshev--Lobatto crossover down to 4.097321 K +and channel-specific endpoint forms below that. Their maximum sampled XYZ +relative errors are `8.16e-5` and `6.69e-5`. + +The Wolfram generation/validation scripts that produced this table, together +with the bounded performance and accuracy evidence, are preserved on the +experiment branch `codex/blackbody-cost-experiment` (tag +`blackbody-cost-experiment-2026-09-12`) and are not part of the mainline tree. diff --git a/assets/blackbody/cie1931_2deg_xyz_1024.grbblut b/assets/blackbody/cie1931_2deg_xyz_1024.grbblut new file mode 100644 index 0000000000000000000000000000000000000000..7d45f2b7041058110b4d69065d0deea03592116b GIT binary patch literal 24736 zcmW(+byQHz(Qq@_jamIeXo?ha`Y5fCI4415U*6{JK_8bLr% zLM0XW{C$7i-JQGV?4I47=b4?EvugTEN?JyS;=}|51f~DJW5h&+g#YUQ$&UXT*uTdA zWP)S=D+maz!-lzQNzWidJEkAK)VzpZlxUhVhEn!lmqx&J0|eV%`>ZS558r6m4; zl79vH0}q;@7+D&jEk?mlmqM_b6Z>2=zn!o@pWbnHTv9y9P3@-7N#un&OSn`P84M{Z-t*JXN{+)O1uiWMi27n-RS(XT zx4%~`J5PT$WiSbw(XAw#T@ApLhxz$y=Ef9bJ2Mf&TA@$0rrn zKm4rx{WcTJvQVSFEEbNDnO$wU>oTXb9w>ne24+ZUr{krZ>+_l! z_0B|=7<<0+gzUy%6SkvZdxE364C5R;@LFIX#bgZFPTj#;us6X~Z@M4c!;%RO%pN|X z#1e)@L@}}#*rZGZ;~L`wOyC!P@lYNeCbQH2E7-XWEB5H*521dDC2P-qfBS_2tBCzl zrigW7M34BCP8rr?ZJaBpZSWMfs?$GuO`sdA;Dh7)hmF{|jY4Oi5LS#if6RX<#A>4MNa`?r#E!RCRIJ-ZUKkxw{+1>`t)h`9hJOmOxj zMa3d^SARd;d}tKoCT-)(YF5T({`kuXF|S~Aey-7al%KHn;mglriY{Y^Zc8WjsMfG6 zjr4y6Iw!E12YZ%(<+ZUFKl;cgG&Zr{M)xl{5YJ+L>48GFrg~UAPu5=6{4eag7LVgp8ccr2G(XHQ7eSm zV$%^}Byyv~xKLW`T*kehm>TUZz0d0oSV8WX$Frv-IQC-G>aFt@Hnf@)oTu)JZEY!r z@g1kY>k>AvF68WD$r1Y&3XeUoBCG3ZF}>8dYn9~lHL1VY)o6izu@-MEYFY>%BBjS& zFY`Zi5+uNHKYBnO1pe5WL6!1*{tP&>Av;u0Ai`;^x*IMy2VqmC?+dgPnecYT(+jE! z$MO8E7hg~Cgke!PspQVsp2q3^1KlrI$Z$JqF==b!NUW>KrKf+F74Ntg)x|PMiGSgi z9mA|*uoUt{mFcK+_}RmXkdi)HoWf*2efmlQ7WP=DFZw+vp4TUudz=3RZnDuL@NFdp z<8zBRa`QZo=P6A&egBUUZ%QEXCFIM%elW1G^*Hk3EWc7@OKVTz85vVV%{tlG$}ku6 zLDxmxX}|4}2N4Ti^AK1kkMb~!p^mA=3}O8B`}ysLx9s@MCBf;TD@B;IJlUAnhA2+6 z`Ha0Gfdj9_dpEI*Wms)G$6=zC1a8)7AzHV`g?pIt;9h6$VePcr72Rj0ahY0E!PhT% zaIN09m)Bw{v9#_1RXKS%eD$jK&(ett_-U0!LFULBZ0K6Evh+(uoIC&aSzbFqyimVX z{~JR+rkpg~D*PGYO~175rW-|Y*N^S?wY`m)xA6R9N5)HdO{|cj=(ZT{5arUQ^Zg0t zW}VT-m#2d3-z9VT6(x!Db|hp}W)|Yc zSNN#g0blK^Udv4C#geN8Lan|S;I*SwBd?j1adV3o85eE)u$q3e0$TU0xB`!$rs6YI zygh~Gk@?sl_FL!PUa7h%?hD@1S@r68c1RKTe#i*ML7m&{9Abe_PzH7%C)LL7Py2nk zck2U2qS5rlp7J_wI`M|{kBu(=K>wGN+3F{ZmMe>%kJcKGN|tn@4AaLgc=NRB7bmf6 z=Q(Jmb!~BdQKPD^X+vDoOIlI5Z5C7at9o5JW{-c0(PTIlXN~ zN6*IaR8juM+8)qvXAk?~U)9+RotGW(@u}=RcCH<4f5PBIU+OKq@UrXk!U<=*QAW|c zvVI@aFpRxdP#lblt^@`j=X1yF=Dek#`w$~(=aq9N4#SO!olkyV^}uPpG>N#L6XHL_ zJ;G0Z3&%TU`ex_$z3@Eu;^j%FV|Y-$3PbhrD7-DP2Uk(^#a#&|8}vg+@w=CdO=+xR z@Kf^PA|qq|_<5G<`IpBjar&n9i=BVs@Zh)K6=)v^;*B(6YN~X!IDy?W-ydHR@k874 z7OjQbIDyMqCXP?^_=CI@ZV%Q}e1@sUX!2vFRfMYxMh-#kCXH3$u&jb zme{uH<;YX`pEdV{L&HqGF2JDW<6;!v@+odLcjpW)*PED17@v)^&dQ{Yo5tea94}h0 zd}GJAUw-w2{aoDQ`cMTeLjqn_a3`+&F$ZpY>Qyw=qdRzL>fqn?)5-Y8o7KCcD%`jL zN5?Dex+1)jZAIH$J`GP-eLm1F#*3@_j}^TOEyew(3O8gDLO*1v53^?p&*W!W`xJ$*%CL zkDnl%g72E;hwqK}zC6KpHeP|(H|L-BC%J?RWu1Q`|LZZnsM6)3kXwl#pL%*$=e`QQ zuDW{NK&S=ZA-d}_Em(~wJuKR`D7lQYe{=qNvHTehLw4kyD}55?OW9k|jj z%_FzYCfq|%%f~j@5Pw$Q&Q5pg6)w7e%O-jD2~OD6J9b6E821-2Cy1fz#(8Z2E8LQN zijUQ9{YPtThTjixu}qWi#ZL`bm>%>#$1lHW`Tk|h68C3(@6lcO2ER>K7UEy{65qXY zbD-wybzC|mN`gbSA8*Uv_&jaVj_aP()qHW*2B+~9O`~QS#IH2v%~~>a;v^@Zac!5` z;rnvJ=}%aPalUfz<=b1Y@S`%@`z*eWIQbBj@w?FxJdu!K#FVrLcMvpL`*GI=rx7|? z_^NpnAJav4PTH?=3%$peu1MX)k6z{8H0k+>3oriJjq87llm7lFSy$wVi)`8`7LtC( zhy5Sd5hlFD)4VR`_0M?YeO+p#+aVLU&(m(6H!4H;0@EXD5+^^L^FWrCUvwIed>mlY zx;l)9Z|z64z%AUTvog3oeHI_=(3QcSzsIKrM}q{sgYj)FtT|Wk3tk{h9rpL}2fVU3 zj(X}J{wFd=*>%?EapnM=e`{t8_sSPycRL<|2hZrFpC12;)BmTE5{!@IH59bG70OX~ z=-w#t?ad|p&Dlh8&bUKYjIIlz-|t|K~Pln>Axg=8sarZoY(Pby<@bgr{;0ulz`Y)jx?N~`@)q4p-r4& zoVbPI=K{_b^VM*OIRih((jhna<|ocd!L=;CyNK8Ik5MISWa0jqmgO69zj1gsc-SHM z4KKX6Lmluk2k&0mWm3JljZa6b`}E|j;^EQSZ+Zgq@w@44)1+a$_@0XRlUnID+>yCe z#UrW^AMXuftNm{uw;%l@RFtrR+hT*t=vxV{b)idsL+AjPIi0(ycIGF3?}QtZ-={MC z)P?wz_oGMnS$Cz6h3|gh0^#}VjK9k9%{__>24;jP*f7#Czj6z|tthH*t8pLqrWM$Y z`bLC~|JZHttJ%S&=xPNQ7arix-e7*cDDE9Wy8l2*Ur^9as|)v;ZWg)M=2*zSAJ9ze-9~=5~RQFn*S27 zxOZ33nwJqZ%QYYTT_8s?Ay4%eecEx47i7iEb0^UU5|?EUQ7UxOjO3<;Mkj74eN9*A zzf&mXaf2t@BN}w<9WGy|*M;A$zY<>8aT=)wc z0b&)Fe5mzda#SGSIn+6?;dG_^Gp?g(-+yS#k6QQl&M-%Fpxev%c)sZ*K7P1+d@w)& zv6`-d&H*QS7%t5A#ApUzPZP-P>l8xMj5qNMncQe2qp_SrXAW0h?g?`#6hS!>Bm*Ot zc~JO$%M6@t9+&nF%_dwIMZV2u+QsMikc6Vz3fJkcxVa@p{nktz)%AtZld$n4$19(O zPN^;7=llg+!{;PWaI@D)9TGs#+2DL<%{LsKek@B6ErnjY3oAOO2%(0StVfbs-|?U| z;~WJw8N_~gY_u>+1YNUh!hQ*^<26>|}El9wV$Pn$Bt@L139pq;xOj1Xu ziiVw&5PX}h&s|kO^0Lq{PF7OxU4;` z#A8}S5oFC+re%QY7dQKeeKirhnx^_=LXXJbDGRXF8=_Uyz2R$-+9)+WIOfmZ3AB&1 z^BZ2gibNUG`2w1CkZ?!&%2pdA%3s+TYHu?}1?^{tuafDZM7!Y7gIXqZ+G8eR#oZJY z1iVs38v3Y;@6g%0ml-KgsknH^nWIMnvb~yj4bW|;&&Qi-Sx|N|g#sD3C0Y`zR-L6c zLazqQeFqcRkn)X!>`*}~G#&8L^46WJh0q_IKgDcuyq#?6t)J3$VyuiU6@ z^Zlg{7Pe?gQr)t1%o0iNgg5?v#e=M(mLgNz?NC?pQ|ZBNDcL8N=+_5c5Knzp#`B6R-!0nrWX6g?v6}atl3gr?a}%f zFX5g!F+_DNzJ`(31AWWCBg?<%fNCVIn39qZr^&diQCUhUYwgb>3@`wI!cfHV{k|+;$UJ1$i9PboQ2}HO5GAe$S_d&v< z*3Og!Dxj`J^NLLFIqdla*9G-1-Wop1q_pfpwQ|(*ZYTl$n3AkFAMM-s?5khFZ^g9AsB^^FgxM_siek}-xo< zGe+}d)Ktvsi71J){;{`9G_u0|$a^+TP!b=v)#_>z+N8IPwCIdMB7>tf7rvUIa#I)C z|0+_@+-IWMzXEaSaQK`9(SZdTIb#3ySw0P2uq(>DBsih&xra9fKd$wrry7GPh%xt8UG!#94DIJ;KmGo9h zwMT@m5`WN04thzZVY{K7ftZex2ws6BiWi#U{K%1qrYOl%AKPalLmTh>#uLuShlbIJ zC^8?-a6Ma{OUXiQk4A<_kGY~biCBUB?K>!UIy6DPHyfp!^pt%(<&M0z?{4_rE<~7P z*Uu!nT=d(R_75MQ2U=BDnhRqtM*Z?`=s1ge-}-86$k$=C`KH*C1>J(`5|AL)Je9d-vIFDJng&Klv~(0M!qd zVMYD-kOO^ucMeM#!hbR4jF<+ZcNV!zGiw#-Pb}kpG2>m-#(YDjP9_+&V%wY%XDbn_ z?@!&dGv%oM2FDXB=@66_rz*!|^#IM9zGXTlbPvtC$ewV%7=~nKxhs6?s!-&IiLQK| z3ba9#8BC-dj>P^r$+6PZAb(Ob$(_*q=-z>UsbzQsvLAD_9tnAfRK&5$ce8nRXPLh=zo`i^R>E&@_>;rNq^GM6%>drfr>!60&?|ubMtaR!{#f zypw7`Hg_Us>CI9RN0(rasqGWA;+}TNn6nY(4+uF1KTSi`=jx501vMkuci#T$r<+iP z<^Nn!kCFGfi^5&POmxS<=uzJAQ#27&GiQDA37Y=aLn94Y zNPFi6)iHtB zDBS*%*sf?d`uR$^x#mYDn*L;}g6qCPyU)7LpHt{TZeB?=%{EmieEisN*?em^bwy$r(dekHp!|U>{H|P4b5o zwN^y&RH5>C$tRS1+%VWh>LW@!-QqLxJPC7QI68VipnipU_U%fD@Q-X&!XzU3w5VNrjevnez*6BZWMlfgW{vY9CGnz=533dL3c!d zWn_``q6JdtTk+d-=$6{U9zCyFJvndC z-6&0irI2~_hKVy}YiSMNgIs_^rKx zz8qg-kB(bI@9?eTakS%zB0BaK@#Y5F`O8xQ-RmgqSR!G^wF#7yEY$ZPe-pK<6}%~p z*+3$^FDzMJPND%yhr865eXA-$MjSLx{U>1LWq9@&SMdGuY$h!Q! zL{ZfY8g)t-ZfpL9n5|;Ip3V7*c%Dw~MikDX*PkDR+xq`T>!z9=*Y$rP+f);0a<@5D ztsLp@k8Pm?SH_U@Cx0VbJ<-*ts$b9=b5Roy)i!z^PH@ZfaZ3)%gh)LBA+()Ol94LKc_Ym_C z$FXwBWfa7!^t0Og4-)yQ8elHGkCKO7d6!MUq59JY%!Iyw(aD|FyN-f?&@&NqIUr#L z`O*ILj}JLO{b~d&qvC(jVb_Dr?%q`d?n*^LX@|(l@Ni;W=>U~;6rEN(@dMp(7nPu^ zIYO!IS#r&0hv>2VoVc&^8k%Snh(Zno;AePnP#$!I?wOucc{j9X(fq7DeaeEEXF zral4K+3oBv-Q7S8WUtP9Sa||&UtE%TkeIkM%-KSR)CVnCE z4NWa6$zyPhtg7(Czy2_VZtnm7`;9m#oGxpAI|c=@MXfn<$3RwMlXSCs3r$ZfWWK$7 z9EA5R<5qtigLX)K=2^Op2#?XYe!M{fnVoU6yN$tz-v9?3LA~^N@nu-R#G6 z{v=S6s5e~lU>Eh4%YQt(NeZR^OyfHxQuteyXz=^R9x@$hJzmj427R!$JI6o|w&l#pZ1l2uSk0pSakoJ)R(D7M{|SY3$<(jDI9Oy^O;`=IJ46P-thKQ!jm2Qg}R zLI0|oA(09sQf+C=k_q6Y?w8u=^Ga?A;o{zs*L6vjlptHsYXmySO!=CQ?-V=HlXQE;I;6x99E-@^8=ROYKLdtI+ zs+|B#s!zA#fF6t|!XAVXk$_0pzBTFe3E;>iv$>mn0&e%TQIIZ^z>MUq^~C@Nu<)?= zS`cLb9~UZ{mmf$0weuE)GciJ-U;E$45eC4&e(?X&M+P5WyGC(5V+387oSA%oMo{or zifDdI4mMn=KRK_Tgvz{|COF4QcvX=k6*NKt)8_mqiz%7l_Mpu1l);lg@$pfo_dF$7 zpCu2ddBFr$+l31y5lpc5Zu9fP0ToQTXA#?aoC3xZCRbDCPJy}6^JYd)8YoP_4rVVf zL+#Tj7XO`7K;MYYY3kB~b<8(Xg=J2Y$x}S78ChZFDgoh@Z5D_sjLa-m?9r9?tYCLC*jMu&1NbR; zJg|Dp28yK(1~koVP^h+2RBO!$#B3Go&uZCWalB@vs+b)<{alsS;XMie{rLH?Kj$o@ z$%Z!J$!8(2t}Vs<=Sc_&zUMZsbPj6k3@`7Xb1>9Q)qnVk3DjIjR@Gjf1A~?yP0s(F z1Nrw{`*ww=;M>UkhcD$gAR&O@9XZYcl}FMfZ9&YSDX`bB@qz=sLaFaUCkL1Wr*z6# zpN4~oP5VQP6T(%3Iw&+afn=>>nMn5x^cFMGEO&E)+3e_6(kLe^NgA?TQDK3cs`RU; zbh#k+LvVbMH5YuX9$(f}VukoC^7MsMT<|#f75CCFE+D_Q(J!RN25dU_@fde*IQrsT zHXFwc7f*JI@f)(kSRyauBGGw}p39=R#c>`ME0_GdoX!Gw$i!yfo%8U5r>?N$<$0ix zs`M$obPggps$4Hi@c^H~)$Tq^9ylSAv^4Ya9GGe-sk?sQfre0Jt|~%a2%VD_%XjC1 z%T64(F9+}f8#{lDOBF9%vg5tENWlqI&!dFa*!iG>GJTEPhz~Nd*`yX~IKg+y{;~TT zKFB6`n!I%E0`O<-7e3PC0+{y6nee{=Uv9f}cR#%VEAvS*_J6s6wu~@Ym7gCXrzdzv zJorKVv;`Z_LvGk;+YalP<%j2i9{y(qFM`NZ*IDWt=OJaJGx1u%MIfBdz2E-hBBWR! zJi_v-`b2N_b%} zGM)VSf*^RFBWyO;5Q3|I8hiE5d_X7REtyv@1WmMp+0<;pFn)9~LqYxm(EU8>+X@o~ zpFEMzKEH)Q_mEAtk)9u}Uz;WV#V-QIZq%>!9Yx^l^TyK>Yy7Z$6nf#mToG{95&QG` zlL(OIoc_8nbP*Q1hq)F>M8WX9k^H!(C={5MjVwJA03Gs(F`ocYXf0_Nwd@du`6jW_ z>s5l#&NPR7zKcSe++Ps|fkku6mCq8;QQBen8YJP| zP5|4rATih=Qa#RoT@sk~Xa3IjOG4_QsbDs(IH+76NILmV5`Gm+-kCR(0_Gs{-KI)$ zaCQi4pYV|a%jdR6Bdby{UByQ_cu4~6TCF03h^3+T0jukaP-&nemc8(2UIM-y+Zi~@ zl7^w|C6Q_-8Hf+A^wf@%1PNl|ChrR}uytyOqPkuN+Wp_lOo&JU1HFUMmv$MrKD_SR zqb&>N#`LVur=*~ylIZ7Eb6Fr@drU~VCJR$99vmE{NCVSs!uaH_?7wvi>PM!_LAuqz zu~*76@bIU#UueD@q@6f*R$fLPuHVToI(;Mq1W&z{G!*55#;#Kw`yvm>Q4R zA%R)GwG#qabz2A3F$CowX-Hxv5HMEEGsnC|K)LzxBEJGaf5ZejQ-gq@K;hagL4YN- z+bJ7&0ebreor1Ulj0%JElrjMRBjEg8La7963^_V2r2v`sjs(lQ0E^S>mYAOskkz6u z<7buNG_lF2n>Uo8D^kyY5@;9u?n^zo1X7&*=d`Zp$gdm#iqdMJ8o9 z8c83Ot5OEzhKW}{UMs^J!?H`V@ygKTvLf^6tO~qrYuC#$Qvn|$cA4x0W!MVw*dEJM zf!Bf`TcXxgfaj^Z2Z4{uKYTnb8e>oeir|6Oo>)~V-RkfU|E>bf6lpyK>8j9Gg%P`O zsez=RWQXQmZX`1Kf)qsn6O5;tzWpJg8y!8FCIt+<@J}v2R8T@LA zM$dDrLl@myOFe51XmjdCUEeQ*8^xQ0_jT%!oG~D;fKHMRkwa!P}ybhU2D?^@i=47TX%F}*HNI}kxLKqzP~x!ENK8<*UqWbt?NR2 z_`pbekRCAb9}^s^H30n;8eTb7JJ~X5TXVT#TKD0T88k5NhyWg~7AC^6``gz$C&(E?Hv$C%HN}h#Ia!`rC(h zotzEev)X`0AhjV3^qe%JZZw&xlJGCI=X$VLoT$QWG7>F9UABPASLJW<- zO`ng3uyTDiyU5T4nvLGJK2J3S8Rc)WB8o-YIS6w^MRX7u zgO=z{&U&{oyqXkY5|6QfiY$TTY7b-R*^Oh#mNtP8!3)=B`7Pn|E$oPI))*vuf5WG0 z6Ht7bs~tRI33H!fijbZOw1t*6ww^VGHt~x)w70K;@TZSLTdzzY@#*PPsyU{>euFOX z9={c&UEdDUkTC_b=%dpCG-e?G{q&*IgcY1z%+XPIGtR*e|{6GI#enREk^wCsASs$8R5G@PwN~&k=>;#L;!&?4#HfJ7W&M z#pegk5?H{eC4%R!*KWYXeWtC29COg@tfm|dvw#mzh1Kr-cLT)ONaC-Yw16;|*MC%r zE#Y@-eU;-mYj{Y%6X%_60X7v8I>E7)&=g#mIT>jUf9CS~j89nte^bFzf9h+%&B^~_ zf7u$&grB=`bjK3ZkK(^6WM6~IwF~3DIyO+yp%CPD?iy%nx@QlvTS0rbdEMowHb6E+ zQFOoZ8Wf+S@DQ)Ig1*6Xd*Ph6U^7kOQYT^sJ}|-5%yJz>q;G}>q}W0U;mKEVFRVb! ztF{BPxejlLUXY9awFM#$-cDBT>!2qjmieUSIuO#$*wnk&LCt!KnzhSyU{Akh&$@LT zw1r75|4iD!8{70RIZfALjLVwI7rz01c3Hhly7rJ#nHAOg_d2Xi-3{qTxB=3Jbhk*m z?IA8K?d%8D8z7_hERcBY1{C#EX^hA_fRY3G{V(Y^K(N(yQGweULX%H01wVFx<{C=M zrKua>Jo>2Urk6EDh$Orq7j}fz7d%P}{MK-DwJ^A-!x~1_?z4o|I08wGBR&~m4e$HB ztX|XGz)qd^nS~2ZAh~0*Gt_Ghv%%?{Kdo(`la8cVuG$I8bGNQ_ow9*;G5@=)Pi)}T zz3J3A0cZGlsn+GGqYVh%c=zItF0R>s><_zx3L4=zH*5j zI(x`5EGRYEa)WxZ#|KsGYoO=;C5#8j`>M@5U6svv)*G5j|aLr zf3e<#3=PwX+MD*kWm~!$%jEzv$3`6gG~I;NpewG7efE%cW8Gge+yS1gR(#6i@qjD! zWNPWW4!|TE%DO+}08aF8L$rN7U|x*djwjjyHo0y$kKvB+i&Tq`{gnrt`V!v$Vcy{% z&z|BMDs=>pAuD5Q22aqH&q zd4cqaKy&&5Cup?f%D7&W{V~5W0==>S!W(sde&R*o^ zx#JAH*R)%B99)1vfh(fJ+8Yd?A5^aqx7M4R{yW!~b&RYrNxrb2 zeZM9u-VIJ{#^ZtZ?qEO9X>j3(FZ91|G^^Qj1Marh4;MbUgEE8lZ#o@6pwfLoMCIZR zuel=ri`BXbw{|P+1#` z2b7zAe-;$%4`bh~nu~gF0;NpWQx{_oh%y*-N}2ZuIYS%9+w2~2&gl6>OQr`f*zHOG z)C>UIJ13R3%{;(>a(M6FkO%zA^nO(JBmms4b)DLCJYX}*nNWhp6Baganzx+01+PLy zBgfwV(-$4KIXf*+2su#hWlO#V^{lIDne?9j;!M%8@`;{ME#?$5vwsU_mtLiB>3D*l zpz{2ax1K=Q?@o2(6bNs0PemW5ctWiP@%#r0FW7y2^fzQY5Ee#^oF5H%f_RU&Z=t3a z{C@a(;)YfbROlQBa#}ASrqzw7NcIASg0~XYZ9#C{zk`fg*9$h5EvPvLy`Yt$Hq%TZ z7&0-1ju+`(Ae`sCF-!kXZ>a2|e5-?D)*f?a8ufzc(^ok(jl4m;Yck`2<2DG=J|Jd4 zcWs)d2p z%b|qyBOmZPv0Hf*=L;fb8DDMU!(e{+QyuqZUm!3tEJ_>kg}!iQtiLY|JTJBNxTX8T z8HO{mA6Wh1pyk-y4ox_min{gQZNe8e&d|9GIQT()dI`g_emI;~VgGoZ#}B3iSh;RD z`N1ox&n(5c;b2&^IGE$<2L`>$dS8zGph4daXPF8Ijx^!Uj+cHg`?vAyq1HdW9H!~I z&J_Xwk{orF6#l?u_;~339e-$j>ac#(IRaV=J#y4d{6VmxmB{*=KNQoPCe>(+0D;Zf zpa=K;K}_AzwOS?s8t0ztZtO<@=kKoT-CO>UUVQtoHZA}ZttRPhl_SA27^#)11i-tP z^(wc|0Z>YGEbV?uBoyVAiz?;>fVAw4)e8SDI5qx=?%C%^aA$Yryt5Plx9yAR&)>cU zr{|157oLv-UIR7OXR^0o+@ztNedrc&dDsdYctycsA_HeZ(k*z&IxO^wD-b$6luV}D zqF{SzP6}pkK}W-n+xh{4P?V&XuRs$G$uzr`4UP;3mit+NyIKXi`S#jMs02!(iyQQLQ-8h`Dm<#KYTg zkNAMVR4W9ig0mUw@5RBc*iCI2@(>uikn|5TMMypV&zj z3eP&lpD-H7!|{73wD-D0fN-&pUcfOFLcB6tN*m%~x9Ngz&iPQFb4Yz+`#Kap7Jq$RMbj=l>A~447=64+;n6 zp$EG3l}XUj&UVf8TsWxdZT+404hJ{)*u>xD$uRC$EU;}54mKVCoq19j4w$~PO`lsb zL?(re62^tYis+lNor!RG?Oo5FJCY1lt{M};t>Mr#6-8S~7XfFVO%R*Pq(H&QvmS@d zaQIG8I`>K?0tQ_(O;41hfCOtd>nolJxTr#wb0sJOu65qG?;=Ts^Cg8gsTL9c?AiBf zwT1}Dcck_5xS0wzgbz91WkkS_)YjrhOA+9F;dke$[*RWKF&CISe43vR@)M?!F) z(a;n}8Yn#0V=p0y1mDwg>;p!T@YUL?s?InKjCp@MOUOjRm<`8+=;TQFNbW~>l$`c2 z9{5c|(IXO!3|y(}x+5X+AggMyD-Be0##Np_hy?M2@9)$O|KW3Ql$qmi8uY(uny;CQ zgo)^OCj;py`2Nt1=(l+KKfn0?!xJn~U|FbpW7s1K_@!GY)Lqlz&IA-U&a-H&R6pVWF`5&H*hEfL+ zb#vwnSSO_u)H)swEv5TQ8s^ae4z6--`Waw$1{*mB(U2yfAy1VX4R6NucH}fnSw7-*QuG7W<+kr_n_3=!QsJbmfyb}#H zITwmfXvaXyVC#nj`1fDasG?LN1`c!hvL7bIz~x6pKT<<7p-@JnL>MlU7%a1;NJq3Lj*$ z;Pu{$?pNVh;CnDj+Z7!PZO&Vf1OZv_@z&Yx?>GPP!Psrrr``W>UKzd9+>!tlg6@7o9uNgR|(U&Z!zvcMuc^P$gXEU2z>y`fc$gPS9zN4+xH@Fu?Z)0S8qfLi7U z@3=VlJ3sM-D=-`AePnmHeBwZqha}YCZ5-5GYfYtY%?AD?!+hVCI8e`ito4X09*$C_ z>!J^{VOo7j!*@RpUUWfEtWG@8??o*pK@PO;aTjhW#>1t&m)=tu@et2{ZjLKD2V4v) zS$D$YK}F{NRQ<pokFFTY(6X9ZR%%dWnJa~}wv*_tRUil&GImZx?2rHSDp_1-- z|DKu9fxc@Zm^nGg50YPQ73lR%bZrls~`GGsI#i8Ars0hh*?RlD5(c<7-%=}EujfAN8T z&y&6H0RJjcwyH-ma0hfhC+tWDA^decXZsxx&its&-kJ;*`^z_cs8ZnP<_k4C>H=7( zXLJc8O938iJIT!?1tNU72=mPfAj*Oz%FQSRxcnaV|G1X|87(P^FCP>@f@8nP;@uRu znR}ygYcmC`F79N9>=eK!id6E|jTE5R_&OMkr$X&O$Ka%DA%w-~RR;f)dD_Zbijz|z z3h!9{lv4=U$At2NlvL1vm;M(`rb5;`>2&Xrz@txg6DuiY(f#3KK$xjiKoH-H>24`yEHH`e>_AuRRn>mE!BpO zX>jjubJ^*XH0Y2nyUQk04E<(FY#QllVEx50x~Dk}I+@Q0Pec~OGD|;w-19ULa;g0I zY$gr#bhPLcJ{7~~-@G4b7SiCf+5YS=vUH%2D?5P%OCbDu;j?M#bQt|5z2PI44g<@= z(#sJgAl2bXxhtIxTAtCmiWcc`LO5DX@^cB~$|$m!Tc?9czws6Am~{BI^m&qjxfBXQ z=VZu|(}6h2N2;Pe9T=1NVfIof4E~q;_^2fv2r_;esZ69paHR?dfk!F4vVT1d^XX8w zD7P<2k^z0I;h($NRKuL9Qr~I2zK*FvOm2w&Xc+Z_|+#&;- zm(IwWY?ng1v}E z4X3gV;&?(7DpW*;q>^~EMkW1AvSg2>BqBvByOgp-L_&7i_ua7``#Cs{ec$&eg-DA? z)O$WX*EMs`Joks^n)~<6Ju|9bX5!RC-$H3x+ROEHA{v~JaS12C#a!sqoBQhCVt>{_ z4>hYqoN(0}FYczr-3D8hli%K=rjXrF#69o=JKad zNPd~1zq~07FW)T<3DhS;A(8_~@Ghtw9Gb}x~9tP&eWB0P!QeYC$oU@B390MB>d=dYJLlE&-C+Sm>2}8~6 z8sWI8BxZ%%;qVQ!|7&+D6;U%AxxYNZF-wT}ekME|%Hr*HcWhE2Chd6iQ)xKnHOSlj zTEe08<0V^JdMdKK7r8_i!XdZ5&V77S1XliifYoFwtP3`V8V^R`<4HCve(4C9_XJU0 zxzo@&m)cTp5P{$Ye&f5A5fC>_Vf?0*hKo^#s`PInKzFfElOa9=Q_P9CDZXiNJGWAp zNQpq*BbkVweG%~Z;i&SbHVqw`&37W#Be3z5^j3~F64D>t_};FkA!u6=`IBTMj1Jt2 z?KvNbb|*iL*D~qoy7gEvz$y|)-0gPNJ&uI6NVme=qja#hAN=B<6bVc279D0nBzT_D z8D(arL-=IelMf>_Kh~KgNuMKOK`*ST{T&L9%D__a)*ZF9b0APB zWshTX6tMMWK&In^Nq$*{)3Zt$#7ol#a%Km(SH&%SXfT z!KTfhv@oDkQoSEIMPqJZcTmThOjyNK`l)6`^F`N_FVV0P5ohRU%|h>jcT&m0clgL+%O;OzO#Xdtc6`-3Q#Jr0K>W*+S#){wpzv9Z^;- ze4UMw63dIEpct?&<%(#~o+na`PpNifV>wtrTDCI={NZmfUTW%e9I?K-uw zV@E8WsIL5^u*Ra_>hi^z&KnN_S0B@8HaVvg5(qPxzIL$ z@F}D#4s1P-kGM9*;YFUt;iwaNaDDSps%mRISk8M2<}t+MD~Elatz#Z44)RmpAB#u- zn;8}@p?FN+ohdD^%)@)`+t+0;$0N7?Ur#(PjAsgJx5Q&YA@-|LP(H*Yrq4_)$D>64S^)P#JSu4D(^s62L;m)@ceoRf zp=j04z9Rt=B(k>ma6U!^3e5ygC*UIAZC-it1bpt*dH?8NK6DfZn0DSwfcjmQoG9G{ z=vF!ygdQS*vC}kL!ku<5g*jCmb4)-<%IfkD6#_=%1^j)J6EI@I@O&X60mma(R{b9m zV6f%b`(2F*n9?ygFPKL4gTcHHZpB($xJD&jSrQ6J`6A4&XF%^ASL4beLr{1O= ziFjgVdhYJ#L?o&PjuYP#@SOkr><`IA%r6eHuL>oi?pMuEdFBG>em-$i*B}w&obO61 zR1&dyI4b|ISONHbE8I-zH6+5WOKC!(BoUA4-fCwg7eI|KCTV^?5zpnW?(G{+q}}uG6~EqEfWwoF z&#$p3p=Dh0a@@Z}2!FHqtMan|Hv-oV{gt7~S@H^F4<%un)t+^SeT7&|o$4OCk%Z4L zJpFo9k}%@NVLhQ(h>#$MhA$pTSX=$OBMA6 z(R}_JrIPlwXIvL_0-aUWYl#H{EA9voLDS~XSml%m91l+K;Vz zr9f-wq{GpBDToW)Zetx%jOU*wlSJ}TkT>Am86A{@jm96dfuv%DyV=@D4yHg#oYz~q zJcV|irAhw!QjCLZe8J~8r(%A=dq8YD1=6p!hR5tJ!HILS=Z+jt#n)7qEt;&Us3Bc7 zdZ|5ONMQ{i$dKu$2B1e;1# zbmg;B;XTiJoy|8DuMTug_qCVc@doESM}I0l?6HpEDon*he#zd}^%Be;WV=#^Y9jtQn);g<#0Hj28lB~Lt{*7SdCTFdvUE4Sv~%un+?(+=3`+NC6R{D zaVFoDd`hwRdvjzk3R(ek%6O-G+=to_=_bO6FJ(YDdqwyTO+PcWW*Pgg48jlp9;kB2 zz=wYUW(-OhP!wl*pVMB3EwT9zJd!g&JZ~}k-Yx^=-+$^D>B?c|zW<1C4^6&Aq5GGR zfs^;In5RmXL*-Y(+pf)-D2viwx!#!pf;bab;GJ@q-uwISxknRD#={Wk+7Mmv$Z z$Z}}jxJo=~nu&a0+r)m+Oq|fVy{yniTR$(oRW>*ikv(S@e%{E$ftsC-Tj(m#t?Q%k zsx}h}htd?T_-DfPy}$~WbOp>^64Uv9WJ2@%*^fUfGO?*PtKaTk1zy{p4PF(@LPN`j z@W6Z~EIlMDcw#HCb&IR&ZoMq94k{S_;mv~jO2%*2o(e3Tnc;Kv&Vs*TL#@f>ESyqw z{G`fUNh_$iXTnR&LX$Kxf17I-t~lQx?v<~^BLA(MjB{C7dyr`JD<=!LM0$>9*;T^q zqwgK5eY9aZ1nPDH}cy= zXf7pxOFHS(y_P0nX0N#;>az$O*023B=O)bmK^X68~oY$ zSp^e?xa&fmIVk5jlwj7%`a~_9@MQ=bR@e?JgPl0au*{&8zX?%gm_6 zmmEC5khZHOB?o<`k&0EZ)i`8wv`0uF7c;_@2AMSe_?D}B?)yMBxGMdOO!RYM#C2Yx zlPwn)en>SqvDbh=o$37OmkS$jp`T~Y=VIbUw1kFQ4SaN-dlZs$k>g7Gw&a)#mwTtD z9K33vNMdPhUde?l@5r++S-EKUQvv~YbTXY^67M(lZT+YKMOI_t;SrUeB`B-f4$%BIF z&e?2_JXqICGJETjz#z6wSDi5*&x^xObC%^njr8wDkUa^vOnJ}C3FO1-sQv!BuXz|! zX<2LtCZU+koHpPmA1_?Uwl_KRF~vCNaF{>>k@*cgujRx4>r01C|K(#V;j3zX2MPH* zo@Tvw$cMM-LM>G_AD*`Pv1W@TjIyal+zHM{a|FxFqC5H6C62@kY-H@erz6dgmyfTG zGK&o!`MC9Y#asP28GqD|<$q|&ho~woO%zM>_t)x~hZ-4+e-dTJKht=MaZM1o70!&e6;*70p1VMQgN&n50^KYF%k)o&K^2_#G8ON zb_N?QwORynKJJt6C!j1#wY(yUfP$j|KiKcrVpc!*ckE9B?DiIPJCO-6|HAch)9YH~ zSd45E<}85hHmRcP69n+r9=NqHw-y;jF1(}SG=A}t2(Z{7AT`vy(7dA-ZaXi!5VZ>6 z)1w|3&07E-$>rkFrCKQ7mDZ?!SOBMHv$++y0thDTi*e_ugITUu<3d0IBKL>3J~b$S z>l2&UD^hh>8qepslv{x8tkV|9oC?sV-uYEezYb^DY#p213c%Fwog20QDRBA!3ttxN>xhxSM(*-m3V{u5c8hDeHF{4}Cqxz7z8XX?(QC>pvPJ zQHZ^IMptAG*Mpn2+u}w*Avg~Pzqi(*@snU-m%BzXC`4k~!+HPULa-ar^JfIsqfkL};mcAX#_x9OjTRK*Y}o-*ucCU4Zx;yS zqj8gk_oW#9oc_TmZW<@ zt-FeFR!A5M0u;=B5vU_O7C~gXIHgg#2;0lHXG~nA;F_B2yVQsxOzpXTXGy;ZW4ZiU z1WO8pNzFdQ@*<4S^)M+q7lDp()V=r(1zvmTt5-&gu>Qymx_tcM%Uns_kS=Pc= zn+9;Yt`#m*3HcU7?yt_`x9SGCj0$Mq zTP}uy|)3ooh~DNY=HHTD`XWKe@PByJ`?OH#slpa1$C@cs6P)*xh7u%Iu7et zi|@tIRAlG=BT0pC!&jL|lM?WWZaX%-y9E8bTt_XgQxP83n#=560uITy*Ph6gU_E}H z0JA3*?~8hr(lbjiu3&1eW>kXonKg`OQX#VIkF0Y?32Mf8yj9&xz)Bcvn(d~-!0<(l z^Uo3(-pc&xoKgapu?wrm{!o#5Y-3mE-cnpWr##o+RD!JLf2@gujdypNuVcyjdyQM(MNoYkci#ZgO<5%rLc8sVjNd&!qMrKJW*95 zlvft1fj>&|=DfoLwuep73UX@Xw;_VbU9B>CFAFCD#nkGMdt@PjC zX2>U-N_^d3h67q3#AXX=da1V*)r-yWnQr#}AX5g$t#Qe#!$g!g}yC2Soc{euW_j#|je|cr_yzIWu zP`eBr-><9?e>Nk(BmYR_U>PR=%ZN6QD&$}uIAusQNVIc(kwNY!p`#X?rJ#D}hO_;hv? z3HQsv8Y99TdZHDcY0IHU*2{5ZFdKKUo#Y=zN{(j$q7E1;FS{L!R5E)}h`_^UrH|Fi;h3PK0OI4h7o@`u(W(F%=3 zw-DpB3iR-qPhXLz=}S1PFYax_z2z*AvhE6)aB$fcSX98{F_W&3avL)048N=Wt$^~8 zqimOgD^L^=CArnP4L=twj&dEV#JArgD;!l7;9K~X*_hUbF?S|QUj0f8^0qnb`CNfV zAt4OZfi@@|rLG%zR6^5jp;?==5{g^y`!g`N!#^b8aA$5MT0VQ9rkt(BLsET!yKFno z{frZJAE^XU`!uVKWhG`xH4M()YlpFMgu#90Ds&n%dgz5!;$-iBnbT42;8%IH7$jYV z&F!*NwzZWAwaRpGZfS?gw=yq!* zxC879bMtjP)o`k`;xJFF!uqi+{wI_UC|(GiRadEoylh6%SVt911FCqj+JWXAHuYW3 z)sWrqdLn3}3U_?}b|&xcM0fbfTU=??u=3|@e=brD#aI>()&Dw?IiP2KwZ9r>1WsY+a}_p`;puDW_@} zT{~(2H)%I~%xkcueOl?w$7*=pajlUE@5E(aB5kHx4T2iidop&_V1&nM&7IJR)gy(W zQVliOa{crfIz<}SvQ$O>Xz4`8Du34JUp2Hk9%56dT@B3Lc&IODI>9pQIJ3Y>Lf+;O zszp)__Rs}opQYPdgzo*pI9>aQQaGY%zT>ez+C`H?PsH3+B~SSI1w(ckU2b4j%L zYZyMY(uI&q73(28$jIs+Jncgzp@hwKdV{l@_CI;V-cFc|qQK^DZIdMIS!S=7KhX_4 z%PbAM3uIJ$E4VGOMnc`ka`VK?-7vYYR~%wQhJ8=Q10p9G(K~dPMecQDtlpw)z=;gf z>or|tQ8M1;evKRU>W1$o!`3swWQ55aI&WT(hDzr7+xDU^)f5t$j5wQlTj z@*y1FAVYRp=U7M{8NYKy>K(aykdiC*^dxsJCd9N}($n~OP28!INxBDz=H?DsN!HT( zwXDZOK9F(zQ>$E!P7k~vz#h~;i!2>)r z{iooGgn%AeUy9h`pjR!9wy&7&k*-BNYotqCb`NX~r-~vIYO&P6OI2K}7QOj>|M9o< zz-c1r#IveeT)V@+@xZDU+)Ky*i<<300{P-=@1a_}xyzH$;a!V@Vp95P=3Wd!VWnWX z7E{x5k98AkVaWWgmElM)^zW@k(bM?%r2Im}8qME{t2bJgFZM!0>rjfbSRDqR4ogyK z{&=Yq1KZ7ep@XbqnX$deIBT{F?Yqi#nW&=kTMj*Fk$@e4s3; z7h!_MUkYE;;l6Zj+#k_8WM9?(?{Zx)=!S)_R3_A+Db8V1PqhwV5lse6)4kv@`_3>| zT?eb<3y*v6)Io#kw)gzzJ`i^t5S1UT1NGs-NGH!ae0P34w|cM-s^Tka#lP$D-e1{D zA-)bT;xkJByU>Th`}8qaxaw*9Z{2^K#@W}5yZ8fd_hJ4Y{rSDp_2|{t8@nM};esYi558RH+#dT1C}8q6~F!{MjoKJnIi=-m_{ zH(J-z`bajCE(!M|sQ*{-*jzns&)rI>^s9$ga__HB)qaFAgv1lJQV=B;-`Af}4;Ra6 zFR%Okc)(u7nsbx_Ru`++T#fbEVt32iJfI)WUGgHGmnb;V<-YoCrXG!KR8gD!e*8?+ z7v;5}z)EgHiH(s$tGgZ`q<8nDzt>bTz?*_ij8*Bif)uQr3sYuY?#K69#zf&%3bgrx zHR8@wVB@ZF^c~j#EkAtFd4WQ~n`@l+l1wOI*9v@TaB2Y0&Wx-MeWZX=Dk=UX7~c9;=ApzO{`g+(=xu6%VNP4Kg>eH^ z=lqSs^#@Uc?VW7%4d7;vFCBQ^fbUr`{}x;ZX?fn=sx~z4X6CGBvx#rO^BDmxySPCR z3|MXdiBgdk8o5b|+KvvmecH8VA?$k`hy|bZW zhCb5oBs~?+dUaCwZy&;qTTVN(-%`QLH0S4hh>F?kOzXyyLuk2e!SJb=iqV_8zo#!! z!O?p~sM}x&j!Hy=`XCirV!_NV7F5Uu=o{HT8-k%-k;dd26{osV7f<=p-VfW=V-x=Nj>pBW2gwHY%;2gqz^HG=vAAQ)Uuw zH^Rtis>bar6}fg_oTv5-<9d6_of+>&Sj@z5e`RmP5l8zc&ho=Z>K~rHn%M|d>Tx|^ zsYY5~MfkFU*)TkLrP6=3(K;DRg|juUHp0ro?~t6|Fl?DV57p8*{osWCmGdr*2$1Qy zXjU+c^&I=s%H2%}QI+yCiD^X3ZiJJEhOtv9{_zL7CKR0Qb;uz#LiEuueQmlCP{y6F zox0hC4?#>f4}WL`GrL(Gr|<~;>GcDvJ(|!%-{{cF)P$A5q~^_{{c1Ac_07) literal 0 HcmV?d00001 diff --git a/src/blackbody_common.c b/src/blackbody_common.c new file mode 100644 index 0000000..30fb3d2 --- /dev/null +++ b/src/blackbody_common.c @@ -0,0 +1,128 @@ +#include "blackbody_internal.h" + +#include + +void blackbody_rayleigh_jeans_xyz(double temperature_K, int terms, + double xyz[3]) { + /* CSV-derived coefficients multiply T, 1, and T^-1. */ + static const double coefficients[3][3] = { + {9.75239167492079072e3, -1.3438558214879372e8, + 6.31051886434512943e11}, + {9.51500076632225143e3, -1.25823411343233966e8, + 5.58059336995288075e11}, + {2.13884044303540753e4, -3.43402928472968325e8, + 1.84249991169126495e12}}; + xyz[0] = xyz[1] = xyz[2] = 0.0; + if (!isfinite(temperature_K) || temperature_K <= 0.0 || terms <= 0) + return; + if (terms > 3) terms = 3; + for (int channel = 0; channel < 3; ++channel) { + xyz[channel] = coefficients[channel][0] * temperature_K; + if (terms >= 2) xyz[channel] += coefficients[channel][1]; + if (terms >= 3) xyz[channel] += coefficients[channel][2] / temperature_K; + } +} + +static double chebyshev_clenshaw(const double coefficients[9], double x) { + double b1 = 0.0, b2 = 0.0; + for (int degree = 8; degree >= 1; --degree) { + const double b0 = 2.0 * x * b1 - b2 + coefficients[degree]; + b2 = b1; + b1 = b0; + } + return coefficients[0] + x * b1 - b2; +} + +int blackbody_low_hybrid_log_xyz(double temperature_K, double log_xyz[3]) { + static const double crossover_K = 4.09732109813541501; + static const double maximum_K = 670.146556; + static const double q_over_endpoint[3] = { + 17334.66117474619, 17334.66117474619, 22135.028884675903}; + static const double endpoint_delta[3] = { + 6.27453428003507278e-6, 6.27523604457500139e-6, + 9.72024116337713046e-5}; + /* The 1 nm-linear endpoint makes these polynomials exactly cubic. Z has + * zero constant term because its 649--650 nm segment ends at zero. */ + static const double endpoint[3][4] = { + {1.81136495194567529e-2, 6.61028930839653126e-5, + 1.48916636199141013e-8, 1.27817061855715528e-12}, + {6.54117161006302892e-3, 2.38707969128581198e-5, + 5.37761358258185981e-9, 4.61567454126977889e-13}, + {0.0, 1.18022769903543601e-3, + 2.13277824065096828e-7, 1.44529622149769957e-11}}; + static const double edges[7] = { + 4.09732109813541501, 9.58175285991856125, 22.4073207028704393, + 52.4004353296953668, 122.540559808645089, 286.566107776330601, + 670.146556}; + static const double coefficients[6][3][9] = { + { + {-2.19805342404174431,-0.427989370548201992,0.0457830231212280024,-0.00653024253368962181,0.00104792258851775496,-0.000179376495463268954,0.0000320194392828742251,-6.07285820395394865e-6,1.09768907792163728e-6}, + {-3.21660305217802182,-0.427989317567608284,0.0457830153418609421,-0.00653024167298543334,0.00104792254862916667,-0.000179376508222333492,0.0000320194445897039005,-6.07285961285998572e-6,1.09768936710106915e-6}, + {-3.15707121695679764,-0.842260091928597102,0.0900358213706993461,-0.0130901684838339702,0.00216325774923957118,-0.000378541811375747837,0.0000677130171872672182,-0.0000125764470117522098,2.20543446698214866e-6} + },{ + {-1.31770949444839666,-0.441299064080369115,0.0486743388559153732,-0.00715733212557194619,0.00118377405373184402,-0.000208808548120484381,0.0000384098393642844184,-7.51938342275489625e-6,1.39706879054438949e-6}, + {-2.33625922712095605,-0.441299020893064677,0.0486743393223627113,-0.00715733306978699914,0.00118377427098721058,-0.000208808583754120261,0.0000384098442083096037,-7.51938398617148892e-6,1.39706884379154304e-6}, + {-1.41692623695155927,-0.878068580096058675,0.0982372095233778281,-0.0142171838603874451,0.00222410258776171213,-0.000361979109489253779,0.000060753034151493841,-0.0000107563963609356147,1.83055417864655942e-6} + },{ + {-0.3921912858481657,-0.475886418801351195,0.0565482685596524812,-0.00894861708837536552,0.00159121925979185938,-0.000301479141093084493,0.0000595321929174356658,-0.0000125546605490984269,2.48503406235378882e-6}, + {-1.41074105985953395,-0.475886424111238694,0.0565482743959804094,-0.00894861785536701427,0.00159121932381919629,-0.000301479143653568406,0.0000595321922909604336,-0.0000125546602334654497,2.48503397406341601e-6}, + {0.405784497738078867,-0.915430672087484161,0.0971591394650644897,-0.0129050910771596898,0.00191184765152995853,-0.000308981746742567431,0.0000538107592092635484,-0.0000103283467419743846,1.93585445253475723e-6} + },{ + {0.661994617272024548,-0.581350758628830698,0.0836597467165894418,-0.0158988753362009282,0.00336718167075529186,-0.000753662101271342622,0.000174585741711791468,-0.0000433048808785784171,9.75843572764750381e-6}, + {-0.356555122545496873,-0.581350782933150466,0.0836597486771623816,-0.015898875146527317,0.00336718160135198675,-0.000753662094742459443,0.000174585746021139991,-0.0000433048853568368614,9.75843776186480067e-6}, + {2.25724030901313405,-0.906149909784837897,0.0947967083167351412,-0.0147236528292987198,0.00294177345982300343,-0.000664178311762264752,0.000153812749558811817,-0.0000358188260289889306,7.26093914036807359e-6} + },{ + {2.24517368371397051,-1.10137638803233863,0.260109466476170278,-0.0728062330999542422,0.0205487830596007146,-0.00547019170643925872,0.00127956756145661419,-0.000220004358600067039,0.0000138071701243527944}, + {1.22677135347875279,-1.10165295376793399,0.260338122908468991,-0.0729732322401341771,0.0206568819007716045,-0.00553248986531515299,0.00131198334055904886,-0.000236511389192968018,0.0000196226106014253749}, + {4.22881792094750399,-1.09273677796148737,0.167321326581346532,-0.0312233937525439919,0.00516074702621717587,-0.000696731180600597722,0.000080598581477044639,-0.000015356852062622951,5.09084046909121171e-6} + },{ + {6.33932556236418281,-3.07589360106486531,0.456226630959950593,-0.00444896560449986224,-0.00803034195544535392,-0.000626248357061336987,0.000421665206117149112,0.0000469192966084690651,-0.0000234956663324828471}, + {5.3740651556541423,-3.15361482486290811,0.488370588919019434,-0.0119565427661776158,-0.0071239164324135692,-0.000725343801153405533,0.000461709147315171223,0.0000420002316640812232,-0.0000256745877581467219}, + {6.97759542644704118,-1.7029539661084809,0.314329784196094464,-0.0805402042898634484,0.0214126542337208624,-0.00404920923437737251,0.0000432174624337086157,0.000376147060473595442,-0.00015590604333367762} + }}; + + if (!isfinite(temperature_K) || temperature_K <= 0.0 || + temperature_K > maximum_K) + return -1; + if (temperature_K <= crossover_K) { + for (int channel = 0; channel < 3; ++channel) { + double polynomial = endpoint[channel][3]; + for (int term = 2; term >= 0; --term) + polynomial = endpoint[channel][term] + temperature_K * polynomial; + log_xyz[channel] = -q_over_endpoint[channel] / temperature_K + + log(temperature_K * polynomial) + + endpoint_delta[channel] * temperature_K / crossover_K; + } + return 0; + } + + int segment = 0; + while (segment < 5 && temperature_K >= edges[segment + 1]) ++segment; + const double inverse_center = + 0.5 * (1.0 / edges[segment] + 1.0 / edges[segment + 1]); + const double inverse_half = + 0.5 * (1.0 / edges[segment] - 1.0 / edges[segment + 1]); + const double x = (1.0 / temperature_K - inverse_center) / inverse_half; + for (int channel = 0; channel < 3; ++channel) + log_xyz[channel] = -q_over_endpoint[channel] / temperature_K + + chebyshev_clenshaw(coefficients[segment][channel], x); + return 0; +} + +int blackbody_low_hybrid_xyz(double temperature_K, double xyz[3]) { + double log_xyz[3]; + if (blackbody_low_hybrid_log_xyz(temperature_K, log_xyz)) { + xyz[0] = xyz[1] = xyz[2] = 0.0; + return -1; + } + for (int channel = 0; channel < 3; ++channel) + xyz[channel] = exp(log_xyz[channel]); + return 0; +} + +LinearRgb blackbody_xyz_to_linear_rgb(const double xyz[3]) { + return (LinearRgb){ + fmax(0.0, 3.24096994 * xyz[0] - 1.53738318 * xyz[1] - 0.49861076 * xyz[2]), + fmax(0.0, -0.96924364 * xyz[0] + 1.87596750 * xyz[1] + 0.04155506 * xyz[2]), + fmax(0.0, 0.05563008 * xyz[0] - 0.20397696 * xyz[1] + 1.05697151 * xyz[2])}; +} diff --git a/src/blackbody_internal.h b/src/blackbody_internal.h new file mode 100644 index 0000000..6267194 --- /dev/null +++ b/src/blackbody_internal.h @@ -0,0 +1,12 @@ +#ifndef BLACKBODY_INTERNAL_H +#define BLACKBODY_INTERNAL_H + +#include "optics.h" + +void blackbody_rayleigh_jeans_xyz(double temperature_K, int terms, + double xyz[3]); +int blackbody_low_hybrid_log_xyz(double temperature_K, double log_xyz[3]); +int blackbody_low_hybrid_xyz(double temperature_K, double xyz[3]); +LinearRgb blackbody_xyz_to_linear_rgb(const double xyz[3]); + +#endif diff --git a/src/blackbody_lut.c b/src/blackbody_lut.c new file mode 100644 index 0000000..c984a18 --- /dev/null +++ b/src/blackbody_lut.c @@ -0,0 +1,175 @@ +#include "blackbody_internal.h" + +#include +#include +#include +#include +#include +#include + +enum { BLACKBODY_LUT_VERSION_CIE1931 = 3, + BLACKBODY_LUT_ELEMENT_F64 = 1, + BLACKBODY_LUT_REFERENCE_CIE1931_1NM_LINEAR = 3 }; + +static const char default_table_path[] = + "assets/blackbody/cie1931_2deg_xyz_1024.grbblut"; + +typedef struct { + char magic[8]; + uint32_t version, header_bytes, endian_marker, element_type; + uint32_t reference_type, wavelength_min_nm, wavelength_max_nm, + wavelength_step_nm; + uint64_t node_count; + double log_temperature_min, log_temperature_max; + double planck_h, planck_c, planck_k; + char generator[64]; + uint64_t payload_checksum; +} BlackbodyLutHeader; + +_Static_assert(sizeof(BlackbodyLutHeader) == 160, + "blackbody LUT header layout changed"); + +typedef struct { + double *xyz; + size_t nodes; + double log_min, log_max, index_scale; + int ready; +} BlackbodyLut; + +static BlackbodyLut table; + +static uint64_t fnv1a64(const void *data, size_t bytes) { + const unsigned char *p = data; + uint64_t hash = UINT64_C(14695981039346656037); + for (size_t i = 0; i < bytes; ++i) { + hash ^= p[i]; + hash *= UINT64_C(1099511628211); + } + return hash; +} + +static int validate_header(const BlackbodyLutHeader *h) { + return memcmp(h->magic, "GRBBLUT3", 8) == 0 && + h->version == BLACKBODY_LUT_VERSION_CIE1931 && + h->header_bytes == sizeof *h && + h->endian_marker == UINT32_C(0x01020304) && + h->element_type == BLACKBODY_LUT_ELEMENT_F64 && + h->reference_type == BLACKBODY_LUT_REFERENCE_CIE1931_1NM_LINEAR && + h->wavelength_min_nm == 360 && h->wavelength_max_nm == 830 && + h->wavelength_step_nm == 1 && h->node_count >= 4 && + h->node_count <= SIZE_MAX / (3 * sizeof(double)) && + isfinite(h->log_temperature_min) && + isfinite(h->log_temperature_max) && + h->log_temperature_max > h->log_temperature_min && + h->planck_h == 6.62607015e-34 && h->planck_c == 299792458.0 && + h->planck_k == 1.380649e-23; +} + +static int read_table(const char *path) { + FILE *file = fopen(path, "rb"); + BlackbodyLutHeader h; + if (file == NULL || fread(&h, sizeof h, 1, file) != 1 || + !validate_header(&h)) { + if (file != NULL) fclose(file); + return -1; + } + double *xyz = malloc(3 * (size_t)h.node_count * sizeof *xyz); + const size_t count = 3 * (size_t)h.node_count; + unsigned char extra; + if (xyz == NULL || fread(xyz, sizeof *xyz, count, file) != count || + fread(&extra, 1, 1, file) != 0 || ferror(file) || + fnv1a64(xyz, count * sizeof *xyz) != h.payload_checksum) { + free(xyz); fclose(file); return -1; + } + fclose(file); + table = (BlackbodyLut){.xyz = xyz, .nodes = (size_t)h.node_count, + .log_min = h.log_temperature_min, + .log_max = h.log_temperature_max, .ready = 1}; + table.index_scale = (table.nodes - 1) / (table.log_max - table.log_min); + return 0; +} + +int blackbody_backend_init(const char *table_path, size_t lut_nodes, + double temperature_min_K, + double temperature_max_K, + const char *write_table_path, FILE *report) { + struct timespec start, finish; + timespec_get(&start, TIME_UTC); + if (table.ready) return -1; + if (lut_nodes != 0 || isfinite(temperature_min_K) || + isfinite(temperature_max_K) || write_table_path != NULL) + return -1; + const char *path = table_path != NULL ? table_path : default_table_path; + if (read_table(path)) return -1; + timespec_get(&finish, TIME_UTC); + if (report != NULL) + fprintf(report, + "Blackbody LUT: %zu CIE 1931 2-deg 1 nm-linear XYZ nodes, " + "T=[%.9g, %.9g] K, loaded %s in %.6f s; " + "payload fnv1a64=%016llx\n", + table.nodes, + exp(table.log_min), exp(table.log_max), + path, + finish.tv_sec - start.tv_sec + + 1e-9 * (finish.tv_nsec - start.tv_nsec), + (unsigned long long)fnv1a64(table.xyz, + 3 * table.nodes * sizeof *table.xyz)); + return 0; +} + +void blackbody_backend_destroy(void) { + free(table.xyz); + table = (BlackbodyLut){0}; +} + +const char *blackbody_backend_name(void) { return "lut"; } + +LinearRgb blackbody_to_linear_rgb(double temperature_K) { + if (!table.ready || !isfinite(temperature_K) || temperature_K <= 0.0) { + return (LinearRgb){0}; + } + const double log_t = log(temperature_K); + if (log_t < table.log_min) { + double xyz[3]; + if (blackbody_low_hybrid_xyz(temperature_K, xyz)) + return (LinearRgb){0}; + return blackbody_xyz_to_linear_rgb(xyz); + } + if (log_t > table.log_max) { + double xyz[3]; + blackbody_rayleigh_jeans_xyz(temperature_K, 3, xyz); + return blackbody_xyz_to_linear_rgb(xyz); + } + const double position = (log_t - table.log_min) * table.index_scale; + size_t low = (size_t)position; + if (low >= table.nodes - 1) low = table.nodes - 2; + const double fraction = position - low; + double xyz[3]; + { + size_t base; + double x; + if (low == 0) { + base = 0; + x = fraction; + } else if (low >= table.nodes - 2) { + base = table.nodes - 4; + x = 2.0 + fraction; + } else { + base = low - 1; + x = 1.0 + fraction; + } + const double weights[4] = { + -(x - 1.0) * (x - 2.0) * (x - 3.0) / 6.0, + x * (x - 2.0) * (x - 3.0) / 2.0, + -x * (x - 1.0) * (x - 3.0) / 2.0, + x * (x - 1.0) * (x - 2.0) / 6.0}; + for (int channel = 0; channel < 3; ++channel) { + xyz[channel] = 0.0; + for (int point = 0; point < 4; ++point) + xyz[channel] += weights[point] * + table.xyz[3 * (base + (size_t)point) + channel]; + xyz[channel] = fmax(0.0, xyz[channel]); + } + } + return blackbody_xyz_to_linear_rgb(xyz); +} diff --git a/src/main.c b/src/main.c index 87711db..13a573d 100644 --- a/src/main.c +++ b/src/main.c @@ -12,6 +12,7 @@ #include #include #include +#include #include #include #include @@ -52,6 +53,7 @@ typedef struct { double slab_duration; double minkowski_proper_acceleration; int catalog_load_workers; + const char *blackbody_table_path; RefinementConfig refinement; } Settings; @@ -319,6 +321,8 @@ static int parse_args(int argc, char **argv, Settings *s, &s->minkowski_proper_acceleration)) { } else if (!strcmp(argv[i], "--catalog-load-workers") && i + 1 < argc && !parse_int(argv[++i], &s->catalog_load_workers)) { + } else if (!strcmp(argv[i], "--blackbody-table") && i + 1 < argc) { + s->blackbody_table_path = argv[++i]; } else if (!strcmp(argv[i], "--write-minkowski-accel-track") && i + 1 < argc) s->write_minkowski_accel_track_path = argv[++i]; @@ -377,6 +381,10 @@ static void print_help(const char *program) { " --psf-min-y Y Skip events below this linear HDR luminance (default: 0, disabled)\n" " --psf-direct Disable the PSF lookup cache (default: disabled)\n" " --catalog-load-workers N All-sky catalog loader workers (default: 4)\n" + " --blackbody-table FILE Explicit GRBBLUT3 table\n" + " (default: assets/blackbody/cie1931_2deg_xyz_1024.grbblut)\n", + stdout); + fputs( "\nAdaptive lens mesh:\n" " --coarse-cell-pixels N Initial mesh cell size in pixels (default: 16)\n" " --refine-max-level N Maximum refinement level (default: 0)\n" @@ -1111,6 +1119,18 @@ int main(int argc, char **argv) { } fprintf(stderr, "Created test catalog: %s\n", settings.catalog_path); } + if (blackbody_backend_init(settings.blackbody_table_path, + 0, NAN, NAN, NULL, stderr)) { + fprintf(stderr, + "Blackbody backend '%s' initialization failed; the LUT loads the " + "repository CIE GRBBLUT3 table by default or a valid explicit " + "table.\n", + blackbody_backend_name()); + catalog_destroy(&catalog); + spacetime_destroy(&spacetime); + return 1; + } + fprintf(stderr, "Blackbody backend: %s\n", blackbody_backend_name()); if (!settings.psf_direct && psf_kernel_cache_init(&settings.psf_cache, &settings.psf, settings.psf_relative_tail)) { #ifdef PSF_BACKEND_DUMMY @@ -1118,6 +1138,7 @@ int main(int argc, char **argv) { stderr); catalog_destroy(&catalog); spacetime_destroy(&spacetime); + blackbody_backend_destroy(); return 1; #else fputs("PSF cache construction failed; using direct evaluator.\n", stderr); @@ -1128,6 +1149,7 @@ int main(int argc, char **argv) { const int result = render_lens_map(&settings, &catalog); catalog_destroy(&catalog); psf_kernel_cache_destroy(&settings.psf_cache); + blackbody_backend_destroy(); return result == 0 ? 0 : 1; } int result = settings.frames_dir != NULL @@ -1137,5 +1159,6 @@ int main(int argc, char **argv) { spacetime_destroy(&spacetime); catalog_destroy(&catalog); psf_kernel_cache_destroy(&settings.psf_cache); + blackbody_backend_destroy(); return result == 0 ? 0 : 1; } diff --git a/src/optics.c b/src/optics.c index 75bdf3c..08a6970 100644 --- a/src/optics.c +++ b/src/optics.c @@ -21,68 +21,6 @@ static double clamp(double value, double low, double high) return value < low ? low : value > high ? high : value; } -/* Wyman, Sloan, and Shirley (2013), Eq. 4: analytic fits to the CIE 1931 - * 2-degree color-matching functions. Wavelength is in nanometres. */ -static void cie_1931_2deg(double wavelength_nm, double xyz[3]) -{ - const double x1 = (wavelength_nm - 442.0) * - (wavelength_nm < 442.0 ? 0.0624 : 0.0374); - const double x2 = (wavelength_nm - 599.8) * - (wavelength_nm < 599.8 ? 0.0264 : 0.0323); - const double x3 = (wavelength_nm - 501.1) * - (wavelength_nm < 501.1 ? 0.0490 : 0.0382); - const double y1 = (wavelength_nm - 568.8) * - (wavelength_nm < 568.8 ? 0.0213 : 0.0247); - const double y2 = (wavelength_nm - 530.9) * - (wavelength_nm < 530.9 ? 0.0613 : 0.0322); - const double z1 = (wavelength_nm - 437.0) * - (wavelength_nm < 437.0 ? 0.0845 : 0.0278); - const double z2 = (wavelength_nm - 459.0) * - (wavelength_nm < 459.0 ? 0.0385 : 0.0725); - xyz[0] = 0.362 * exp(-0.5 * x1 * x1) + - 1.056 * exp(-0.5 * x2 * x2) - 0.065 * exp(-0.5 * x3 * x3); - xyz[1] = 0.821 * exp(-0.5 * y1 * y1) + - 0.286 * exp(-0.5 * y2 * y2); - xyz[2] = 1.217 * exp(-0.5 * z1 * z1) + - 0.681 * exp(-0.5 * z2 * z2); -} - -static double planck_radiance_lambda(double wavelength_m, double temperature_K) -{ - const double h = 6.62607015e-34; - const double c = 299792458.0; - const double k = 1.380649e-23; - const double exponent = h * c / (wavelength_m * k * temperature_K); - return 2.0 * h * c * c / - (pow(wavelength_m, 5.0) * expm1(exponent)); -} - -LinearRgb blackbody_to_linear_rgb(double temperature_K) -{ - /* Integrate Planck spectral radiance from 380 to 780 nm into CIE XYZ, - * then transform XYZ to linear sRGB. Results are W m^-2 sr^-1 in each - * linear-primary channel, before the catalog amplitude and exposure. */ - double xyz[3] = {0.0, 0.0, 0.0}; - const double wavelength_step_m = 5e-9; - if (!isfinite(temperature_K) || temperature_K <= 0.0) - return (LinearRgb){0.0, 0.0, 0.0}; - for (int wavelength_nm = 380; wavelength_nm <= 780; wavelength_nm += 5) { - double matching[3]; - const double radiance = - planck_radiance_lambda(wavelength_nm * 1e-9, temperature_K); - cie_1931_2deg(wavelength_nm, matching); - for (int channel = 0; channel < 3; ++channel) - xyz[channel] += radiance * matching[channel] * wavelength_step_m; - } - return (LinearRgb){ - clamp(3.24096994 * xyz[0] - 1.53738318 * xyz[1] - 0.49861076 * xyz[2], - 0.0, INFINITY), - clamp(-0.96924364 * xyz[0] + 1.87596750 * xyz[1] + 0.04155506 * xyz[2], - 0.0, INFINITY), - clamp(0.05563008 * xyz[0] - 0.20397696 * xyz[1] + 1.05697151 * xyz[2], - 0.0, INFINITY)}; -} - enum { PSF_PHASE_RESOLUTION = 64, PSF_QUADRATURE_ORDER = 4, diff --git a/src/optics.h b/src/optics.h index b511bae..7f1a6bf 100644 --- a/src/optics.h +++ b/src/optics.h @@ -50,6 +50,16 @@ extern "C" { /* Integrate a Planck spectrum into absolute linear-sRGB spectral radiance * (W m^-2 sr^-1), before catalog amplitude and display exposure. */ LinearRgb blackbody_to_linear_rgb(double temperature_K); +/* Backend initialization happens before renderer workers start. Integral and + * fixed backends reject LUT-only options. The LUT backend loads the sole + * repository CIE table by default, or validates an explicit equivalent file; + * legacy runtime-generation arguments are rejected. */ +int blackbody_backend_init(const char *table_path, size_t lut_nodes, + double temperature_min_K, + double temperature_max_K, + const char *write_table_path, FILE *report); +void blackbody_backend_destroy(void); +const char *blackbody_backend_name(void); int psf_kernel_cache_init(PsfKernelCache *cache, const PointSpreadFunction *psf, double relative_tail_fraction); diff --git a/tests/capture_psf.c b/tests/capture_psf.c index ef9cff2..48e6edd 100644 --- a/tests/capture_psf.c +++ b/tests/capture_psf.c @@ -1,6 +1,7 @@ /* Standalone bounded diagnostic. Include the production producer so its * private query/mapping path stays identical; renderer builds omit the hook. */ #include "optics.h" +#include #include #include #include @@ -100,6 +101,7 @@ int main(int argc, char **argv) { if (lens_map_read(argv[1],&map) || map.frame_count!=1 || last>map.frames[0].mesh.triangle_count || catalog_load_csv(&catalog,argv[2]) || catalog.count>32768) return 2; + if (blackbody_backend_init(NULL, 0, NAN, NAN, NULL, stderr)) return 1; const int indexed=!strcmp(argv[6],"--diagnose-indexed") || !strcmp(argv[6],"--diagnose-indexed-serial"); if(indexed && index_subset(&catalog))return 1; const char *kind=indexed ? "indexed-subset" : "memory"; @@ -126,7 +128,7 @@ int main(int argc, char **argv) { } printf("DIAGNOSTIC ONLY, no HDR: workers=%d catalog=%s stars=%zu triangles=[%zu,%zu) events=%zu cached=%zu wing=%zu direct=%zu discarded=%zu checksum=%016llx producer_wall=%.9f capped_or_failed=%d\n", workers,kind,catalog.count,first,last,total_events,totals[0]+totals[2],totals[2],totals[1],totals[3],(unsigned long long)total_checksum,omp_get_wtime()-start,failed); - psf_kernel_cache_destroy(&cache);catalog_destroy(&catalog);lens_map_destroy(&map); + psf_kernel_cache_destroy(&cache);catalog_destroy(&catalog);lens_map_destroy(&map);blackbody_backend_destroy(); return failed ? 1 : 0; } PsfEventSink sink={.cache=&cache}; @@ -146,6 +148,6 @@ int main(int argc, char **argv) { } if (ferror(f) || fclose(f)) return 1; } - psf_kernel_cache_destroy(&cache);catalog_destroy(&catalog);lens_map_destroy(&map); + psf_kernel_cache_destroy(&cache);catalog_destroy(&catalog);lens_map_destroy(&map);blackbody_backend_destroy(); return 0; } diff --git a/tests/test_frame.c b/tests/test_frame.c index ccd08b3..832f5e3 100644 --- a/tests/test_frame.c +++ b/tests/test_frame.c @@ -74,7 +74,8 @@ int main(void) { FrameLensMesh mesh = {0}; double *hdr = calloc((size_t)width * height * 3, sizeof *hdr); int result = 1; - if (hdr == NULL || spacetime_create_minkowski(&spacetime, 10.0) || + if (blackbody_backend_init(NULL, 0, NAN, NAN, NULL, stderr) || + hdr == NULL || spacetime_create_minkowski(&spacetime, 10.0) || frame_lens_mesh_build_coarse(&mesh, width, height, 20, 30.0) || frame_lens_mesh_trace(&mesh, &spacetime, &observer, &trace)) goto done; @@ -507,6 +508,7 @@ int main(void) { done: frame_lens_mesh_destroy(&mesh); spacetime_destroy(&spacetime); + blackbody_backend_destroy(); free(hdr); return result; }