Mercurial > public > ostc4

comparison Common/Src/calc_crush.c @ 38:5f11787b4f42

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
include in ostc4 repository
author heinrichsweikamp
date Sat, 28 Apr 2018 11:52:34 +0200
parents
children 8f8ea3a32e82
comparison
equal deleted inserted replaced
37:ccc45c0e1ea2 38:5f11787b4f42
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
727 float critical_radius_n2_microns,
728 critical_radius_he_microns;
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;
761 }
762
763 for (i = 0; i < 16; ++i) {
764 initial_critical_radius_n2[i] = critical_radius_n2_microns * 1e-6f;
765 initial_critical_radius_he[i] = critical_radius_he_microns * 1e-6f;
766 }
767
768
769
770 if( (surface_time > 0)
771 && (!pVpm->repetitive_variables_not_valid) )
772 //&& (pVpm->decomode_vpm_plus_conservatism_last_dive > 0)
773 //&& (pVpm->decomode_vpm_plus_conservatism_last_dive - 1 == pVpm->vpm_conservatism))
774 {
775 vpm_repetitive_algorithm(pVpm, &surface_time,initial_critical_radius_he, initial_critical_radius_n2);
776 }
777 else
778 {
779 //Kein gültiger Wiederholungstauchgang
780 for (i = 0; i < 16; ++i) {
781 pVpm->adjusted_critical_radius_n2[i] = initial_critical_radius_n2[i];
782 pVpm->adjusted_critical_radius_he[i] = initial_critical_radius_he[i];
783 }
784 pVpm->repetitive_variables_not_valid = 0;
785 }
786 for (i = 0; i < 16; ++i) {
787 pVpm->max_crushing_pressure_he[i] = 0.0f;
788 pVpm->max_crushing_pressure_n2[i] = 0.0f;
789 pVpm->max_actual_gradient[i] = 0.0f;
790 pVpm->adjusted_crushing_pressure_he[i] = 0.0f;
791 pVpm->adjusted_crushing_pressure_n2[i] = 0.0f;
792 pVpm->initial_allowable_gradient_he[i] = 0.0f;
793 pVpm->initial_allowable_gradient_n2[i] = 0.0f;
794 }
795 pVpm->max_first_stop_depth_save = 0;
796 pVpm->depth_start_of_deco_zone_save = 0;
797 pVpm->run_time_start_of_deco_zone_save = 0;
798 pVpm->deco_zone_reached = 0;
799 }
800
801 /* =============================================================================== */
802 /* SUBROUTINE VPM_REPETITIVE_ALGORITHM */
803 /* Purpose: This subprogram implements the VPM Repetitive Algorithm that was */
804 /* envisioned by Professor David E. Yount only months before his passing. */
805 /* =============================================================================== */
806 int vpm_repetitive_algorithm(SVpm* pVpm, float *surface_interval_time, float* initial_critical_radius_he, float* initial_critical_radius_n2)
807 {
808 /* Local variables */
809 static float max_actual_gradient_pascals;
810 //static float initial_allowable_grad_n2_pa, initial_allowable_grad_he_pa;
811 static short i;
812 static float adj_crush_pressure_n2_pascals,
813 new_critical_radius_n2,
814 adj_crush_pressure_he_pascals,
815 new_critical_radius_he;
816
817 /* loop */
818 /* =============================================================================== */
819 /* CALCULATIONS */
820
821 /* by hw 160215 */
822 /* IN: */
823 /* pVpm->max_actual_gradient[i] */
824 /* pVpm->initial_allowable_gradient_n2[i] */
825 /* pVpm->initial_allowable_gradient_he[i] */
826 /* pVpm->adjusted_crushing_pressure_he[i] */
827 /* pVpm->adjusted_crushing_pressure_n2[i] */
828 /* OUT: */
829 /* pVpm->adjusted_critical_radius_n2[i] */
830 /* pVpm->adjusted_critical_radius_he[i] */
831 /* =============================================================================== */
832
833 for (i = 0; i < 16; ++i) {
834 max_actual_gradient_pascals = pVpm->max_actual_gradient[i] / UNITS_FACTOR * 101325.0f;
835 adj_crush_pressure_he_pascals = pVpm->adjusted_crushing_pressure_he[i] / UNITS_FACTOR * 101325.0f;
836 adj_crush_pressure_n2_pascals = pVpm->adjusted_crushing_pressure_n2[i] / UNITS_FACTOR * 101325.0f;
837 /*
838 initial_allowable_grad_he_pa =
839 pVpm->initial_allowable_gradient_he[i] / UNITS_FACTOR * 101325.0f;
840 initial_allowable_grad_n2_pa =
841 pVpm->initial_allowable_gradient_n2[i] / UNITS_FACTOR * 101325.0f;
842 */
843 if (pVpm->max_actual_gradient[i] > pVpm->initial_allowable_gradient_n2[i])
844 {
845 new_critical_radius_n2 =
846 SURFACE_TENSION_GAMMA * 2.0f *
847 (SKIN_COMPRESSION_GAMMAC - SURFACE_TENSION_GAMMA) /
848 (max_actual_gradient_pascals * SKIN_COMPRESSION_GAMMAC -
849 SURFACE_TENSION_GAMMA * adj_crush_pressure_n2_pascals);
850
851 pVpm->adjusted_critical_radius_n2[i] =
852 initial_critical_radius_n2[i]
853 + (initial_critical_radius_n2[i] - new_critical_radius_n2)
854 * exp(-(*surface_interval_time) / REGENERATION_TIME_CONSTANT);
855
856 } else {
857 pVpm->adjusted_critical_radius_n2[i] =
858 initial_critical_radius_n2[i];
859 }
860 if (pVpm->max_actual_gradient[i] > pVpm->initial_allowable_gradient_he[i])
861 {
862 new_critical_radius_he =
863 SURFACE_TENSION_GAMMA * 2.0f *
864 (SKIN_COMPRESSION_GAMMAC - SURFACE_TENSION_GAMMA) /
865 (max_actual_gradient_pascals * SKIN_COMPRESSION_GAMMAC -
866 SURFACE_TENSION_GAMMA * adj_crush_pressure_he_pascals);
867
868 pVpm->adjusted_critical_radius_he[i] =
869 initial_critical_radius_he[i]
870 + ( initial_critical_radius_he[i] - new_critical_radius_he)
871 * exp(-(*surface_interval_time) / REGENERATION_TIME_CONSTANT);
872 } else {
873 pVpm->adjusted_critical_radius_he[i] =
874 initial_critical_radius_he[i];
875 }
876 }
877 return 0;
878 } /* vpm_repetitive_algorithm */