nx_stability_rh.nx source
↩ module page · 285 lines · 11663 B
1// nx_stability_rh.nx -- the MOISTURE axis of stability, and the machinery to
2// SEPARATE it from the temperature axis.
3//
4// ===== WHY THIS EXISTS ============================================
5//
6// nx_stability recovers an activation energy from ICH's declaration that
7// 6 months at 40 C / 75% RH supports 24 months at 25 C / 60% RH. Those two
8// conditions differ in temperature AND in humidity, so what comes out is a
9// LUMPED effective parameter carrying the moisture contribution inside it.
10// nx_stability says so honestly but cannot do anything about it, because
11// two conditions differing in two variables are underdetermined -- one
12// equation, two unknowns, and no amount of arithmetic fixes that.
13//
14// This file adds the second axis and, with it, the ability to separate them
15// when the experiment actually supports separation.
16//
17// ===== THE MODEL ==================================================
18//
19// The moisture-modified Arrhenius used by accelerated-stability practice:
20//
21// ln k = ln A - Ea/(RT) + B*RH
22//
23// so between two conditions, in base ten (this suite's convention):
24//
25// log10(k2/k1) = E10*(1/T1 - 1/T2) + B10*(RH2 - RH1)
26//
27// where E10 = Ea/(R ln10) and B10 = B/ln10. Two unknowns, so two
28// independent contrasts -- three conditions -- determine both.
29//
30// ===== THE REFUSAL IS THE FEATURE =================================
31//
32// Three conditions are NECESSARY but not SUFFICIENT: they must not be
33// collinear in (1/T, RH) space. Three points on a line give a singular
34// system, and a solver that "handles" that by picking one of the infinitely
35// many solutions would emit a confident, specific, meaningless Ea. So the
36// determinant is computed and a zero REFUSES.
37//
38// srh_ich_pair_alone_is_separable() returns 0 as a FUNCTION, so this file
39// mechanically demonstrates the limitation nx_stability could only describe.
40//
41// ===== NO 1/T IS EVER MATERIALISED ================================
42//
43// Reciprocal temperature is held as 1/T x 10^7 per Kelvin, which keeps four
44// significant figures for any ambient temperature. This is not a stylistic
45// choice: nx_arrhenius computed 1/T in Q10 fixed point, truncated 0.00335 to
46// the integer 3, and therefore returned "temperature has no effect on rate"
47// for every input it was ever given. See the fix note in that file.
48//
49// grounded: ich_q1a_r2 + moisture_modified_arrhenius_accelerated_stability
50// genealogy_id: nishi_food_chem_lane_stability_rh_2026_07
51// license_tier: ORIGINAL
52
53import "nx_syscalls.nx"
54import "nx_pow10.nx"
55import "nx_stability.nx"
56const SRH_MAGIC_10000000: i64 = 10000000
57const SRH_MAGIC_10000: i64 = 10000
58const SRH_MAGIC_1000000: i64 = 1000000
59
60const SRH_INVALID: i64 = 0 - 1
61// A DISTINCT code from INVALID on purpose: the inputs are well-formed, the
62// experiment simply cannot answer the question. Collapsing the two would
63// let a caller treat a structural limitation as a bad argument and retry.
64const SRH_UNDERDETERMINED: i64 = 0 - 2
65
66// Reciprocal temperature is carried as 1/T x 10^7, in per-Kelvin.
67const SRH_INV_T_SCALE: i64 = 1000000000
68
69// Published moisture-sensitivity range for pharmaceutical degradation,
70// expressed as B10 (base-ten) in milli-decades per %RH. Used ONLY to sanity
71// check a recovered value, never as an input to one.
72const SRH_PUB_B10_LO_MILLI: i64 = 5
73const SRH_PUB_B10_HI_MILLI: i64 = 25
74
75// === Reciprocal temperature =======================================
76
77// 1/T x 10^7 per Kelvin, from whole Celsius. Four significant figures at
78// ambient (25 C -> 33540), which is what makes a 15 C contrast resolvable.
79func srh_inv_t_e7(t_c: i64) -> i64 {
80 var tk: i64 = 0
81 tk = stab_c_to_centi_kelvin(t_c)
82 if tk <= 0 { return SRH_INVALID }
83 return SRH_INV_T_SCALE / tk
84}
85
86// === Forward model ================================================
87
88// log10 of the rate ratio k(cond2)/k(cond1), in milli-decades.
89// Positive = condition 2 degrades faster.
90func srh_log10_rate_ratio_milli(e10_mk: i64, b10_milli: i64, t1_c: i64, rh1: i64, t2_c: i64, rh2: i64) -> i64 {
91 var u1: i64 = 0
92 var u2: i64 = 0
93 var thermal: i64 = 0
94 var moisture: i64 = 0
95 u1 = srh_inv_t_e7(t1_c)
96 u2 = srh_inv_t_e7(t2_c)
97 if u1 == SRH_INVALID { return SRH_INVALID }
98 if u2 == SRH_INVALID { return SRH_INVALID }
99 // E10 is milli-Kelvin and (u1-u2) is 1e-7 per Kelvin, so their product
100 // is 1e-10 dimensionless; 1e7 brings it to milli-decades.
101 thermal = e10_mk * (u1 - u2) / SRH_MAGIC_10000000
102 moisture = b10_milli * (rh2 - rh1)
103 return thermal + moisture
104}
105
106// Lifetime at condition 2 given a known lifetime at condition 1.
107// Lifetime is inversely proportional to rate, so it divides by the ratio.
108func srh_lifetime_months(ref_months: i64, e10_mk: i64, b10_milli: i64, t1_c: i64, rh1: i64, t2_c: i64, rh2: i64) -> i64 {
109 var lg: i64 = 0
110 var f: i64 = 0
111 if ref_months <= 0 { return SRH_INVALID }
112 lg = srh_log10_rate_ratio_milli(e10_mk, b10_milli, t1_c, rh1, t2_c, rh2)
113 if lg == SRH_INVALID { return SRH_INVALID }
114 f = ipow10_q3(0 - lg)
115 return ref_months * f / 1000
116}
117
118// === Separability =================================================
119
120// Determinant of the 2x2 system built from three conditions. Zero means the
121// three points are collinear in (1/T, RH) space and the two axes cannot be
122// told apart no matter how good the lifetime data is.
123func srh_determinant(t1_c: i64, rh1: i64, t2_c: i64, rh2: i64, t3_c: i64, rh3: i64) -> i64 {
124 var u1: i64 = 0
125 var u2: i64 = 0
126 var u3: i64 = 0
127 var a1: i64 = 0
128 var a2: i64 = 0
129 var b1: i64 = 0
130 var b2: i64 = 0
131 u1 = srh_inv_t_e7(t1_c)
132 u2 = srh_inv_t_e7(t2_c)
133 u3 = srh_inv_t_e7(t3_c)
134 if u1 == SRH_INVALID { return 0 }
135 if u2 == SRH_INVALID { return 0 }
136 if u3 == SRH_INVALID { return 0 }
137 a1 = u1 - u2
138 a2 = u1 - u3
139 b1 = (rh2 - rh1) * 1000
140 b2 = (rh3 - rh1) * 1000
141 return a1 * b2 - a2 * b1
142}
143
144// 1 = these three conditions can separate temperature from humidity.
145func srh_is_separable(t1_c: i64, rh1: i64, t2_c: i64, rh2: i64, t3_c: i64, rh3: i64) -> i64 {
146 if srh_determinant(t1_c, rh1, t2_c, rh2, t3_c, rh3) == 0 { return 0 }
147 return 1
148}
149
150// A design is well conditioned when the two contrasts point in genuinely
151// different directions in (1/T, RH) space. Returned in permil: 1000 means
152// the contrasts are as independent as they can be, near 0 means they are
153// nearly parallel and the separation, while arithmetically possible, will
154// amplify any error in the measured lifetimes enormously.
155//
156// ★NON-SINGULAR IS NOT THE SAME AS USABLE, and the gap between them is where
157// a plausible-looking activation energy comes from. Stepping temperature in
158// equal increments alongside humidity LOOKS collinear and is not, because
159// 1/T is nonlinear in T -- such a design passes the singularity test with a
160// determinant ~36x smaller than a good one. A caller that only asked "is it
161// separable" would get a yes and a meaningless answer.
162func srh_design_conditioning_permil(t1_c: i64, rh1: i64, t2_c: i64, rh2: i64, t3_c: i64, rh3: i64) -> i64 {
163 var u1: i64 = 0
164 var u2: i64 = 0
165 var u3: i64 = 0
166 var p1: i64 = 0
167 var p2: i64 = 0
168 var det: i64 = 0
169 var scale: i64 = 0
170 u1 = srh_inv_t_e7(t1_c)
171 u2 = srh_inv_t_e7(t2_c)
172 u3 = srh_inv_t_e7(t3_c)
173 if u1 == SRH_INVALID { return 0 }
174 if u2 == SRH_INVALID { return 0 }
175 if u3 == SRH_INVALID { return 0 }
176 p1 = (u1 - u2) * ((rh3 - rh1) * 1000)
177 p2 = (u1 - u3) * ((rh2 - rh1) * 1000)
178 det = p1 - p2
179 if det < 0 { det = 0 - det }
180 if p1 < 0 { p1 = 0 - p1 }
181 if p2 < 0 { p2 = 0 - p2 }
182 scale = p1 + p2
183 if scale == 0 { return 0 }
184 return det * 1000 / scale
185}
186
187// Minimum conditioning worth acting on. Below this the design is reported
188// as unusable even though the algebra would produce numbers.
189const SRH_MIN_CONDITIONING_PERMIL: i64 = 100
190
191func srh_design_is_well_conditioned(t1_c: i64, rh1: i64, t2_c: i64, rh2: i64, t3_c: i64, rh3: i64) -> i64 {
192 if srh_design_conditioning_permil(t1_c, rh1, t2_c, rh2, t3_c, rh3) < SRH_MIN_CONDITIONING_PERMIL { return 0 }
193 return 1
194}
195
196// ICH's storage design gives TWO conditions. Two conditions are one
197// contrast, and one contrast cannot resolve two axes. A function, not a
198// comment, so a planner cannot route around it and claim a pure Ea.
199func srh_ich_pair_alone_is_separable() -> i64 {
200 return 0
201}
202
203func srh_conditions_required() -> i64 {
204 return 3
205}
206
207// === The solver ===================================================
208
209// Recover BOTH axes from three (temperature, humidity, lifetime) points.
210// out[0] = E10 in milli-Kelvin, out[1] = B10 in milli-decades per %RH.
211// Returns 0 on success, SRH_UNDERDETERMINED if the design cannot separate.
212//
213// Scaling: E10 is solved in units of 10^4 milli-Kelvin so the widest
214// intermediate stays around 5e10, seven orders inside the i64 ceiling. The
215// cost is that E10 lands on a 10 Kelvin grid, which is 0.02% at the
216// magnitudes involved and is stated rather than hidden.
217func srh_separate(t1_c: i64, rh1: i64, m1: i64, t2_c: i64, rh2: i64, m2: i64, t3_c: i64, rh3: i64, m3: i64, out: *i64) -> i64 {
218 var u1: i64 = 0
219 var u2: i64 = 0
220 var u3: i64 = 0
221 var a1: i64 = 0
222 var a2: i64 = 0
223 var b1: i64 = 0
224 var b2: i64 = 0
225 var c1: i64 = 0
226 var c2: i64 = 0
227 var det: i64 = 0
228 var xs: i64 = 0
229 if m1 <= 0 { return SRH_INVALID }
230 if m2 <= 0 { return SRH_INVALID }
231 if m3 <= 0 { return SRH_INVALID }
232 u1 = srh_inv_t_e7(t1_c)
233 u2 = srh_inv_t_e7(t2_c)
234 u3 = srh_inv_t_e7(t3_c)
235 if u1 == SRH_INVALID { return SRH_INVALID }
236 if u2 == SRH_INVALID { return SRH_INVALID }
237 if u3 == SRH_INVALID { return SRH_INVALID }
238 // Observed log rate ratios, in milli-decades. Rate ratio = lifetime ratio
239 // inverted, so log10(k_i/k_1) = log10(m1/m_i).
240 c1 = ilog10_milli(m1 * 1000 / m2)
241 c2 = ilog10_milli(m1 * 1000 / m3)
242 if c1 == IP10_INVALID { return SRH_INVALID }
243 if c2 == IP10_INVALID { return SRH_INVALID }
244 a1 = u1 - u2
245 a2 = u1 - u3
246 b1 = (rh2 - rh1) * 1000
247 b2 = (rh3 - rh1) * 1000
248 det = a1 * b2 - a2 * b1
249 if det == 0 { return SRH_UNDERDETERMINED }
250 // Solve with E10 carried as X' = E10 / 10^4.
251 xs = (c1 * 1000 * b2 - c2 * 1000 * b1) / det
252 out[0] = xs * SRH_MAGIC_10000
253 out[1] = (a1 * c2 * 1000 - a2 * c1 * 1000) / det
254 return 0
255}
256
257// Activation energy in J/mol from a recovered E10. Inverse of
258// stab_e10_milli_kelvin, so the two files round-trip.
259func srh_ea_from_e10(e10_mk: i64) -> i64 {
260 var rl_milli: i64 = 0
261 if e10_mk <= 0 { return SRH_INVALID }
262 rl_milli = STAB_R_MILLI * STAB_LN10_MILLI / 1000
263 return e10_mk * rl_milli / SRH_MAGIC_1000000
264}
265
266// Does a recovered moisture sensitivity sit in the published range?
267func srh_b10_in_published_band(b10_milli: i64) -> i64 {
268 if b10_milli < SRH_PUB_B10_LO_MILLI { return 0 }
269 if b10_milli > SRH_PUB_B10_HI_MILLI { return 0 }
270 return 1
271}
272
273// What share of an observed acceleration (permil) is the moisture axis
274// rather than the temperature axis? This is the number that says how badly
275// a lumped parameter misattributes -- and the answer for ICH's own pair is
276// large enough that the lumping matters commercially.
277func srh_moisture_share_permil(e10_mk: i64, b10_milli: i64, t1_c: i64, rh1: i64, t2_c: i64, rh2: i64) -> i64 {
278 var total: i64 = 0
279 var moisture: i64 = 0
280 total = srh_log10_rate_ratio_milli(e10_mk, b10_milli, t1_c, rh1, t2_c, rh2)
281 if total == SRH_INVALID { return SRH_INVALID }
282 if total == 0 { return SRH_INVALID }
283 moisture = b10_milli * (rh2 - rh1)
284 return moisture * 1000 / total
285}