|
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
|
|
|
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 */
|
