view Common/Src/calc_crush.c @ 488:9eeab3fead8f

Added "I2C_DeInit();" in hardware detection routines. It's the recommended way to do and solves potential issues with older OSTC4 hardware generations. Increased RTE and RTE required numbers updated "Current build"
author heinrichsweikamp
date Sat, 27 Jun 2020 18:40:32 +0200
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 */