Mercurial > public > ostc4
annotate Common/Src/calc_crush.c @ 310:95928ef3986f cleanup-4
Make dive mode detection more advanced
In commit a09b1855d656, a RTE function was factored out, that was
used to detect dive mode vs. surface mode.
In the hunt for the time difference bug between stopwatch and dive time
some progress was made, but its still not totally right. Add some old
logic back, that seems reasonable.
A pool test after this commit shows some improvement, but still there
is some random difference between stopwatch and dive time. Things like
perfectly in sync for 20 min. surface shortly, and descend again in
the same dive. Now, the was a 1 or 2 second difference.
Signed-off-by: Jan Mulder <jlmulder@xs4all.nl>
| author | Jan Mulder <jlmulder@xs4all.nl> |
|---|---|
| date | Sun, 26 May 2019 10:09:22 +0200 |
| 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 */ |