Mercurial > public > ostc4
annotate Common/Src/calc_crush.c @ 120:6347a86caa18 FlipDisplay
Cleanup warning and disable printf calls
author | Ideenmodellierer |
---|---|
date | Tue, 08 Jan 2019 22:31:17 +0100 |
parents | 8f8ea3a32e82 |
children |
rev | line source |
---|---|
38 | 1 /////////////////////////////////////////////////////////////////////////////// |
2 /// -*- coding: UTF-8 -*- | |
3 /// | |
4 /// \file Common/Src/calc_crush.c | |
5 /// \brief VPM Desaturation code | |
6 /// \author Heinrichs Weikamp | |
7 /// \date 2018 | |
8 /// | |
9 /// $Id$ | |
10 /////////////////////////////////////////////////////////////////////////////// | |
11 /// \par Copyright (c) 2014-2018 Heinrichs Weikamp gmbh | |
12 /// | |
13 /// This program is free software: you can redistribute it and/or modify | |
14 /// it under the terms of the GNU General Public License as published by | |
15 /// the Free Software Foundation, either version 3 of the License, or | |
16 /// (at your option) any later version. | |
17 /// | |
18 /// This program is distributed in the hope that it will be useful, | |
19 /// but WITHOUT ANY WARRANTY; without even the implied warranty of | |
20 /// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
21 /// GNU General Public License for more details. | |
22 /// | |
23 /// You should have received a copy of the GNU General Public License | |
24 /// along with this program. If not, see <http://www.gnu.org/licenses/>. | |
25 ////////////////////////////////////////////////////////////////////////////// | |
26 | |
27 #include "calc_crush.h" | |
28 | |
29 #include "decom.h" | |
30 #include "math.h" | |
31 #include "vpm.h" | |
32 | |
33 /* Common Block Declarations */ | |
34 //#pragma warning(disable:1035) | |
35 | |
36 const float SURFACE_TENSION_GAMMA = 0.0179f; //!Adj. Range: 0.015 to 0.065 N/m | |
37 const float SKIN_COMPRESSION_GAMMAC = 0.257f; //!Adj. Range: 0.160 to 0.290 N/m | |
38 const float UNITS_FACTOR = 10.1325f; | |
39 const float WATER_VAPOR_PRESSURE = 0.493f; // (Schreiner value) based on respiratory quotien | |
40 const float CRIT_VOLUME_PARAMETER_LAMBDA = 7500.0f; //!Adj. Range: 6500 to 8300 fsw-min | |
41 const float GRADIENT_ONSET_OF_IMPERM_ATM = 8.2f; //!Adj. Range: 5.0 to 10.0 atm | |
42 const float REGENERATION_TIME_CONSTANT = 20160.0f; //!Adj. Range: 10080 to 51840 min | |
43 const float PRESSURE_OTHER_GASES_MMHG = 102.0f; //!Constant value for PO2 up to 2 atm | |
44 const float CONSTANT_PRESSURE_OTHER_GASES = 102.0f * 10.1325f / 760.0f; // PRESSURE_OTHER_GASES_MMHG / 760. * UNITS_FACTOR; | |
45 | |
46 const float HELIUM_TIME_CONSTANT[16] = {3.68695308808482E-001f, | |
47 2.29518933960247E-001f, | |
48 1.46853216220327E-001f, | |
49 9.91626867753856E-002f, | |
50 6.78890480470074E-002f, | |
51 4.78692804254106E-002f, | |
52 3.37626488338989E-002f, | |
53 2.38113081607676E-002f, | |
54 1.68239606932026E-002f, | |
55 1.25592893741610E-002f, | |
56 9.80544886914621E-003f, | |
57 7.67264977374303E-003f, | |
58 6.01220557342307E-003f, | |
59 4.70185307665137E-003f, | |
60 3.68225234041620E-003f, | |
61 2.88775228329769E-003f}; | |
62 | |
63 const float NITROGEN_TIME_CONSTANT[16] = {1.38629436111989E-001f, | |
64 8.66433975699932E-002f, | |
65 5.54517744447956E-002f, | |
66 3.74674151654024E-002f, | |
67 2.56721177985165E-002f, | |
68 1.80978376125312E-002f, | |
69 1.27651414467762E-002f, | |
70 9.00191143584345E-003f, | |
71 6.35914844550409E-003f, | |
72 4.74758342849278E-003f, | |
73 3.70666941475907E-003f, | |
74 2.90019740820061E-003f, | |
75 2.27261370675392E-003f, | |
76 1.77730046297422E-003f, | |
77 1.39186180835330E-003f, | |
78 1.09157036308653E-003f}; | |
79 | |
80 int onset_of_impermeability(SGas* pGas, float *starting_ambient_pressure, | |
81 float *ending_ambient_pressure, | |
82 float *rate, | |
83 float* amb_pressure_onset_of_imperm, | |
84 float* gas_tension_onset_of_imperm, | |
85 float* initial_helium_pressure, | |
86 float* initial_nitrogen_pressure, | |
87 short i); | |
88 int radius_root_finder (float *a, float *b, float *c, float *low_bound, float *high_bound, float *ending_radius); | |
89 //void get_inert_gases_(SBuehlmann* input, ,short gas_id, float ambient_pressure_bar, float* fraction_nitrogen,float* fraction_helium ); | |
90 int vpm_repetitive_algorithm(SVpm* pVpm, float *surface_interval_time, float* initial_critical_radius_he, float* initial_critical_radius_n2); | |
91 | |
92 | |
93 | |
94 /* =============================================================================== */ | |
95 /* NOTE ABOUT PRESSURE UNITS USED IN CALCULATIONS: */ | |
96 /* It is the convention in decompression calculations to compute all gas */ | |
97 /* loadings, absolute pressures, partial pressures, etc., in the units of */ | |
98 /* depth pressure that you are diving - either feet of seawater (fsw) or */ | |
99 /* meters of seawater (msw). This program follows that convention with the */ | |
100 /* the exception that all VPM calculations are performed in SI units (by */ | |
101 /* necessity). Accordingly, there are several conversions back and forth */ | |
102 /* between the diving pressure units and the SI units. */ | |
103 /* =============================================================================== */ | |
104 /* =============================================================================== */ | |
105 /* FUNCTION SUBPROGRAM FOR GAS LOADING CALCULATIONS - ASCENT AND DESCENT */ | |
106 /* =============================================================================== */ | |
107 | |
108 float schreiner_equation__2(float *initial_inspired_gas_pressure, | |
109 float *rate_change_insp_gas_pressure, | |
110 float *interval_time_minutes, | |
111 const float *gas_time_constant, | |
112 float *initial_gas_pressure) | |
113 { | |
114 /* System generated locals */ | |
115 float ret_val; | |
116 float time_null_pressure = 0.0f; | |
117 float time_rest = 0.0f; | |
118 float time = *interval_time_minutes; | |
119 /* =============================================================================== */ | |
120 /* Note: The Schreiner equation is applied when calculating the uptake or */ | |
121 /* elimination of compartment gases during linear ascents or descents at a */ | |
122 /* constant rate. For ascents, a negative number for rate must be used. */ | |
123 /* =============================================================================== */ | |
124 if( *rate_change_insp_gas_pressure < 0.0f) | |
125 { | |
126 time_null_pressure = -1.0f * *initial_inspired_gas_pressure / *rate_change_insp_gas_pressure; | |
127 if(time > time_null_pressure ) | |
128 { | |
129 time_rest = time - time_null_pressure; | |
130 time = time_null_pressure; | |
131 } | |
132 } | |
133 ret_val = | |
134 *initial_inspired_gas_pressure + | |
135 *rate_change_insp_gas_pressure * | |
136 (time - 1.f / *gas_time_constant) - | |
137 (*initial_inspired_gas_pressure - | |
138 *initial_gas_pressure - | |
139 *rate_change_insp_gas_pressure / *gas_time_constant) * | |
140 expf(-(*gas_time_constant) * time); | |
141 | |
142 if(time_rest > 0.0f) | |
143 { | |
144 ret_val = ret_val * expf(-(*gas_time_constant) * time_rest); | |
145 } | |
146 | |
147 | |
148 return ret_val; | |
149 }; /* schreiner_equation__2 */ | |
150 | |
151 /* =============================================================================== */ | |
152 /* SUBROUTINE CALC_CRUSHING_PRESSURE */ | |
153 /* Purpose: Compute the effective "crushing pressure" in each compartment as */ | |
154 /* a result of descent segment(s). The crushing pressure is the gradient */ | |
155 /* (difference in pressure) between the outside ambient pressure and the */ | |
156 /* gas tension inside a VPM nucleus (bubble seed). This gradient acts to */ | |
157 /* reduce (shrink) the radius smaller than its initial value at the surface. */ | |
158 /* This phenomenon has important ramifications because the smaller the radius */ | |
159 /* of a VPM nucleus, the greater the allowable supersaturation gradient upon */ | |
160 /* ascent. Gas loading (uptake) during descent, especially in the fast */ | |
161 /* compartments, will reduce the magnitude of the crushing pressure. The */ | |
162 /* crushing pressure is not cumulative over a multi-level descent. It will */ | |
163 /* be the maximum value obtained in any one discrete segment of the overall */ | |
164 /* descent. Thus, the program must compute and store the maximum crushing */ | |
165 /* pressure for each compartment that was obtained across all segments of */ | |
166 /* the descent profile. */ | |
167 | |
168 /* The calculation of crushing pressure will be different depending on */ | |
169 /* whether or not the gradient is in the VPM permeable range (gas can diffuse */ | |
170 /* across skin of VPM nucleus) or the VPM impermeable range (molecules in */ | |
171 /* skin of nucleus are squeezed together so tight that gas can no longer */ | |
172 /* diffuse in or out of nucleus; the gas becomes trapped and further resists */ | |
173 /* the crushing pressure). The solution for crushing pressure in the VPM */ | |
174 /* permeable range is a simple linear equation. In the VPM impermeable */ | |
175 /* range, a cubic equation must be solved using a numerical method. */ | |
176 | |
177 /* Separate crushing pressures are tracked for helium and nitrogen because */ | |
178 /* they can have different critical radii. The crushing pressures will be */ | |
179 /* the same for helium and nitrogen in the permeable range of the model, but */ | |
180 /* they will start to diverge in the impermeable range. This is due to */ | |
181 /* the differences between starting radius, radius at the onset of */ | |
182 /* impermeability, and radial compression in the impermeable range. */ | |
183 /* =============================================================================== */ | |
184 int calc_crushing_pressure(SLifeData* lifeData, SVpm* vpm, float * initial_helium_pressure, float * initial_nitrogen_pressure, float starting_ambient_pressure, | |
185 float rate ) | |
186 { | |
187 /* System generated locals */ | |
188 static float r1, r2; | |
189 | |
190 static float low_bound_n2, | |
191 ending_radius_n2, | |
192 gradient_onset_of_imperm_pa; | |
193 static float low_bound_he, | |
194 ending_radius_he, | |
195 high_bound_n2, | |
196 crushing_pressure_n2; | |
197 short i; | |
198 static float crushing_pressure_pascals_n2, | |
199 gradient_onset_of_imperm, | |
200 starting_gas_tension, | |
201 high_bound_he, | |
202 crushing_pressure_he, | |
203 amb_press_onset_of_imperm_pa, | |
204 crushing_pressure_pascals_he, | |
205 radius_onset_of_imperm_n2, | |
206 starting_gradient, | |
207 radius_onset_of_imperm_he, | |
208 ending_gas_tension; | |
209 | |
210 static float ending_ambient_pressure_pa, | |
211 a_n2, | |
212 b_n2, | |
213 c_n2, | |
214 ending_gradient, | |
215 gas_tension_onset_of_imperm_pa, | |
216 a_he, | |
217 b_he, c_he; | |
218 static float amb_pressure_onset_of_imperm[16]; | |
219 static float gas_tension_onset_of_imperm[16]; | |
220 | |
221 static float helium_pressure_crush[16]; | |
222 static float nitrogen_pressure_crush[16]; | |
223 | |
224 | |
225 static float ending_ambient_pressure = 0; | |
226 | |
227 ending_ambient_pressure = lifeData->pressure_ambient_bar * 10; | |
228 for( i = 0; i < 16; i++) | |
229 { | |
230 helium_pressure_crush[i] = lifeData->tissue_helium_bar[i] * 10; | |
231 nitrogen_pressure_crush[i] = lifeData->tissue_nitrogen_bar[i] * 10; | |
232 } | |
233 | |
234 | |
235 | |
236 | |
237 | |
238 /* loop */ | |
239 /* =============================================================================== */ | |
240 /* CALCULATIONS */ | |
241 /* First, convert the Gradient for Onset of Impermeability from units of */ | |
242 /* atmospheres to diving pressure units (either fsw or msw) and to Pascals */ | |
243 /* (SI units). The reason that the Gradient for Onset of Impermeability is */ | |
244 /* given in the program settings in units of atmospheres is because that is */ | |
245 /* how it was reported in the original research papers by Yount and */ | |
246 /* colleauges. */ | |
247 /* =============================================================================== */ | |
248 | |
249 gradient_onset_of_imperm = GRADIENT_ONSET_OF_IMPERM_ATM * UNITS_FACTOR; | |
250 gradient_onset_of_imperm_pa = GRADIENT_ONSET_OF_IMPERM_ATM * 101325.0f; | |
251 | |
252 /* =============================================================================== */ | |
253 /* Assign values of starting and ending ambient pressures for descent segment */ | |
254 /* =============================================================================== */ | |
255 | |
256 //starting_ambient_pressure = *starting_depth; | |
257 //ending_ambient_pressure = *ending_depth; | |
258 | |
259 /* =============================================================================== */ | |
260 /* MAIN LOOP WITH NESTED DECISION TREE */ | |
261 /* For each compartment, the program computes the starting and ending */ | |
262 /* gas tensions and gradients. The VPM is different than some dissolved gas */ | |
263 /* algorithms, Buhlmann for example, in that it considers the pressure due to */ | |
264 /* oxygen, carbon dioxide, and water vapor in each compartment in addition to */ | |
265 /* the inert gases helium and nitrogen. These "other gases" are included in */ | |
266 /* the calculation of gas tensions and gradients. */ | |
267 /* =============================================================================== */ | |
268 | |
269 crushing_pressure_he = 0.0f; | |
270 crushing_pressure_n2 = 0.0f; | |
271 | |
272 for (i = 0; i < 16; ++i) { | |
273 starting_gas_tension = initial_helium_pressure[i] + initial_nitrogen_pressure[i] + CONSTANT_PRESSURE_OTHER_GASES; | |
274 starting_gradient = starting_ambient_pressure - starting_gas_tension; | |
275 ending_gas_tension = helium_pressure_crush[i] + nitrogen_pressure_crush[i] + CONSTANT_PRESSURE_OTHER_GASES; | |
276 ending_gradient = ending_ambient_pressure - ending_gas_tension; | |
277 | |
278 /* =============================================================================== */ | |
279 /* Compute radius at onset of impermeability for helium and nitrogen */ | |
280 /* critical radii */ | |
281 /* =============================================================================== */ | |
282 | |
283 radius_onset_of_imperm_he = 1.0f / ( | |
284 gradient_onset_of_imperm_pa / | |
285 ((SKIN_COMPRESSION_GAMMAC - | |
286 SURFACE_TENSION_GAMMA) * 2.0f) + | |
287 1.0f / vpm->adjusted_critical_radius_he[i]); | |
288 radius_onset_of_imperm_n2 = 1.0f / ( | |
289 gradient_onset_of_imperm_pa / | |
290 ((SKIN_COMPRESSION_GAMMAC - | |
291 SURFACE_TENSION_GAMMA) * 2.0f) + | |
292 1.0f / vpm->adjusted_critical_radius_n2[i]); | |
293 | |
294 /* =============================================================================== */ | |
295 /* FIRST BRANCH OF DECISION TREE - PERMEABLE RANGE */ | |
296 /* Crushing pressures will be the same for helium and nitrogen */ | |
297 /* =============================================================================== */ | |
298 | |
299 if (ending_gradient <= gradient_onset_of_imperm) { | |
300 crushing_pressure_he = ending_ambient_pressure - ending_gas_tension; | |
301 crushing_pressure_n2 = ending_ambient_pressure - ending_gas_tension; | |
302 } | |
303 | |
304 /* =============================================================================== */ | |
305 /* SECOND BRANCH OF DECISION TREE - IMPERMEABLE RANGE */ | |
306 /* Both the ambient pressure and the gas tension at the onset of */ | |
307 /* impermeability must be computed in order to properly solve for the ending */ | |
308 /* radius and resultant crushing pressure. The first decision block */ | |
309 /* addresses the special case when the starting gradient just happens to be */ | |
310 /* equal to the gradient for onset of impermeability (not very likely!). */ | |
311 /* =============================================================================== */ | |
312 | |
313 if (ending_gradient > gradient_onset_of_imperm) { | |
314 if (starting_gradient == gradient_onset_of_imperm) { | |
315 amb_pressure_onset_of_imperm[i] = starting_ambient_pressure; | |
316 gas_tension_onset_of_imperm[i] = starting_gas_tension; | |
317 } | |
318 | |
319 /* =============================================================================== */ | |
320 /* In most cases, a subroutine will be called to find these values using a */ | |
321 /* numerical method. */ | |
322 /* =============================================================================== */ | |
323 | |
324 if (starting_gradient < gradient_onset_of_imperm) { | |
325 onset_of_impermeability(&(lifeData->actualGas), &starting_ambient_pressure, &ending_ambient_pressure, &rate, | |
326 amb_pressure_onset_of_imperm, gas_tension_onset_of_imperm, | |
327 initial_helium_pressure, initial_nitrogen_pressure, i); | |
328 } | |
329 | |
330 /* =============================================================================== */ | |
331 /* Next, using the values for ambient pressure and gas tension at the onset */ | |
332 /* of impermeability, the equations are set up to process the calculations */ | |
333 /* through the radius root finder subroutine. This subprogram will find the */ | |
334 /* root (solution) to the cubic equation using a numerical method. In order */ | |
335 /* to do this efficiently, the equations are placed in the form */ | |
336 /* Ar^3 - Br^2 - C = 0, where r is the ending radius after impermeable */ | |
337 /* compression. The coefficients A, B, and C for helium and nitrogen are */ | |
338 /* computed and passed to the subroutine as arguments. The high and low */ | |
339 /* bounds to be used by the numerical method of the subroutine are also */ | |
340 /* computed (see separate page posted on Deco List ftp site entitled */ | |
341 /* "VPM: Solving for radius in the impermeable regime"). The subprogram */ | |
342 /* will return the value of the ending radius and then the crushing */ | |
343 /* pressures for helium and nitrogen can be calculated. */ | |
344 /* =============================================================================== */ | |
345 | |
346 ending_ambient_pressure_pa = ending_ambient_pressure / UNITS_FACTOR * 101325.0f; | |
347 amb_press_onset_of_imperm_pa = amb_pressure_onset_of_imperm[i] / UNITS_FACTOR * 101325.0f; | |
348 gas_tension_onset_of_imperm_pa = gas_tension_onset_of_imperm[i] / UNITS_FACTOR * 101325.0f; | |
349 b_he = (SKIN_COMPRESSION_GAMMAC - SURFACE_TENSION_GAMMA) * 2.0f; | |
350 a_he = ending_ambient_pressure_pa - amb_press_onset_of_imperm_pa + gas_tension_onset_of_imperm_pa | |
351 + (SKIN_COMPRESSION_GAMMAC - SURFACE_TENSION_GAMMA) * 2.0f / radius_onset_of_imperm_he; | |
352 /* Computing 3rd power */ | |
353 r1 = radius_onset_of_imperm_he; | |
354 c_he = gas_tension_onset_of_imperm_pa * (r1 * (r1 * r1)); | |
355 high_bound_he = radius_onset_of_imperm_he; | |
356 low_bound_he = b_he / a_he; | |
357 radius_root_finder(&a_he, &b_he, &c_he, &low_bound_he, &high_bound_he, &ending_radius_he); | |
358 /* Computing 3rd power */ | |
359 r1 = radius_onset_of_imperm_he; | |
360 /* Computing 3rd power */ | |
361 r2 = ending_radius_he; | |
362 crushing_pressure_pascals_he = | |
363 gradient_onset_of_imperm_pa + | |
364 ending_ambient_pressure_pa - | |
365 amb_press_onset_of_imperm_pa + | |
366 gas_tension_onset_of_imperm_pa * | |
367 (1.0f - r1 * (r1 * r1) / (r2 * (r2 * r2))); | |
368 crushing_pressure_he = | |
369 crushing_pressure_pascals_he / 101325.0f * UNITS_FACTOR; | |
370 b_n2 = (SKIN_COMPRESSION_GAMMAC - SURFACE_TENSION_GAMMA) * 2.0f; | |
371 a_n2 = ending_ambient_pressure_pa - | |
372 amb_press_onset_of_imperm_pa + | |
373 gas_tension_onset_of_imperm_pa + | |
374 (SKIN_COMPRESSION_GAMMAC - SURFACE_TENSION_GAMMA) * | |
375 2.0f / radius_onset_of_imperm_n2; | |
376 /* Computing 3rd power */ | |
377 r1 = radius_onset_of_imperm_n2; | |
378 c_n2 = gas_tension_onset_of_imperm_pa * (r1 * (r1 * r1)); | |
379 high_bound_n2 = radius_onset_of_imperm_n2; | |
380 low_bound_n2 = b_n2 / a_n2; | |
381 radius_root_finder(&a_n2, | |
382 &b_n2, | |
383 &c_n2, | |
384 &low_bound_n2, | |
385 &high_bound_n2, | |
386 &ending_radius_n2); | |
387 | |
388 /* Computing 3rd power */ | |
389 r1 = radius_onset_of_imperm_n2; | |
390 /* Computing 3rd power */ | |
391 r2 = ending_radius_n2; | |
392 crushing_pressure_pascals_n2 = | |
393 gradient_onset_of_imperm_pa + | |
394 ending_ambient_pressure_pa - | |
395 amb_press_onset_of_imperm_pa + | |
396 gas_tension_onset_of_imperm_pa * (1.0f - r1 * | |
397 (r1 * r1) / (r2 * (r2 * r2))); | |
398 crushing_pressure_n2 = crushing_pressure_pascals_n2 / 101325.0f * UNITS_FACTOR; | |
399 } | |
400 | |
401 /* =============================================================================== */ | |
402 /* UPDATE VALUES OF MAX CRUSHING PRESSURE IN GLOBAL ARRAYS */ | |
403 /* =============================================================================== */ | |
404 | |
405 /* Computing MAX */ | |
406 r1 = vpm->max_crushing_pressure_he[i]; | |
407 vpm->max_crushing_pressure_he[i] = fmaxf(r1, crushing_pressure_he); | |
408 /* Computing MAX */ | |
409 r1 = vpm->max_crushing_pressure_n2[i]; | |
410 vpm->max_crushing_pressure_n2[i] = fmaxf(r1, crushing_pressure_n2); | |
411 } | |
412 return 0; | |
413 } /* calc_crushing_pressure */ | |
414 | |
415 /* =============================================================================== */ | |
416 /* SUBROUTINE ONSET_OF_IMPERMEABILITY */ | |
417 /* Purpose: This subroutine uses the Bisection Method to find the ambient */ | |
418 /* pressure and gas tension at the onset of impermeability for a given */ | |
419 /* compartment. Source: "Numerical Recipes in Fortran 77", */ | |
420 /* Cambridge University Press, 1992. */ | |
421 /* =============================================================================== */ | |
422 | |
423 int onset_of_impermeability(SGas* pGas, float *starting_ambient_pressure, | |
424 float *ending_ambient_pressure, | |
425 float *rate, | |
426 float* amb_pressure_onset_of_imperm, | |
427 float* gas_tension_onset_of_imperm, | |
428 float* initial_helium_pressure, | |
429 float* initial_nitrogen_pressure, | |
430 short i) | |
431 { | |
432 /* Local variables */ | |
433 float time, last_diff_change, mid_range_nitrogen_pressure; | |
434 short j; | |
435 float gas_tension_at_mid_range, | |
436 initial_inspired_n2_pressure, | |
437 gradient_onset_of_imperm, | |
438 starting_gas_tension, | |
439 low_bound, | |
440 initial_inspired_he_pressure, | |
441 high_bound_nitrogen_pressure, | |
442 nitrogen_rate, | |
443 function_at_mid_range, | |
444 function_at_low_bound, | |
445 high_bound, | |
446 mid_range_helium_pressure, | |
447 mid_range_time, | |
448 ending_gas_tension, | |
449 function_at_high_bound; | |
450 | |
451 float mid_range_ambient_pressure, | |
452 high_bound_helium_pressure, | |
453 helium_rate, | |
454 differential_change; | |
455 float fraction_helium_begin; | |
456 float fraction_helium_end; | |
457 float fraction_nitrogen_begin; | |
458 float fraction_nitrogen_end; | |
459 /* loop */ | |
460 /* =============================================================================== */ | |
461 /* CALCULATIONS */ | |
462 /* First convert the Gradient for Onset of Impermeability to the diving */ | |
463 /* pressure units that are being used */ | |
464 /* =============================================================================== */ | |
465 | |
466 gradient_onset_of_imperm = GRADIENT_ONSET_OF_IMPERM_ATM * UNITS_FACTOR; | |
467 | |
468 /* =============================================================================== */ | |
469 /* ESTABLISH THE BOUNDS FOR THE ROOT SEARCH USING THE BISECTION METHOD */ | |
470 /* In this case, we are solving for time - the time when the ambient pressure */ | |
471 /* minus the gas tension will be equal to the Gradient for Onset of */ | |
472 /* Impermeabliity. The low bound for time is set at zero and the high */ | |
473 /* bound is set at the elapsed time (segment time) it took to go from the */ | |
474 /* starting ambient pressure to the ending ambient pressure. The desired */ | |
475 /* ambient pressure and gas tension at the onset of impermeability will */ | |
476 /* be found somewhere between these endpoints. The algorithm checks to */ | |
477 /* make sure that the solution lies in between these bounds by first */ | |
478 /* computing the low bound and high bound function values. */ | |
479 /* =============================================================================== */ | |
480 | |
481 /*initial_inspired_he_pressure = | |
482 (*starting_ambient_pressure - water_vapor_pressure) * fraction_helium[mix_number - 1]; | |
483 initial_inspired_n2_pressure = | |
484 (*starting_ambient_pressure - water_vapor_pressure) * fraction_nitrogen[mix_number - 1]; | |
485 helium_rate = *rate * fraction_helium[mix_number - 1]; | |
486 nitrogen_rate = *rate * fraction_nitrogen[mix_number - 1];*/ | |
487 low_bound = 0.; | |
488 high_bound = (*ending_ambient_pressure - *starting_ambient_pressure) / *rate; | |
489 | |
490 //New | |
491 decom_get_inert_gases( *starting_ambient_pressure / 10.0f, pGas, &fraction_nitrogen_begin, &fraction_helium_begin ); | |
492 decom_get_inert_gases(*ending_ambient_pressure / 10.0f, pGas, &fraction_nitrogen_end, &fraction_helium_end ); | |
493 initial_inspired_he_pressure = (*starting_ambient_pressure - WATER_VAPOR_PRESSURE) * fraction_helium_begin; | |
494 initial_inspired_n2_pressure = (*starting_ambient_pressure - WATER_VAPOR_PRESSURE) * fraction_nitrogen_begin; | |
495 helium_rate = ((*ending_ambient_pressure - WATER_VAPOR_PRESSURE)* fraction_helium_end - initial_inspired_he_pressure)/high_bound; | |
496 nitrogen_rate = ((*ending_ambient_pressure - WATER_VAPOR_PRESSURE)* fraction_nitrogen_end - initial_inspired_n2_pressure)/high_bound; | |
497 | |
498 starting_gas_tension = | |
499 initial_helium_pressure[i] + | |
500 initial_nitrogen_pressure[i] + | |
501 CONSTANT_PRESSURE_OTHER_GASES; | |
502 function_at_low_bound = | |
503 *starting_ambient_pressure - | |
504 starting_gas_tension - | |
505 gradient_onset_of_imperm; | |
506 high_bound_helium_pressure = | |
507 schreiner_equation__2(&initial_inspired_he_pressure, | |
508 &helium_rate, | |
509 &high_bound, | |
510 &HELIUM_TIME_CONSTANT[i], | |
511 &initial_helium_pressure[i]); | |
512 high_bound_nitrogen_pressure = | |
513 schreiner_equation__2(&initial_inspired_n2_pressure, | |
514 &nitrogen_rate, | |
515 &high_bound, | |
516 &NITROGEN_TIME_CONSTANT[i], | |
517 &initial_nitrogen_pressure[i]); | |
518 ending_gas_tension = | |
519 high_bound_helium_pressure + | |
520 high_bound_nitrogen_pressure + | |
521 CONSTANT_PRESSURE_OTHER_GASES; | |
522 function_at_high_bound = | |
523 *ending_ambient_pressure - | |
524 ending_gas_tension - | |
525 gradient_onset_of_imperm; | |
526 if (function_at_high_bound * function_at_low_bound >= 0.0f) { | |
527 //printf("\nERROR! ROOT IS NOT WITHIN BRACKETS"); | |
528 } | |
529 | |
530 /* =============================================================================== */ | |
531 /* APPLY THE BISECTION METHOD IN SEVERAL ITERATIONS UNTIL A SOLUTION WITH */ | |
532 /* THE DESIRED ACCURACY IS FOUND */ | |
533 /* Note: the program allows for up to 100 iterations. Normally an exit will */ | |
534 /* be made from the loop well before that number. If, for some reason, the */ | |
535 /* program exceeds 100 iterations, there will be a pause to alert the user. */ | |
536 /* =============================================================================== */ | |
537 | |
538 if (function_at_low_bound < 0.0f) { | |
539 time = low_bound; | |
540 differential_change = high_bound - low_bound; | |
541 } else { | |
542 time = high_bound; | |
543 differential_change = low_bound - high_bound; | |
544 } | |
545 for (j = 1; j <= 100; ++j) { | |
546 last_diff_change = differential_change; | |
547 differential_change = last_diff_change * 0.5f; | |
548 mid_range_time = time + differential_change; | |
549 mid_range_ambient_pressure = *starting_ambient_pressure + *rate * mid_range_time; | |
550 mid_range_helium_pressure = | |
551 schreiner_equation__2(&initial_inspired_he_pressure, | |
552 &helium_rate, | |
553 &mid_range_time, | |
554 &HELIUM_TIME_CONSTANT[i], | |
555 &initial_helium_pressure[i]); | |
556 mid_range_nitrogen_pressure = | |
557 schreiner_equation__2(&initial_inspired_n2_pressure, | |
558 &nitrogen_rate, | |
559 &mid_range_time, | |
560 &NITROGEN_TIME_CONSTANT[i], | |
561 &initial_nitrogen_pressure[i]); | |
562 gas_tension_at_mid_range = | |
563 mid_range_helium_pressure + | |
564 mid_range_nitrogen_pressure + | |
565 CONSTANT_PRESSURE_OTHER_GASES; | |
566 function_at_mid_range = | |
567 mid_range_ambient_pressure - | |
568 gas_tension_at_mid_range - | |
569 gradient_onset_of_imperm; | |
570 if (function_at_mid_range <= 0.0f) { | |
571 time = mid_range_time; | |
572 } | |
573 if (fabs(differential_change) < .001f || | |
574 function_at_mid_range == 0.0f) { | |
575 goto L100; | |
576 } | |
577 } | |
578 //printf("\nERROR! ROOT SEARCH EXCEEDED MAXIMUM ITERATIONS"); | |
579 | |
580 /* =============================================================================== */ | |
581 /* When a solution with the desired accuracy is found, the program jumps out */ | |
582 /* of the loop to Line 100 and assigns the solution values for ambient */ | |
583 /* pressure and gas tension at the onset of impermeability. */ | |
584 /* =============================================================================== */ | |
585 | |
586 L100: | |
587 amb_pressure_onset_of_imperm[i] = mid_range_ambient_pressure; | |
588 gas_tension_onset_of_imperm[i] = gas_tension_at_mid_range; | |
589 return 0; | |
590 } /* onset_of_impermeability */ | |
591 | |
592 | |
593 /* =============================================================================== */ | |
594 /* SUBROUTINE RADIUS_ROOT_FINDER */ | |
595 /* Purpose: This subroutine is a "fail-safe" routine that combines the */ | |
596 /* Bisection Method and the Newton-Raphson Method to find the desired root. */ | |
597 /* This hybrid algorithm takes a bisection step whenever Newton-Raphson would */ | |
598 /* take the solution out of bounds, or whenever Newton-Raphson is not */ | |
599 /* converging fast enough. Source: "Numerical Recipes in Fortran 77", */ | |
600 /* Cambridge University Press, 1992. */ | |
601 /* =============================================================================== */ | |
602 | |
603 int radius_root_finder (float *a, | |
604 float *b, | |
605 float *c, | |
606 float *low_bound, | |
607 float *high_bound, | |
608 float *ending_radius) | |
609 { | |
610 /* System generated locals */ | |
611 float r1, r2; | |
612 | |
613 /* Local variables */ | |
614 float radius_at_low_bound, | |
615 last_diff_change, | |
616 function, | |
617 radius_at_high_bound; | |
618 short i; | |
619 float function_at_low_bound, | |
620 last_ending_radius, | |
621 function_at_high_bound, | |
622 derivative_of_function, | |
623 differential_change; | |
624 | |
625 /* loop */ | |
626 /* =============================================================================== */ | |
627 /* BEGIN CALCULATIONS BY MAKING SURE THAT THE ROOT LIES WITHIN BOUNDS */ | |
628 /* In this case we are solving for radius in a cubic equation of the form, */ | |
629 /* Ar^3 - Br^2 - C = 0. The coefficients A, B, and C were passed to this */ | |
630 /* subroutine as arguments. */ | |
631 /* =============================================================================== */ | |
632 | |
633 function_at_low_bound = | |
634 *low_bound * (*low_bound * (*a * *low_bound - *b)) - *c; | |
635 function_at_high_bound = | |
636 *high_bound * (*high_bound * (*a * *high_bound - *b)) - *c; | |
637 if (function_at_low_bound > 0.0f && function_at_high_bound > 0.0f) { | |
638 // printf("\nERROR! ROOT IS NOT WITHIN BRACKETS"); | |
639 | |
640 } | |
641 | |
642 /* =============================================================================== */ | |
643 /* Next the algorithm checks for special conditions and then prepares for */ | |
644 /* the first bisection. */ | |
645 /* =============================================================================== */ | |
646 | |
647 if (function_at_low_bound < 0.0f && function_at_high_bound < 0.0f) { | |
648 //printf("\nERROR! ROOT IS NOT WITHIN BRACKETS"); | |
649 | |
650 } | |
651 if (function_at_low_bound == 0.0f) { | |
652 *ending_radius = *low_bound; | |
653 return 0; | |
654 } else if (function_at_high_bound == 0.0f) { | |
655 *ending_radius = *high_bound; | |
656 return 0; | |
657 } else if (function_at_low_bound < 0.0f) { | |
658 radius_at_low_bound = *low_bound; | |
659 radius_at_high_bound = *high_bound; | |
660 } else { | |
661 radius_at_high_bound = *low_bound; | |
662 radius_at_low_bound = *high_bound; | |
663 } | |
664 *ending_radius = (*low_bound + *high_bound) * .5f; | |
665 last_diff_change = (r1 = *high_bound - *low_bound, fabs(r1)); | |
666 differential_change = last_diff_change; | |
667 | |
668 /* =============================================================================== */ | |
669 /* At this point, the Newton-Raphson Method is applied which uses a function */ | |
670 /* and its first derivative to rapidly converge upon a solution. */ | |
671 /* Note: the program allows for up to 100 iterations. Normally an exit will */ | |
672 /* be made from the loop well before that number. If, for some reason, the */ | |
673 /* program exceeds 100 iterations, there will be a pause to alert the user. */ | |
674 /* When a solution with the desired accuracy is found, exit is made from the */ | |
675 /* loop by returning to the calling program. The last value of ending */ | |
676 /* radius has been assigned as the solution. */ | |
677 /* =============================================================================== */ | |
678 | |
679 function = | |
680 *ending_radius * (*ending_radius * (*a * *ending_radius - *b)) - *c; | |
681 derivative_of_function = | |
682 *ending_radius * (*ending_radius * 3.0f * *a - *b * 2.0f); | |
683 for (i = 1; i <= 100; ++i) { | |
684 if (((*ending_radius - radius_at_high_bound) * derivative_of_function - function) * | |
685 ((*ending_radius - radius_at_low_bound) * derivative_of_function - function) >= 0.0f | |
686 || (r1 = function * 2.0f, fabs(r1)) > | |
687 (r2 = last_diff_change * derivative_of_function, fabs(r2))) { | |
688 last_diff_change = differential_change; | |
689 differential_change = | |
690 (radius_at_high_bound - radius_at_low_bound) * .5f; | |
691 *ending_radius = radius_at_low_bound + differential_change; | |
692 if (radius_at_low_bound == *ending_radius) { | |
693 return 0; | |
694 } | |
695 } else { | |
696 last_diff_change = differential_change; | |
697 differential_change = function / derivative_of_function; | |
698 last_ending_radius = *ending_radius; | |
699 *ending_radius -= differential_change; | |
700 if (last_ending_radius == *ending_radius) { | |
701 return 0; | |
702 } | |
703 } | |
704 if (fabs(differential_change) < 1e-12) { | |
705 return 0; | |
706 } | |
707 function = | |
708 *ending_radius * (*ending_radius * (*a * *ending_radius - *b)) - *c; | |
709 derivative_of_function = | |
710 *ending_radius * (*ending_radius * 3.0f * *a - *b * 2.0f); | |
711 if (function < 0.0f) { | |
712 radius_at_low_bound = *ending_radius; | |
713 } else { | |
714 radius_at_high_bound = *ending_radius; | |
715 } | |
716 } | |
717 // printf("\nERROR! ROOT SEARCH EXCEEDED MAXIMUM ITERATIONS"); | |
718 return 0; | |
719 } /* radius_root_finder */ | |
720 | |
721 | |
722 | |
723 | |
724 void vpm_init(SVpm* pVpm, short conservatism, short repetitive_dive, long seconds_since_last_dive) | |
725 { | |
726 | |
51
8f8ea3a32e82
Resolved warnings pointing to possible invalid memory access
Ideenmodellierer
parents:
38
diff
changeset
|
727 float critical_radius_n2_microns = 0.82; /* be conservative in case of an unexpected parameter value */ |
8f8ea3a32e82
Resolved warnings pointing to possible invalid memory access
Ideenmodellierer
parents:
38
diff
changeset
|
728 float critical_radius_he_microns = 0.72; |
38 | 729 float initial_critical_radius_n2[16]; |
730 float initial_critical_radius_he[16]; | |
731 int i = 0; | |
732 float surface_time = seconds_since_last_dive / 60; | |
733 pVpm->repetitive_variables_not_valid = !repetitive_dive; | |
734 //pVpm->vpm_conservatism = conservatism; | |
735 switch(conservatism) | |
736 { | |
737 case 0: | |
738 critical_radius_n2_microns=0.55; //!Adj. Range: 0.2 to 1.35 microns | |
739 critical_radius_he_microns=0.45; //!Adj. Range: 0.2 to 1.35 microns | |
740 break; | |
741 case 1: | |
742 critical_radius_n2_microns=0.58; | |
743 critical_radius_he_microns=0.48; | |
744 break; | |
745 case 2: | |
746 critical_radius_n2_microns=0.62; | |
747 critical_radius_he_microns=0.52; | |
748 break; | |
749 case 3: | |
750 critical_radius_n2_microns=0.68; | |
751 critical_radius_he_microns=0.58; | |
752 break; | |
753 case 4: | |
754 critical_radius_n2_microns=0.75; | |
755 critical_radius_he_microns=0.65; | |
756 break; | |
757 case 5: | |
758 critical_radius_n2_microns=0.82; | |
759 critical_radius_he_microns=0.72; | |
760 break; | |
51
8f8ea3a32e82
Resolved warnings pointing to possible invalid memory access
Ideenmodellierer
parents:
38
diff
changeset
|
761 default: |
8f8ea3a32e82
Resolved warnings pointing to possible invalid memory access
Ideenmodellierer
parents:
38
diff
changeset
|
762 critical_radius_n2_microns=0.82; |
8f8ea3a32e82
Resolved warnings pointing to possible invalid memory access
Ideenmodellierer
parents:
38
diff
changeset
|
763 critical_radius_he_microns=0.72; |
8f8ea3a32e82
Resolved warnings pointing to possible invalid memory access
Ideenmodellierer
parents:
38
diff
changeset
|
764 break; |
38 | 765 } |
766 | |
767 for (i = 0; i < 16; ++i) { | |
768 initial_critical_radius_n2[i] = critical_radius_n2_microns * 1e-6f; | |
769 initial_critical_radius_he[i] = critical_radius_he_microns * 1e-6f; | |
770 } | |
771 | |
772 | |
773 | |
774 if( (surface_time > 0) | |
775 && (!pVpm->repetitive_variables_not_valid) ) | |
776 //&& (pVpm->decomode_vpm_plus_conservatism_last_dive > 0) | |
777 //&& (pVpm->decomode_vpm_plus_conservatism_last_dive - 1 == pVpm->vpm_conservatism)) | |
778 { | |
779 vpm_repetitive_algorithm(pVpm, &surface_time,initial_critical_radius_he, initial_critical_radius_n2); | |
780 } | |
781 else | |
782 { | |
783 //Kein gültiger Wiederholungstauchgang | |
784 for (i = 0; i < 16; ++i) { | |
785 pVpm->adjusted_critical_radius_n2[i] = initial_critical_radius_n2[i]; | |
786 pVpm->adjusted_critical_radius_he[i] = initial_critical_radius_he[i]; | |
787 } | |
788 pVpm->repetitive_variables_not_valid = 0; | |
789 } | |
790 for (i = 0; i < 16; ++i) { | |
791 pVpm->max_crushing_pressure_he[i] = 0.0f; | |
792 pVpm->max_crushing_pressure_n2[i] = 0.0f; | |
793 pVpm->max_actual_gradient[i] = 0.0f; | |
794 pVpm->adjusted_crushing_pressure_he[i] = 0.0f; | |
795 pVpm->adjusted_crushing_pressure_n2[i] = 0.0f; | |
796 pVpm->initial_allowable_gradient_he[i] = 0.0f; | |
797 pVpm->initial_allowable_gradient_n2[i] = 0.0f; | |
798 } | |
799 pVpm->max_first_stop_depth_save = 0; | |
800 pVpm->depth_start_of_deco_zone_save = 0; | |
801 pVpm->run_time_start_of_deco_zone_save = 0; | |
802 pVpm->deco_zone_reached = 0; | |
803 } | |
804 | |
805 /* =============================================================================== */ | |
806 /* SUBROUTINE VPM_REPETITIVE_ALGORITHM */ | |
807 /* Purpose: This subprogram implements the VPM Repetitive Algorithm that was */ | |
808 /* envisioned by Professor David E. Yount only months before his passing. */ | |
809 /* =============================================================================== */ | |
810 int vpm_repetitive_algorithm(SVpm* pVpm, float *surface_interval_time, float* initial_critical_radius_he, float* initial_critical_radius_n2) | |
811 { | |
812 /* Local variables */ | |
813 static float max_actual_gradient_pascals; | |
814 //static float initial_allowable_grad_n2_pa, initial_allowable_grad_he_pa; | |
815 static short i; | |
816 static float adj_crush_pressure_n2_pascals, | |
817 new_critical_radius_n2, | |
818 adj_crush_pressure_he_pascals, | |
819 new_critical_radius_he; | |
820 | |
821 /* loop */ | |
822 /* =============================================================================== */ | |
823 /* CALCULATIONS */ | |
824 | |
825 /* by hw 160215 */ | |
826 /* IN: */ | |
827 /* pVpm->max_actual_gradient[i] */ | |
828 /* pVpm->initial_allowable_gradient_n2[i] */ | |
829 /* pVpm->initial_allowable_gradient_he[i] */ | |
830 /* pVpm->adjusted_crushing_pressure_he[i] */ | |
831 /* pVpm->adjusted_crushing_pressure_n2[i] */ | |
832 /* OUT: */ | |
833 /* pVpm->adjusted_critical_radius_n2[i] */ | |
834 /* pVpm->adjusted_critical_radius_he[i] */ | |
835 /* =============================================================================== */ | |
836 | |
837 for (i = 0; i < 16; ++i) { | |
838 max_actual_gradient_pascals = pVpm->max_actual_gradient[i] / UNITS_FACTOR * 101325.0f; | |
839 adj_crush_pressure_he_pascals = pVpm->adjusted_crushing_pressure_he[i] / UNITS_FACTOR * 101325.0f; | |
840 adj_crush_pressure_n2_pascals = pVpm->adjusted_crushing_pressure_n2[i] / UNITS_FACTOR * 101325.0f; | |
841 /* | |
842 initial_allowable_grad_he_pa = | |
843 pVpm->initial_allowable_gradient_he[i] / UNITS_FACTOR * 101325.0f; | |
844 initial_allowable_grad_n2_pa = | |
845 pVpm->initial_allowable_gradient_n2[i] / UNITS_FACTOR * 101325.0f; | |
846 */ | |
847 if (pVpm->max_actual_gradient[i] > pVpm->initial_allowable_gradient_n2[i]) | |
848 { | |
849 new_critical_radius_n2 = | |
850 SURFACE_TENSION_GAMMA * 2.0f * | |
851 (SKIN_COMPRESSION_GAMMAC - SURFACE_TENSION_GAMMA) / | |
852 (max_actual_gradient_pascals * SKIN_COMPRESSION_GAMMAC - | |
853 SURFACE_TENSION_GAMMA * adj_crush_pressure_n2_pascals); | |
854 | |
855 pVpm->adjusted_critical_radius_n2[i] = | |
856 initial_critical_radius_n2[i] | |
857 + (initial_critical_radius_n2[i] - new_critical_radius_n2) | |
858 * exp(-(*surface_interval_time) / REGENERATION_TIME_CONSTANT); | |
859 | |
860 } else { | |
861 pVpm->adjusted_critical_radius_n2[i] = | |
862 initial_critical_radius_n2[i]; | |
863 } | |
864 if (pVpm->max_actual_gradient[i] > pVpm->initial_allowable_gradient_he[i]) | |
865 { | |
866 new_critical_radius_he = | |
867 SURFACE_TENSION_GAMMA * 2.0f * | |
868 (SKIN_COMPRESSION_GAMMAC - SURFACE_TENSION_GAMMA) / | |
869 (max_actual_gradient_pascals * SKIN_COMPRESSION_GAMMAC - | |
870 SURFACE_TENSION_GAMMA * adj_crush_pressure_he_pascals); | |
871 | |
872 pVpm->adjusted_critical_radius_he[i] = | |
873 initial_critical_radius_he[i] | |
874 + ( initial_critical_radius_he[i] - new_critical_radius_he) | |
875 * exp(-(*surface_interval_time) / REGENERATION_TIME_CONSTANT); | |
876 } else { | |
877 pVpm->adjusted_critical_radius_he[i] = | |
878 initial_critical_radius_he[i]; | |
879 } | |
880 } | |
881 return 0; | |
882 } /* vpm_repetitive_algorithm */ |