Mercurial > public > ostc4
view Common/Src/calc_crush.c @ 155:4fd8bbc7d841 FlipDisplay
Do hard sync after communication timeout > 1 second
author | Ideenmodellierer |
---|---|
date | Sun, 03 Mar 2019 17:46:50 +0100 |
parents | 8f8ea3a32e82 |
children |
line wrap: on
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 */