line source
///////////////////////////////////////////////////////////////////////////////+ −
/// -*- coding: UTF-8 -*-+ −
///+ −
/// \file Common/Src/calc_crush.c+ −
/// \brief VPM Desaturation code+ −
/// \author Heinrichs Weikamp+ −
/// \date 2018+ −
///+ −
/// $Id$+ −
///////////////////////////////////////////////////////////////////////////////+ −
/// \par Copyright (c) 2014-2018 Heinrichs Weikamp gmbh+ −
///+ −
/// This program is free software: you can redistribute it and/or modify+ −
/// it under the terms of the GNU General Public License as published by+ −
/// the Free Software Foundation, either version 3 of the License, or+ −
/// (at your option) any later version.+ −
///+ −
/// This program is distributed in the hope that it will be useful,+ −
/// but WITHOUT ANY WARRANTY; without even the implied warranty of+ −
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the+ −
/// GNU General Public License for more details.+ −
///+ −
/// You should have received a copy of the GNU General Public License+ −
/// along with this program. If not, see <http://www.gnu.org/licenses/>.+ −
//////////////////////////////////////////////////////////////////////////////+ −
+ −
#include "calc_crush.h"+ −
+ −
#include "decom.h"+ −
#include "math.h"+ −
#include "vpm.h"+ −
+ −
/* Common Block Declarations */+ −
//#pragma warning(disable:1035)+ −
+ −
const float SURFACE_TENSION_GAMMA = 0.0179f; //!Adj. Range: 0.015 to 0.065 N/m+ −
const float SKIN_COMPRESSION_GAMMAC = 0.257f; //!Adj. Range: 0.160 to 0.290 N/m+ −
const float UNITS_FACTOR = 10.1325f;+ −
const float WATER_VAPOR_PRESSURE = 0.493f; // (Schreiner value) based on respiratory quotien+ −
const float CRIT_VOLUME_PARAMETER_LAMBDA = 7500.0f; //!Adj. Range: 6500 to 8300 fsw-min+ −
const float GRADIENT_ONSET_OF_IMPERM_ATM = 8.2f; //!Adj. Range: 5.0 to 10.0 atm+ −
const float REGENERATION_TIME_CONSTANT = 20160.0f; //!Adj. Range: 10080 to 51840 min+ −
const float PRESSURE_OTHER_GASES_MMHG = 102.0f; //!Constant value for PO2 up to 2 atm+ −
const float CONSTANT_PRESSURE_OTHER_GASES = 102.0f * 10.1325f / 760.0f; // PRESSURE_OTHER_GASES_MMHG / 760. * UNITS_FACTOR;+ −
+ −
const float HELIUM_TIME_CONSTANT[16] = {3.68695308808482E-001f,+ −
2.29518933960247E-001f,+ −
1.46853216220327E-001f,+ −
9.91626867753856E-002f,+ −
6.78890480470074E-002f,+ −
4.78692804254106E-002f,+ −
3.37626488338989E-002f,+ −
2.38113081607676E-002f,+ −
1.68239606932026E-002f,+ −
1.25592893741610E-002f,+ −
9.80544886914621E-003f,+ −
7.67264977374303E-003f,+ −
6.01220557342307E-003f,+ −
4.70185307665137E-003f,+ −
3.68225234041620E-003f,+ −
2.88775228329769E-003f};+ −
+ −
const float NITROGEN_TIME_CONSTANT[16] = {1.38629436111989E-001f,+ −
8.66433975699932E-002f,+ −
5.54517744447956E-002f,+ −
3.74674151654024E-002f,+ −
2.56721177985165E-002f,+ −
1.80978376125312E-002f,+ −
1.27651414467762E-002f,+ −
9.00191143584345E-003f,+ −
6.35914844550409E-003f,+ −
4.74758342849278E-003f,+ −
3.70666941475907E-003f,+ −
2.90019740820061E-003f,+ −
2.27261370675392E-003f,+ −
1.77730046297422E-003f,+ −
1.39186180835330E-003f,+ −
1.09157036308653E-003f};+ −
+ −
int onset_of_impermeability(SGas* pGas, float *starting_ambient_pressure,+ −
float *ending_ambient_pressure,+ −
float *rate,+ −
float* amb_pressure_onset_of_imperm,+ −
float* gas_tension_onset_of_imperm,+ −
float* initial_helium_pressure,+ −
float* initial_nitrogen_pressure,+ −
short i);+ −
int radius_root_finder (float *a, float *b, float *c, float *low_bound, float *high_bound, float *ending_radius);+ −
//void get_inert_gases_(SBuehlmann* input, ,short gas_id, float ambient_pressure_bar, float* fraction_nitrogen,float* fraction_helium );+ −
int vpm_repetitive_algorithm(SVpm* pVpm, float *surface_interval_time, float* initial_critical_radius_he, float* initial_critical_radius_n2);+ −
+ −
+ −
+ −
/* =============================================================================== */+ −
/* NOTE ABOUT PRESSURE UNITS USED IN CALCULATIONS: */+ −
/* It is the convention in decompression calculations to compute all gas */+ −
/* loadings, absolute pressures, partial pressures, etc., in the units of */+ −
/* depth pressure that you are diving - either feet of seawater (fsw) or */+ −
/* meters of seawater (msw). This program follows that convention with the */+ −
/* the exception that all VPM calculations are performed in SI units (by */+ −
/* necessity). Accordingly, there are several conversions back and forth */+ −
/* between the diving pressure units and the SI units. */+ −
/* =============================================================================== */+ −
/* =============================================================================== */+ −
/* FUNCTION SUBPROGRAM FOR GAS LOADING CALCULATIONS - ASCENT AND DESCENT */+ −
/* =============================================================================== */+ −
+ −
float schreiner_equation__2(float *initial_inspired_gas_pressure,+ −
float *rate_change_insp_gas_pressure,+ −
float *interval_time_minutes,+ −
const float *gas_time_constant,+ −
float *initial_gas_pressure)+ −
{+ −
/* System generated locals */+ −
float ret_val;+ −
float time_null_pressure = 0.0f;+ −
float time_rest = 0.0f;+ −
float time = *interval_time_minutes;+ −
/* =============================================================================== */+ −
/* Note: The Schreiner equation is applied when calculating the uptake or */+ −
/* elimination of compartment gases during linear ascents or descents at a */+ −
/* constant rate. For ascents, a negative number for rate must be used. */+ −
/* =============================================================================== */+ −
if( *rate_change_insp_gas_pressure < 0.0f)+ −
{+ −
time_null_pressure = -1.0f * *initial_inspired_gas_pressure / *rate_change_insp_gas_pressure;+ −
if(time > time_null_pressure )+ −
{+ −
time_rest = time - time_null_pressure;+ −
time = time_null_pressure;+ −
}+ −
}+ −
ret_val =+ −
*initial_inspired_gas_pressure ++ −
*rate_change_insp_gas_pressure *+ −
(time - 1.f / *gas_time_constant) -+ −
(*initial_inspired_gas_pressure -+ −
*initial_gas_pressure -+ −
*rate_change_insp_gas_pressure / *gas_time_constant) *+ −
expf(-(*gas_time_constant) * time);+ −
+ −
if(time_rest > 0.0f)+ −
{+ −
ret_val = ret_val * expf(-(*gas_time_constant) * time_rest);+ −
}+ −
+ −
+ −
return ret_val;+ −
}; /* schreiner_equation__2 */+ −
+ −
/* =============================================================================== */+ −
/* SUBROUTINE CALC_CRUSHING_PRESSURE */+ −
/* Purpose: Compute the effective "crushing pressure" in each compartment as */+ −
/* a result of descent segment(s). The crushing pressure is the gradient */+ −
/* (difference in pressure) between the outside ambient pressure and the */+ −
/* gas tension inside a VPM nucleus (bubble seed). This gradient acts to */+ −
/* reduce (shrink) the radius smaller than its initial value at the surface. */+ −
/* This phenomenon has important ramifications because the smaller the radius */+ −
/* of a VPM nucleus, the greater the allowable supersaturation gradient upon */+ −
/* ascent. Gas loading (uptake) during descent, especially in the fast */+ −
/* compartments, will reduce the magnitude of the crushing pressure. The */+ −
/* crushing pressure is not cumulative over a multi-level descent. It will */+ −
/* be the maximum value obtained in any one discrete segment of the overall */+ −
/* descent. Thus, the program must compute and store the maximum crushing */+ −
/* pressure for each compartment that was obtained across all segments of */+ −
/* the descent profile. */+ −
+ −
/* The calculation of crushing pressure will be different depending on */+ −
/* whether or not the gradient is in the VPM permeable range (gas can diffuse */+ −
/* across skin of VPM nucleus) or the VPM impermeable range (molecules in */+ −
/* skin of nucleus are squeezed together so tight that gas can no longer */+ −
/* diffuse in or out of nucleus; the gas becomes trapped and further resists */+ −
/* the crushing pressure). The solution for crushing pressure in the VPM */+ −
/* permeable range is a simple linear equation. In the VPM impermeable */+ −
/* range, a cubic equation must be solved using a numerical method. */+ −
+ −
/* Separate crushing pressures are tracked for helium and nitrogen because */+ −
/* they can have different critical radii. The crushing pressures will be */+ −
/* the same for helium and nitrogen in the permeable range of the model, but */+ −
/* they will start to diverge in the impermeable range. This is due to */+ −
/* the differences between starting radius, radius at the onset of */+ −
/* impermeability, and radial compression in the impermeable range. */+ −
/* =============================================================================== */+ −
int calc_crushing_pressure(SLifeData* lifeData, SVpm* vpm, float * initial_helium_pressure, float * initial_nitrogen_pressure, float starting_ambient_pressure,+ −
float rate )+ −
{+ −
/* System generated locals */+ −
static float r1, r2;+ −
+ −
static float low_bound_n2,+ −
ending_radius_n2,+ −
gradient_onset_of_imperm_pa;+ −
static float low_bound_he,+ −
ending_radius_he,+ −
high_bound_n2,+ −
crushing_pressure_n2;+ −
short i;+ −
static float crushing_pressure_pascals_n2,+ −
gradient_onset_of_imperm,+ −
starting_gas_tension,+ −
high_bound_he,+ −
crushing_pressure_he,+ −
amb_press_onset_of_imperm_pa,+ −
crushing_pressure_pascals_he,+ −
radius_onset_of_imperm_n2,+ −
starting_gradient,+ −
radius_onset_of_imperm_he,+ −
ending_gas_tension;+ −
+ −
static float ending_ambient_pressure_pa,+ −
a_n2,+ −
b_n2,+ −
c_n2,+ −
ending_gradient,+ −
gas_tension_onset_of_imperm_pa,+ −
a_he,+ −
b_he, c_he;+ −
static float amb_pressure_onset_of_imperm[16];+ −
static float gas_tension_onset_of_imperm[16];+ −
+ −
static float helium_pressure_crush[16];+ −
static float nitrogen_pressure_crush[16];+ −
+ −
+ −
static float ending_ambient_pressure = 0;+ −
+ −
ending_ambient_pressure = lifeData->pressure_ambient_bar * 10;+ −
for( i = 0; i < 16; i++)+ −
{+ −
helium_pressure_crush[i] = lifeData->tissue_helium_bar[i] * 10;+ −
nitrogen_pressure_crush[i] = lifeData->tissue_nitrogen_bar[i] * 10;+ −
}+ −
+ −
+ −
+ −
+ −
+ −
/* loop */+ −
/* =============================================================================== */+ −
/* CALCULATIONS */+ −
/* First, convert the Gradient for Onset of Impermeability from units of */+ −
/* atmospheres to diving pressure units (either fsw or msw) and to Pascals */+ −
/* (SI units). The reason that the Gradient for Onset of Impermeability is */+ −
/* given in the program settings in units of atmospheres is because that is */+ −
/* how it was reported in the original research papers by Yount and */+ −
/* colleauges. */+ −
/* =============================================================================== */+ −
+ −
gradient_onset_of_imperm = GRADIENT_ONSET_OF_IMPERM_ATM * UNITS_FACTOR;+ −
gradient_onset_of_imperm_pa = GRADIENT_ONSET_OF_IMPERM_ATM * 101325.0f;+ −
+ −
/* =============================================================================== */+ −
/* Assign values of starting and ending ambient pressures for descent segment */+ −
/* =============================================================================== */+ −
+ −
//starting_ambient_pressure = *starting_depth;+ −
//ending_ambient_pressure = *ending_depth;+ −
+ −
/* =============================================================================== */+ −
/* MAIN LOOP WITH NESTED DECISION TREE */+ −
/* For each compartment, the program computes the starting and ending */+ −
/* gas tensions and gradients. The VPM is different than some dissolved gas */+ −
/* algorithms, Buhlmann for example, in that it considers the pressure due to */+ −
/* oxygen, carbon dioxide, and water vapor in each compartment in addition to */+ −
/* the inert gases helium and nitrogen. These "other gases" are included in */+ −
/* the calculation of gas tensions and gradients. */+ −
/* =============================================================================== */+ −
+ −
crushing_pressure_he = 0.0f;+ −
crushing_pressure_n2 = 0.0f;+ −
+ −
for (i = 0; i < 16; ++i) {+ −
starting_gas_tension = initial_helium_pressure[i] + initial_nitrogen_pressure[i] + CONSTANT_PRESSURE_OTHER_GASES;+ −
starting_gradient = starting_ambient_pressure - starting_gas_tension;+ −
ending_gas_tension = helium_pressure_crush[i] + nitrogen_pressure_crush[i] + CONSTANT_PRESSURE_OTHER_GASES;+ −
ending_gradient = ending_ambient_pressure - ending_gas_tension;+ −
+ −
/* =============================================================================== */+ −
/* Compute radius at onset of impermeability for helium and nitrogen */+ −
/* critical radii */+ −
/* =============================================================================== */+ −
+ −
radius_onset_of_imperm_he = 1.0f / (+ −
gradient_onset_of_imperm_pa /+ −
((SKIN_COMPRESSION_GAMMAC -+ −
SURFACE_TENSION_GAMMA) * 2.0f) ++ −
1.0f / vpm->adjusted_critical_radius_he[i]);+ −
radius_onset_of_imperm_n2 = 1.0f / (+ −
gradient_onset_of_imperm_pa /+ −
((SKIN_COMPRESSION_GAMMAC -+ −
SURFACE_TENSION_GAMMA) * 2.0f) ++ −
1.0f / vpm->adjusted_critical_radius_n2[i]);+ −
+ −
/* =============================================================================== */+ −
/* FIRST BRANCH OF DECISION TREE - PERMEABLE RANGE */+ −
/* Crushing pressures will be the same for helium and nitrogen */+ −
/* =============================================================================== */+ −
+ −
if (ending_gradient <= gradient_onset_of_imperm) {+ −
crushing_pressure_he = ending_ambient_pressure - ending_gas_tension;+ −
crushing_pressure_n2 = ending_ambient_pressure - ending_gas_tension;+ −
}+ −
+ −
/* =============================================================================== */+ −
/* SECOND BRANCH OF DECISION TREE - IMPERMEABLE RANGE */+ −
/* Both the ambient pressure and the gas tension at the onset of */+ −
/* impermeability must be computed in order to properly solve for the ending */+ −
/* radius and resultant crushing pressure. The first decision block */+ −
/* addresses the special case when the starting gradient just happens to be */+ −
/* equal to the gradient for onset of impermeability (not very likely!). */+ −
/* =============================================================================== */+ −
+ −
if (ending_gradient > gradient_onset_of_imperm) {+ −
if (starting_gradient == gradient_onset_of_imperm) {+ −
amb_pressure_onset_of_imperm[i] = starting_ambient_pressure;+ −
gas_tension_onset_of_imperm[i] = starting_gas_tension;+ −
}+ −
+ −
/* =============================================================================== */+ −
/* In most cases, a subroutine will be called to find these values using a */+ −
/* numerical method. */+ −
/* =============================================================================== */+ −
+ −
if (starting_gradient < gradient_onset_of_imperm) {+ −
onset_of_impermeability(&(lifeData->actualGas), &starting_ambient_pressure, &ending_ambient_pressure, &rate,+ −
amb_pressure_onset_of_imperm, gas_tension_onset_of_imperm,+ −
initial_helium_pressure, initial_nitrogen_pressure, i);+ −
}+ −
+ −
/* =============================================================================== */+ −
/* Next, using the values for ambient pressure and gas tension at the onset */+ −
/* of impermeability, the equations are set up to process the calculations */+ −
/* through the radius root finder subroutine. This subprogram will find the */+ −
/* root (solution) to the cubic equation using a numerical method. In order */+ −
/* to do this efficiently, the equations are placed in the form */+ −
/* Ar^3 - Br^2 - C = 0, where r is the ending radius after impermeable */+ −
/* compression. The coefficients A, B, and C for helium and nitrogen are */+ −
/* computed and passed to the subroutine as arguments. The high and low */+ −
/* bounds to be used by the numerical method of the subroutine are also */+ −
/* computed (see separate page posted on Deco List ftp site entitled */+ −
/* "VPM: Solving for radius in the impermeable regime"). The subprogram */+ −
/* will return the value of the ending radius and then the crushing */+ −
/* pressures for helium and nitrogen can be calculated. */+ −
/* =============================================================================== */+ −
+ −
ending_ambient_pressure_pa = ending_ambient_pressure / UNITS_FACTOR * 101325.0f;+ −
amb_press_onset_of_imperm_pa = amb_pressure_onset_of_imperm[i] / UNITS_FACTOR * 101325.0f;+ −
gas_tension_onset_of_imperm_pa = gas_tension_onset_of_imperm[i] / UNITS_FACTOR * 101325.0f;+ −
b_he = (SKIN_COMPRESSION_GAMMAC - SURFACE_TENSION_GAMMA) * 2.0f;+ −
a_he = ending_ambient_pressure_pa - amb_press_onset_of_imperm_pa + gas_tension_onset_of_imperm_pa+ −
+ (SKIN_COMPRESSION_GAMMAC - SURFACE_TENSION_GAMMA) * 2.0f / radius_onset_of_imperm_he;+ −
/* Computing 3rd power */+ −
r1 = radius_onset_of_imperm_he;+ −
c_he = gas_tension_onset_of_imperm_pa * (r1 * (r1 * r1));+ −
high_bound_he = radius_onset_of_imperm_he;+ −
low_bound_he = b_he / a_he;+ −
radius_root_finder(&a_he, &b_he, &c_he, &low_bound_he, &high_bound_he, &ending_radius_he);+ −
/* Computing 3rd power */+ −
r1 = radius_onset_of_imperm_he;+ −
/* Computing 3rd power */+ −
r2 = ending_radius_he;+ −
crushing_pressure_pascals_he =+ −
gradient_onset_of_imperm_pa ++ −
ending_ambient_pressure_pa -+ −
amb_press_onset_of_imperm_pa ++ −
gas_tension_onset_of_imperm_pa *+ −
(1.0f - r1 * (r1 * r1) / (r2 * (r2 * r2)));+ −
crushing_pressure_he =+ −
crushing_pressure_pascals_he / 101325.0f * UNITS_FACTOR;+ −
b_n2 = (SKIN_COMPRESSION_GAMMAC - SURFACE_TENSION_GAMMA) * 2.0f;+ −
a_n2 = ending_ambient_pressure_pa -+ −
amb_press_onset_of_imperm_pa ++ −
gas_tension_onset_of_imperm_pa ++ −
(SKIN_COMPRESSION_GAMMAC - SURFACE_TENSION_GAMMA) *+ −
2.0f / radius_onset_of_imperm_n2;+ −
/* Computing 3rd power */+ −
r1 = radius_onset_of_imperm_n2;+ −
c_n2 = gas_tension_onset_of_imperm_pa * (r1 * (r1 * r1));+ −
high_bound_n2 = radius_onset_of_imperm_n2;+ −
low_bound_n2 = b_n2 / a_n2;+ −
radius_root_finder(&a_n2,+ −
&b_n2,+ −
&c_n2,+ −
&low_bound_n2,+ −
&high_bound_n2,+ −
&ending_radius_n2);+ −
+ −
/* Computing 3rd power */+ −
r1 = radius_onset_of_imperm_n2;+ −
/* Computing 3rd power */+ −
r2 = ending_radius_n2;+ −
crushing_pressure_pascals_n2 =+ −
gradient_onset_of_imperm_pa ++ −
ending_ambient_pressure_pa -+ −
amb_press_onset_of_imperm_pa ++ −
gas_tension_onset_of_imperm_pa * (1.0f - r1 *+ −
(r1 * r1) / (r2 * (r2 * r2)));+ −
crushing_pressure_n2 = crushing_pressure_pascals_n2 / 101325.0f * UNITS_FACTOR;+ −
}+ −
+ −
/* =============================================================================== */+ −
/* UPDATE VALUES OF MAX CRUSHING PRESSURE IN GLOBAL ARRAYS */+ −
/* =============================================================================== */+ −
+ −
/* Computing MAX */+ −
r1 = vpm->max_crushing_pressure_he[i];+ −
vpm->max_crushing_pressure_he[i] = fmaxf(r1, crushing_pressure_he);+ −
/* Computing MAX */+ −
r1 = vpm->max_crushing_pressure_n2[i];+ −
vpm->max_crushing_pressure_n2[i] = fmaxf(r1, crushing_pressure_n2);+ −
}+ −
return 0;+ −
} /* calc_crushing_pressure */+ −
+ −
/* =============================================================================== */+ −
/* SUBROUTINE ONSET_OF_IMPERMEABILITY */+ −
/* Purpose: This subroutine uses the Bisection Method to find the ambient */+ −
/* pressure and gas tension at the onset of impermeability for a given */+ −
/* compartment. Source: "Numerical Recipes in Fortran 77", */+ −
/* Cambridge University Press, 1992. */+ −
/* =============================================================================== */+ −
+ −
int onset_of_impermeability(SGas* pGas, float *starting_ambient_pressure,+ −
float *ending_ambient_pressure,+ −
float *rate,+ −
float* amb_pressure_onset_of_imperm,+ −
float* gas_tension_onset_of_imperm,+ −
float* initial_helium_pressure,+ −
float* initial_nitrogen_pressure,+ −
short i)+ −
{+ −
/* Local variables */+ −
float time, last_diff_change, mid_range_nitrogen_pressure;+ −
short j;+ −
float gas_tension_at_mid_range,+ −
initial_inspired_n2_pressure,+ −
gradient_onset_of_imperm,+ −
starting_gas_tension,+ −
low_bound,+ −
initial_inspired_he_pressure,+ −
high_bound_nitrogen_pressure,+ −
nitrogen_rate,+ −
function_at_mid_range,+ −
function_at_low_bound,+ −
high_bound,+ −
mid_range_helium_pressure,+ −
mid_range_time,+ −
ending_gas_tension,+ −
function_at_high_bound;+ −
+ −
float mid_range_ambient_pressure,+ −
high_bound_helium_pressure,+ −
helium_rate,+ −
differential_change;+ −
float fraction_helium_begin;+ −
float fraction_helium_end;+ −
float fraction_nitrogen_begin;+ −
float fraction_nitrogen_end;+ −
/* loop */+ −
/* =============================================================================== */+ −
/* CALCULATIONS */+ −
/* First convert the Gradient for Onset of Impermeability to the diving */+ −
/* pressure units that are being used */+ −
/* =============================================================================== */+ −
+ −
gradient_onset_of_imperm = GRADIENT_ONSET_OF_IMPERM_ATM * UNITS_FACTOR;+ −
+ −
/* =============================================================================== */+ −
/* ESTABLISH THE BOUNDS FOR THE ROOT SEARCH USING THE BISECTION METHOD */+ −
/* In this case, we are solving for time - the time when the ambient pressure */+ −
/* minus the gas tension will be equal to the Gradient for Onset of */+ −
/* Impermeabliity. The low bound for time is set at zero and the high */+ −
/* bound is set at the elapsed time (segment time) it took to go from the */+ −
/* starting ambient pressure to the ending ambient pressure. The desired */+ −
/* ambient pressure and gas tension at the onset of impermeability will */+ −
/* be found somewhere between these endpoints. The algorithm checks to */+ −
/* make sure that the solution lies in between these bounds by first */+ −
/* computing the low bound and high bound function values. */+ −
/* =============================================================================== */+ −
+ −
/*initial_inspired_he_pressure =+ −
(*starting_ambient_pressure - water_vapor_pressure) * fraction_helium[mix_number - 1];+ −
initial_inspired_n2_pressure =+ −
(*starting_ambient_pressure - water_vapor_pressure) * fraction_nitrogen[mix_number - 1];+ −
helium_rate = *rate * fraction_helium[mix_number - 1];+ −
nitrogen_rate = *rate * fraction_nitrogen[mix_number - 1];*/+ −
low_bound = 0.;+ −
high_bound = (*ending_ambient_pressure - *starting_ambient_pressure) / *rate;+ −
+ −
//New+ −
decom_get_inert_gases( *starting_ambient_pressure / 10.0f, pGas, &fraction_nitrogen_begin, &fraction_helium_begin );+ −
decom_get_inert_gases(*ending_ambient_pressure / 10.0f, pGas, &fraction_nitrogen_end, &fraction_helium_end );+ −
initial_inspired_he_pressure = (*starting_ambient_pressure - WATER_VAPOR_PRESSURE) * fraction_helium_begin;+ −
initial_inspired_n2_pressure = (*starting_ambient_pressure - WATER_VAPOR_PRESSURE) * fraction_nitrogen_begin;+ −
helium_rate = ((*ending_ambient_pressure - WATER_VAPOR_PRESSURE)* fraction_helium_end - initial_inspired_he_pressure)/high_bound;+ −
nitrogen_rate = ((*ending_ambient_pressure - WATER_VAPOR_PRESSURE)* fraction_nitrogen_end - initial_inspired_n2_pressure)/high_bound;+ −
+ −
starting_gas_tension =+ −
initial_helium_pressure[i] ++ −
initial_nitrogen_pressure[i] ++ −
CONSTANT_PRESSURE_OTHER_GASES;+ −
function_at_low_bound =+ −
*starting_ambient_pressure -+ −
starting_gas_tension -+ −
gradient_onset_of_imperm;+ −
high_bound_helium_pressure =+ −
schreiner_equation__2(&initial_inspired_he_pressure,+ −
&helium_rate,+ −
&high_bound,+ −
&HELIUM_TIME_CONSTANT[i],+ −
&initial_helium_pressure[i]);+ −
high_bound_nitrogen_pressure =+ −
schreiner_equation__2(&initial_inspired_n2_pressure,+ −
&nitrogen_rate,+ −
&high_bound,+ −
&NITROGEN_TIME_CONSTANT[i],+ −
&initial_nitrogen_pressure[i]);+ −
ending_gas_tension =+ −
high_bound_helium_pressure ++ −
high_bound_nitrogen_pressure ++ −
CONSTANT_PRESSURE_OTHER_GASES;+ −
function_at_high_bound =+ −
*ending_ambient_pressure -+ −
ending_gas_tension -+ −
gradient_onset_of_imperm;+ −
if (function_at_high_bound * function_at_low_bound >= 0.0f) {+ −
//printf("\nERROR! ROOT IS NOT WITHIN BRACKETS");+ −
}+ −
+ −
/* =============================================================================== */+ −
/* APPLY THE BISECTION METHOD IN SEVERAL ITERATIONS UNTIL A SOLUTION WITH */+ −
/* THE DESIRED ACCURACY IS FOUND */+ −
/* Note: the program allows for up to 100 iterations. Normally an exit will */+ −
/* be made from the loop well before that number. If, for some reason, the */+ −
/* program exceeds 100 iterations, there will be a pause to alert the user. */+ −
/* =============================================================================== */+ −
+ −
if (function_at_low_bound < 0.0f) {+ −
time = low_bound;+ −
differential_change = high_bound - low_bound;+ −
} else {+ −
time = high_bound;+ −
differential_change = low_bound - high_bound;+ −
}+ −
for (j = 1; j <= 100; ++j) {+ −
last_diff_change = differential_change;+ −
differential_change = last_diff_change * 0.5f;+ −
mid_range_time = time + differential_change;+ −
mid_range_ambient_pressure = *starting_ambient_pressure + *rate * mid_range_time;+ −
mid_range_helium_pressure =+ −
schreiner_equation__2(&initial_inspired_he_pressure,+ −
&helium_rate,+ −
&mid_range_time,+ −
&HELIUM_TIME_CONSTANT[i],+ −
&initial_helium_pressure[i]);+ −
mid_range_nitrogen_pressure =+ −
schreiner_equation__2(&initial_inspired_n2_pressure,+ −
&nitrogen_rate,+ −
&mid_range_time,+ −
&NITROGEN_TIME_CONSTANT[i],+ −
&initial_nitrogen_pressure[i]);+ −
gas_tension_at_mid_range =+ −
mid_range_helium_pressure ++ −
mid_range_nitrogen_pressure ++ −
CONSTANT_PRESSURE_OTHER_GASES;+ −
function_at_mid_range =+ −
mid_range_ambient_pressure -+ −
gas_tension_at_mid_range -+ −
gradient_onset_of_imperm;+ −
if (function_at_mid_range <= 0.0f) {+ −
time = mid_range_time;+ −
}+ −
if (fabs(differential_change) < .001f ||+ −
function_at_mid_range == 0.0f) {+ −
goto L100;+ −
}+ −
}+ −
//printf("\nERROR! ROOT SEARCH EXCEEDED MAXIMUM ITERATIONS");+ −
+ −
/* =============================================================================== */+ −
/* When a solution with the desired accuracy is found, the program jumps out */+ −
/* of the loop to Line 100 and assigns the solution values for ambient */+ −
/* pressure and gas tension at the onset of impermeability. */+ −
/* =============================================================================== */+ −
+ −
L100:+ −
amb_pressure_onset_of_imperm[i] = mid_range_ambient_pressure;+ −
gas_tension_onset_of_imperm[i] = gas_tension_at_mid_range;+ −
return 0;+ −
} /* onset_of_impermeability */+ −
+ −
+ −
/* =============================================================================== */+ −
/* SUBROUTINE RADIUS_ROOT_FINDER */+ −
/* Purpose: This subroutine is a "fail-safe" routine that combines the */+ −
/* Bisection Method and the Newton-Raphson Method to find the desired root. */+ −
/* This hybrid algorithm takes a bisection step whenever Newton-Raphson would */+ −
/* take the solution out of bounds, or whenever Newton-Raphson is not */+ −
/* converging fast enough. Source: "Numerical Recipes in Fortran 77", */+ −
/* Cambridge University Press, 1992. */+ −
/* =============================================================================== */+ −
+ −
int radius_root_finder (float *a,+ −
float *b,+ −
float *c,+ −
float *low_bound,+ −
float *high_bound,+ −
float *ending_radius)+ −
{+ −
/* System generated locals */+ −
float r1, r2;+ −
+ −
/* Local variables */+ −
float radius_at_low_bound,+ −
last_diff_change,+ −
function,+ −
radius_at_high_bound;+ −
short i;+ −
float function_at_low_bound,+ −
last_ending_radius,+ −
function_at_high_bound,+ −
derivative_of_function,+ −
differential_change;+ −
+ −
/* loop */+ −
/* =============================================================================== */+ −
/* BEGIN CALCULATIONS BY MAKING SURE THAT THE ROOT LIES WITHIN BOUNDS */+ −
/* In this case we are solving for radius in a cubic equation of the form, */+ −
/* Ar^3 - Br^2 - C = 0. The coefficients A, B, and C were passed to this */+ −
/* subroutine as arguments. */+ −
/* =============================================================================== */+ −
+ −
function_at_low_bound =+ −
*low_bound * (*low_bound * (*a * *low_bound - *b)) - *c;+ −
function_at_high_bound =+ −
*high_bound * (*high_bound * (*a * *high_bound - *b)) - *c;+ −
if (function_at_low_bound > 0.0f && function_at_high_bound > 0.0f) {+ −
// printf("\nERROR! ROOT IS NOT WITHIN BRACKETS");+ −
+ −
}+ −
+ −
/* =============================================================================== */+ −
/* Next the algorithm checks for special conditions and then prepares for */+ −
/* the first bisection. */+ −
/* =============================================================================== */+ −
+ −
if (function_at_low_bound < 0.0f && function_at_high_bound < 0.0f) {+ −
//printf("\nERROR! ROOT IS NOT WITHIN BRACKETS");+ −
+ −
}+ −
if (function_at_low_bound == 0.0f) {+ −
*ending_radius = *low_bound;+ −
return 0;+ −
} else if (function_at_high_bound == 0.0f) {+ −
*ending_radius = *high_bound;+ −
return 0;+ −
} else if (function_at_low_bound < 0.0f) {+ −
radius_at_low_bound = *low_bound;+ −
radius_at_high_bound = *high_bound;+ −
} else {+ −
radius_at_high_bound = *low_bound;+ −
radius_at_low_bound = *high_bound;+ −
}+ −
*ending_radius = (*low_bound + *high_bound) * .5f;+ −
last_diff_change = (r1 = *high_bound - *low_bound, fabs(r1));+ −
differential_change = last_diff_change;+ −
+ −
/* =============================================================================== */+ −
/* At this point, the Newton-Raphson Method is applied which uses a function */+ −
/* and its first derivative to rapidly converge upon a solution. */+ −
/* Note: the program allows for up to 100 iterations. Normally an exit will */+ −
/* be made from the loop well before that number. If, for some reason, the */+ −
/* program exceeds 100 iterations, there will be a pause to alert the user. */+ −
/* When a solution with the desired accuracy is found, exit is made from the */+ −
/* loop by returning to the calling program. The last value of ending */+ −
/* radius has been assigned as the solution. */+ −
/* =============================================================================== */+ −
+ −
function =+ −
*ending_radius * (*ending_radius * (*a * *ending_radius - *b)) - *c;+ −
derivative_of_function =+ −
*ending_radius * (*ending_radius * 3.0f * *a - *b * 2.0f);+ −
for (i = 1; i <= 100; ++i) {+ −
if (((*ending_radius - radius_at_high_bound) * derivative_of_function - function) *+ −
((*ending_radius - radius_at_low_bound) * derivative_of_function - function) >= 0.0f+ −
|| (r1 = function * 2.0f, fabs(r1)) >+ −
(r2 = last_diff_change * derivative_of_function, fabs(r2))) {+ −
last_diff_change = differential_change;+ −
differential_change =+ −
(radius_at_high_bound - radius_at_low_bound) * .5f;+ −
*ending_radius = radius_at_low_bound + differential_change;+ −
if (radius_at_low_bound == *ending_radius) {+ −
return 0;+ −
}+ −
} else {+ −
last_diff_change = differential_change;+ −
differential_change = function / derivative_of_function;+ −
last_ending_radius = *ending_radius;+ −
*ending_radius -= differential_change;+ −
if (last_ending_radius == *ending_radius) {+ −
return 0;+ −
}+ −
}+ −
if (fabs(differential_change) < 1e-12) {+ −
return 0;+ −
}+ −
function =+ −
*ending_radius * (*ending_radius * (*a * *ending_radius - *b)) - *c;+ −
derivative_of_function =+ −
*ending_radius * (*ending_radius * 3.0f * *a - *b * 2.0f);+ −
if (function < 0.0f) {+ −
radius_at_low_bound = *ending_radius;+ −
} else {+ −
radius_at_high_bound = *ending_radius;+ −
}+ −
}+ −
// printf("\nERROR! ROOT SEARCH EXCEEDED MAXIMUM ITERATIONS");+ −
return 0;+ −
} /* radius_root_finder */+ −
+ −
+ −
+ −
+ −
void vpm_init(SVpm* pVpm, short conservatism, short repetitive_dive, long seconds_since_last_dive)+ −
{+ −
+ −
float critical_radius_n2_microns = 0.82; /* be conservative in case of an unexpected parameter value */+ −
float critical_radius_he_microns = 0.72;+ −
float initial_critical_radius_n2[16];+ −
float initial_critical_radius_he[16];+ −
int i = 0;+ −
float surface_time = seconds_since_last_dive / 60;+ −
pVpm->repetitive_variables_not_valid = !repetitive_dive;+ −
//pVpm->vpm_conservatism = conservatism;+ −
switch(conservatism)+ −
{+ −
case 0:+ −
critical_radius_n2_microns=0.55; //!Adj. Range: 0.2 to 1.35 microns+ −
critical_radius_he_microns=0.45; //!Adj. Range: 0.2 to 1.35 microns+ −
break;+ −
case 1:+ −
critical_radius_n2_microns=0.58;+ −
critical_radius_he_microns=0.48;+ −
break;+ −
case 2:+ −
critical_radius_n2_microns=0.62;+ −
critical_radius_he_microns=0.52;+ −
break;+ −
case 3:+ −
critical_radius_n2_microns=0.68;+ −
critical_radius_he_microns=0.58;+ −
break;+ −
case 4:+ −
critical_radius_n2_microns=0.75;+ −
critical_radius_he_microns=0.65;+ −
break;+ −
case 5:+ −
critical_radius_n2_microns=0.82;+ −
critical_radius_he_microns=0.72;+ −
break;+ −
default:+ −
critical_radius_n2_microns=0.82;+ −
critical_radius_he_microns=0.72;+ −
break;+ −
}+ −
+ −
for (i = 0; i < 16; ++i) {+ −
initial_critical_radius_n2[i] = critical_radius_n2_microns * 1e-6f;+ −
initial_critical_radius_he[i] = critical_radius_he_microns * 1e-6f;+ −
}+ −
+ −
+ −
+ −
if( (surface_time > 0)+ −
&& (!pVpm->repetitive_variables_not_valid) )+ −
//&& (pVpm->decomode_vpm_plus_conservatism_last_dive > 0)+ −
//&& (pVpm->decomode_vpm_plus_conservatism_last_dive - 1 == pVpm->vpm_conservatism))+ −
{+ −
vpm_repetitive_algorithm(pVpm, &surface_time,initial_critical_radius_he, initial_critical_radius_n2);+ −
}+ −
else+ −
{+ −
//Kein gültiger Wiederholungstauchgang+ −
for (i = 0; i < 16; ++i) {+ −
pVpm->adjusted_critical_radius_n2[i] = initial_critical_radius_n2[i];+ −
pVpm->adjusted_critical_radius_he[i] = initial_critical_radius_he[i];+ −
}+ −
pVpm->repetitive_variables_not_valid = 0;+ −
}+ −
for (i = 0; i < 16; ++i) {+ −
pVpm->max_crushing_pressure_he[i] = 0.0f;+ −
pVpm->max_crushing_pressure_n2[i] = 0.0f;+ −
pVpm->max_actual_gradient[i] = 0.0f;+ −
pVpm->adjusted_crushing_pressure_he[i] = 0.0f;+ −
pVpm->adjusted_crushing_pressure_n2[i] = 0.0f;+ −
pVpm->initial_allowable_gradient_he[i] = 0.0f;+ −
pVpm->initial_allowable_gradient_n2[i] = 0.0f;+ −
}+ −
pVpm->max_first_stop_depth_save = 0;+ −
pVpm->depth_start_of_deco_zone_save = 0;+ −
pVpm->run_time_start_of_deco_zone_save = 0;+ −
pVpm->deco_zone_reached = 0;+ −
}+ −
+ −
/* =============================================================================== */+ −
/* SUBROUTINE VPM_REPETITIVE_ALGORITHM */+ −
/* Purpose: This subprogram implements the VPM Repetitive Algorithm that was */+ −
/* envisioned by Professor David E. Yount only months before his passing. */+ −
/* =============================================================================== */+ −
int vpm_repetitive_algorithm(SVpm* pVpm, float *surface_interval_time, float* initial_critical_radius_he, float* initial_critical_radius_n2)+ −
{+ −
/* Local variables */+ −
static float max_actual_gradient_pascals;+ −
//static float initial_allowable_grad_n2_pa, initial_allowable_grad_he_pa;+ −
static short i;+ −
static float adj_crush_pressure_n2_pascals,+ −
new_critical_radius_n2,+ −
adj_crush_pressure_he_pascals,+ −
new_critical_radius_he;+ −
+ −
/* loop */+ −
/* =============================================================================== */+ −
/* CALCULATIONS */+ −
+ −
/* by hw 160215 */+ −
/* IN: */+ −
/* pVpm->max_actual_gradient[i] */+ −
/* pVpm->initial_allowable_gradient_n2[i] */+ −
/* pVpm->initial_allowable_gradient_he[i] */+ −
/* pVpm->adjusted_crushing_pressure_he[i] */+ −
/* pVpm->adjusted_crushing_pressure_n2[i] */+ −
/* OUT: */+ −
/* pVpm->adjusted_critical_radius_n2[i] */+ −
/* pVpm->adjusted_critical_radius_he[i] */+ −
/* =============================================================================== */+ −
+ −
for (i = 0; i < 16; ++i) {+ −
max_actual_gradient_pascals = pVpm->max_actual_gradient[i] / UNITS_FACTOR * 101325.0f;+ −
adj_crush_pressure_he_pascals = pVpm->adjusted_crushing_pressure_he[i] / UNITS_FACTOR * 101325.0f;+ −
adj_crush_pressure_n2_pascals = pVpm->adjusted_crushing_pressure_n2[i] / UNITS_FACTOR * 101325.0f;+ −
/*+ −
initial_allowable_grad_he_pa =+ −
pVpm->initial_allowable_gradient_he[i] / UNITS_FACTOR * 101325.0f;+ −
initial_allowable_grad_n2_pa =+ −
pVpm->initial_allowable_gradient_n2[i] / UNITS_FACTOR * 101325.0f;+ −
*/+ −
if (pVpm->max_actual_gradient[i] > pVpm->initial_allowable_gradient_n2[i])+ −
{+ −
new_critical_radius_n2 =+ −
SURFACE_TENSION_GAMMA * 2.0f *+ −
(SKIN_COMPRESSION_GAMMAC - SURFACE_TENSION_GAMMA) /+ −
(max_actual_gradient_pascals * SKIN_COMPRESSION_GAMMAC -+ −
SURFACE_TENSION_GAMMA * adj_crush_pressure_n2_pascals);+ −
+ −
pVpm->adjusted_critical_radius_n2[i] =+ −
initial_critical_radius_n2[i]+ −
+ (initial_critical_radius_n2[i] - new_critical_radius_n2)+ −
* exp(-(*surface_interval_time) / REGENERATION_TIME_CONSTANT);+ −
+ −
} else {+ −
pVpm->adjusted_critical_radius_n2[i] =+ −
initial_critical_radius_n2[i];+ −
}+ −
if (pVpm->max_actual_gradient[i] > pVpm->initial_allowable_gradient_he[i])+ −
{+ −
new_critical_radius_he =+ −
SURFACE_TENSION_GAMMA * 2.0f *+ −
(SKIN_COMPRESSION_GAMMAC - SURFACE_TENSION_GAMMA) /+ −
(max_actual_gradient_pascals * SKIN_COMPRESSION_GAMMAC -+ −
SURFACE_TENSION_GAMMA * adj_crush_pressure_he_pascals);+ −
+ −
pVpm->adjusted_critical_radius_he[i] =+ −
initial_critical_radius_he[i]+ −
+ ( initial_critical_radius_he[i] - new_critical_radius_he)+ −
* exp(-(*surface_interval_time) / REGENERATION_TIME_CONSTANT);+ −
} else {+ −
pVpm->adjusted_critical_radius_he[i] =+ −
initial_critical_radius_he[i];+ −
}+ −
}+ −
return 0;+ −
} /* vpm_repetitive_algorithm */+ −