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
+ − 727 float critical_radius_n2_microns = 0.82; /* be conservative in case of an unexpected parameter value */
+ − 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
+ − 761 default:
+ − 762 critical_radius_n2_microns=0.82;
+ − 763 critical_radius_he_microns=0.72;
+ − 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 */