nx_ice_distribution.nx source
↩ module page · 355 lines · 14028 B
1// nx_ice_distribution.nx -- THE CRYSTAL SIZE DISTRIBUTION, and the coupling
2// that nx_ice_recrystal had to assert because it did not have one.
3//
4// ===== WHY THIS EXISTS ============================================
5//
6// nx_ice_recrystal models cycling damage as
7//
8// V_new = V_old * (1 + phi * G)
9//
10// where phi is the fraction of ice MASS that melts and refreezes on an
11// excursion -- derived honestly from the freezing curve -- and G is a
12// geometric coupling that was ASSERTED. That constant was flagged as the
13// organ's weak point, and it deserved the flag: it converts "this much mass
14// moved" into "the mean diameter grew this much", which is a statement about
15// the crystal size distribution, and the organ had no distribution.
16//
17// ★IT IS NOT A FREE PARAMETER. Given a distribution the whole step is
18// exact arithmetic:
19//
20// 1. Gibbs-Thomson orders the melting -- the smallest crystals have the
21// highest surface curvature, the most depressed local melting point,
22// and therefore go first. The order is physics, not a modelling choice.
23// 2. Remove the smallest crystals until their cumulative MASS reaches
24// phi * M. Count how many crystals that took.
25// 3. That mass refreezes onto the survivors, so total mass is unchanged
26// while the count has fallen.
27// 4. Mean volume is mass over count, so it rises by exactly N / N'.
28//
29// ★★THE WHOLE POINT IS THAT THE NUMBER FRACTION IS NOT THE MASS FRACTION.
30// Small crystals are numerous and light. Removing a given fraction of the
31// MASS removes a much larger fraction of the COUNT, because it takes many
32// small crystals to make up that mass. A model that quietly equates the two
33// -- which is what a coupling of 1.0 does -- therefore UNDER-STATES how much
34// the mean size grows per cycle. Under-stating damage is the direction that
35// says a product is fine when it is not.
36//
37// So this organ computes the coupling instead of asserting it, and
38// dist_coupling_error_permil measures how wrong the assertion was.
39//
40// ===== WHAT REMAINS ASSERTED ======================================
41//
42// The DISTRIBUTION ITSELF. A real one comes from image analysis of a real
43// sample; the built-in reference profile is a plausible shape for well-made
44// ice cream and is flagged by dist_profile_is_reference(). What this organ
45// removes is the coupling constant, not the need for measurement -- callers
46// with real bin counts get a fully derived answer, and callers using the
47// reference profile get a shape-correct one that says so.
48//
49// grounded: gibbs_thomson_curvature + lifshitz_slyozov_wagner
50// + donhowe_hartel_ice_cream_crystal_distributions
51// genealogy_id: nishi_food_chem_lane_ice_distribution_2026_07
52// license_tier: ORIGINAL
53
54import "nx_syscalls.nx"
55const DIST_MAGIC_1200: i64 = 1200
56const DIST_MAGIC_2100: i64 = 2100
57const DIST_MAGIC_2600: i64 = 2600
58const DIST_MAGIC_2200: i64 = 2200
59const DIST_MAGIC_1300: i64 = 1300
60const DIST_MAGIC_3000000: i64 = 3000000
61
62const DIST_INVALID: i64 = 0 - 1
63
64// Diameters in TENTHS OF A MICROMETRE, matching nx_ice_recrystal so the two
65// organs never need a unit conversion between them.
66const DIST_MAXBINS: i64 = 32
67
68// The reference profile: a plausible well-made ice cream, mean near 25 um,
69// with the long small-crystal tail that real samples show. ASSERTED shape.
70const DIST_REF_BINS: i64 = 8
71
72struct IceDist {
73 n: i64,
74 d: *i64, // bin diameter (0.1 um), ascending
75 cnt: *i64, // crystals in that bin
76}
77
78func dist_new(nbins: i64) -> *IceDist {
79 if nbins <= 0 { return 0 as *IceDist }
80 if nbins > DIST_MAXBINS { return 0 as *IceDist }
81 let s: *IceDist = sys_mmap(64) as *IceDist
82 s.n = nbins
83 s.d = sys_mmap(DIST_MAXBINS * 8) as *i64
84 s.cnt = sys_mmap(DIST_MAXBINS * 8) as *i64
85 var i: i64 = 0
86 while i < nbins { s.d[i] = 0; s.cnt[i] = 0; i = i + 1 }
87 return s
88}
89
90func dist_set(s: *IceDist, i: i64, diameter: i64, count: i64) -> i64 {
91 if i < 0 { return DIST_INVALID }
92 if i >= s.n { return DIST_INVALID }
93 if diameter <= 0 { return DIST_INVALID }
94 if count < 0 { return DIST_INVALID }
95 s.d[i] = diameter
96 s.cnt[i] = count
97 return 0
98}
99
100// The built-in reference profile. Counts fall off toward the large end,
101// which is what gives the small bins their disproportionate share of the
102// COUNT and their small share of the MASS -- the whole effect this organ
103// exists to quantify.
104func dist_reference() -> *IceDist {
105 let s: *IceDist = dist_new(DIST_REF_BINS)
106 dist_set(s, 0, 100, DIST_MAGIC_1200)
107 dist_set(s, 1, 150, DIST_MAGIC_2100)
108 dist_set(s, 2, 200, DIST_MAGIC_2600)
109 dist_set(s, 3, 250, DIST_MAGIC_2200)
110 dist_set(s, 4, 300, DIST_MAGIC_1300)
111 dist_set(s, 5, 350, 620)
112 dist_set(s, 6, 400, 240)
113 dist_set(s, 7, 450, 80)
114 return s
115}
116
117func dist_profile_is_reference() -> i64 {
118 return 1
119}
120
121// Ascending order is a PRECONDITION of the melt walk, not a convention:
122// Gibbs-Thomson says the smallest go first, so a mis-ordered set of bins
123// would melt the wrong crystals and silently invert the result. Measured
124// rather than assumed, in the spirit of nx_peptide checking its bisection's
125// monotonicity.
126func dist_is_ascending(s: *IceDist) -> i64 {
127 var i: i64 = 1
128 while i < s.n {
129 if s.d[i] <= s.d[i - 1] { return 0 }
130 i = i + 1
131 }
132 return 1
133}
134
135func dist_total_count(s: *IceDist) -> i64 {
136 var t: i64 = 0
137 var i: i64 = 0
138 while i < s.n { t = t + s.cnt[i]; i = i + 1 }
139 return t
140}
141
142// Total mass in units of (0.1 um)^3 x crystals. Density is constant across
143// crystals so it cancels everywhere and is never introduced.
144func dist_total_mass(s: *IceDist) -> i64 {
145 var t: i64 = 0
146 var i: i64 = 0
147 while i < s.n {
148 let d: i64 = s.d[i]
149 t = t + s.cnt[i] * d * d * d
150 i = i + 1
151 }
152 return t
153}
154
155// Mean diameter, count-weighted (0.1 um).
156func dist_mean_diameter(s: *IceDist) -> i64 {
157 let n: i64 = dist_total_count(s)
158 if n <= 0 { return DIST_INVALID }
159 var t: i64 = 0
160 var i: i64 = 0
161 while i < s.n { t = t + s.cnt[i] * s.d[i]; i = i + 1 }
162 return t / n
163}
164
165// === The melt walk ================================================
166
167// Number of crystals lost when the smallest are melted until their
168// cumulative mass reaches phi (permil) of the total. Walks the bins in
169// ascending order, which is exactly the Gibbs-Thomson ordering.
170//
171// The final bin is taken FRACTIONALLY -- rounding a partial bin up to the
172// whole bin would over-count the crystals lost and inflate the very effect
173// this organ is measuring.
174func dist_count_melted(s: *IceDist, phi_permil: i64) -> i64 {
175 if phi_permil < 0 { return DIST_INVALID }
176 if phi_permil > 1000 { return DIST_INVALID }
177 if dist_is_ascending(s) == 0 { return DIST_INVALID }
178 let total: i64 = dist_total_mass(s)
179 if total <= 0 { return DIST_INVALID }
180 let target: i64 = total / 1000 * phi_permil
181 var acc: i64 = 0
182 var lost: i64 = 0
183 var i: i64 = 0
184 var done: i64 = 0
185 while i < s.n {
186 if done == 0 {
187 let d: i64 = s.d[i]
188 let binmass: i64 = s.cnt[i] * d * d * d
189 if acc + binmass <= target {
190 acc = acc + binmass
191 lost = lost + s.cnt[i]
192 } else {
193 let need: i64 = target - acc
194 if need > 0 {
195 if binmass > 0 { lost = lost + s.cnt[i] * need / binmass }
196 }
197 done = 1
198 }
199 }
200 i = i + 1
201 }
202 return lost
203}
204
205// Number fraction lost, in permil. ★This is the quantity the asserted
206// coupling was standing in for.
207func dist_number_fraction_permil(s: *IceDist, phi_permil: i64) -> i64 {
208 let lost: i64 = dist_count_melted(s, phi_permil)
209 if lost == DIST_INVALID { return DIST_INVALID }
210 let n: i64 = dist_total_count(s)
211 if n <= 0 { return DIST_INVALID }
212 return lost * 1000 / n
213}
214
215// === The coupling, derived ========================================
216
217// Mean VOLUME growth factor per cycle, in permil above 1.000.
218//
219// mass is conserved, count falls from N to N'
220// mean volume = M/N', so it grows by N/N' = 1/(1 - nu)
221//
222// Returned as (factor - 1) * 1000 so it composes in integer arithmetic.
223func dist_volume_growth_permil(s: *IceDist, phi_permil: i64) -> i64 {
224 let nu: i64 = dist_number_fraction_permil(s, phi_permil)
225 if nu == DIST_INVALID { return DIST_INVALID }
226 if nu >= 1000 { return DIST_INVALID }
227 return nu * 1000 / (1000 - nu)
228}
229
230// ★THE DERIVED COUPLING -- what nx_ice_recrystal's IR_COUPLING should be.
231//
232// its model: growth = phi * G / 1000
233// the truth: growth = nu / (1 - nu)
234// therefore: G = 1000 * nu / ((1 - nu) * phi)
235//
236// A coupling of 1000 (i.e. G = 1.0) is exactly the claim that the number
237// fraction equals the mass fraction. It does not.
238func dist_derived_coupling(s: *IceDist, phi_permil: i64) -> i64 {
239 if phi_permil <= 0 { return DIST_INVALID }
240 let growth: i64 = dist_volume_growth_permil(s, phi_permil)
241 if growth == DIST_INVALID { return DIST_INVALID }
242 return growth * 1000 / phi_permil
243}
244
245// How wrong was the assertion, in permil of the derived value? Signed:
246// negative means the asserted constant UNDER-states the damage, which is the
247// unsafe direction and is the answer this returns for any realistic profile.
248func dist_coupling_error_permil(s: *IceDist, phi_permil: i64, asserted: i64) -> i64 {
249 let derived: i64 = dist_derived_coupling(s, phi_permil)
250 if derived == DIST_INVALID { return DIST_INVALID }
251 if derived == 0 { return DIST_INVALID }
252 return (asserted - derived) * 1000 / derived
253}
254
255// Does the asserted constant err toward safety? 1 = it over-states damage
256// (conservative), 0 = it under-states it. A function rather than a comment,
257// because the direction of a modelling error is the part a consumer must not
258// have to infer.
259func dist_asserted_is_conservative(s: *IceDist, phi_permil: i64, asserted: i64) -> i64 {
260 let e: i64 = dist_coupling_error_permil(s, phi_permil, asserted)
261 if e == DIST_INVALID { return 0 }
262 if e >= 0 { return 1 }
263 return 0
264}
265
266// === Evolving the distribution ====================================
267
268// Apply one cycle: melt the smallest crystals accounting for phi of the mass
269// and redeposit onto the survivors, growing each surviving bin's diameter so
270// that total mass is conserved exactly.
271//
272// Redeposition is spread in proportion to existing bin mass, which is the
273// diffusion-limited case: a crystal already holding a larger share of the
274// mass presents more surface to grow on.
275func dist_apply_cycle(s: *IceDist, phi_permil: i64) -> i64 {
276 if dist_is_ascending(s) == 0 { return DIST_INVALID }
277 let total: i64 = dist_total_mass(s)
278 if total <= 0 { return DIST_INVALID }
279 if phi_permil <= 0 { return 0 }
280 if phi_permil > 1000 { return DIST_INVALID }
281 let target: i64 = total / 1000 * phi_permil
282 // Pass 1: remove the smallest crystals up to the mass target.
283 var acc: i64 = 0
284 var i: i64 = 0
285 var done: i64 = 0
286 while i < s.n {
287 if done == 0 {
288 let d: i64 = s.d[i]
289 let binmass: i64 = s.cnt[i] * d * d * d
290 if acc + binmass <= target {
291 acc = acc + binmass
292 s.cnt[i] = 0
293 } else {
294 let need: i64 = target - acc
295 if need > 0 {
296 if binmass > 0 {
297 let drop: i64 = s.cnt[i] * need / binmass
298 s.cnt[i] = s.cnt[i] - drop
299 acc = acc + drop * d * d * d
300 }
301 }
302 done = 1
303 }
304 }
305 i = i + 1
306 }
307 // Pass 2: grow the survivors so the removed mass is exactly restored.
308 let remaining: i64 = dist_total_mass(s)
309 if remaining <= 0 { return DIST_INVALID }
310 // Every surviving crystal's volume scales by total/remaining, so its
311 // diameter scales by the cube root of that -- applied in permil to keep
312 // the arithmetic integral.
313 // ⚠TWO REAL BUGS LIVED IN THIS LOOP AND THE MASS-CONSERVATION CHECK
314 // FOUND BOTH.
315 //
316 // (1) OVERFLOW. The first cut scaled each bin with
317 // v = v/remaining*total + (v%remaining)*total/remaining
318 // which is the standard trick for keeping a ratio exact. Here it is
319 // wrong: `remaining` is ~1e11, so v/remaining is 0 for every bin and
320 // the whole value goes through the second term as v*total. For the
321 // largest bin that is 9.1e7 x 1.5e11 = 1.35e19, PAST the i64 ceiling.
322 // The scale factor is a small number, so take the ratio FIRST and
323 // apply it in permil -- then no intermediate exceeds ~1e11.
324 //
325 // (2) TRUNCATION. The cube root floors, and flooring every surviving
326 // bin loses mass on every cycle -- a slow leak that would read as
327 // crystals mysteriously failing to grow. Rounded to NEAREST, the
328 // same law wave 3 banked when integer floor cost a whole PAC unit.
329 let scale_permil: i64 = total * 1000 / remaining
330 var j: i64 = 0
331 while j < s.n {
332 if s.cnt[j] > 0 {
333 let d0: i64 = s.d[j]
334 let v0: i64 = d0 * d0 * d0
335 let v: i64 = v0 / 1000 * scale_permil + (v0 % 1000) * scale_permil / 1000
336 var lo: i64 = 0
337 var hi: i64 = 1
338 while hi * hi * hi <= v {
339 lo = hi
340 hi = hi * 2
341 if hi > DIST_MAGIC_3000000 { return DIST_INVALID }
342 }
343 while lo + 1 < hi {
344 let mid: i64 = (lo + hi) / 2
345 if mid * mid * mid <= v { lo = mid } else { hi = mid }
346 }
347 // Round to nearest: keep lo+1 when it lands closer to v.
348 let below: i64 = v - lo * lo * lo
349 let above: i64 = (lo + 1) * (lo + 1) * (lo + 1) - v
350 if above < below { s.d[j] = lo + 1 } else { s.d[j] = lo }
351 }
352 j = j + 1
353 }
354 return 0
355}