Mercurial > public > ostc4
annotate Common/Src/calc_crush.c @ 767:14d6f8853b28
Update build for testing
| author | heinrichsweikamp |
|---|---|
| date | Wed, 12 Apr 2023 09:56:00 +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 */ |