Code:
/* =============================================================================== */
/* 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(real *starting_ambient_pressure,
real *ending_ambient_pressure,
real *rate,
integer *i)
{
/* Local variables */
static real time, last_diff_change, mid_range_nitrogen_pressure;
static integer j;
static real 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;
static real mid_range_ambient_pressure,
high_bound_helium_pressure,
helium_rate,
differential_change;
/* 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;
starting_gas_tension =
initial_helium_pressure[*i - 1] +
initial_nitrogen_pressure[*i - 1] +
constant_pressure_other_gases;
function_at_low_bound =
*starting_ambient_pressure -
starting_gas_tension -
gradient_onset_of_imperm;
high_bound_helium_pressure =
schreiner_equation__(&initial_inspired_he_pressure,
&helium_rate,
&high_bound,
&helium_time_constant[*i - 1],
&initial_helium_pressure[*i - 1]);
high_bound_nitrogen_pressure =
schreiner_equation__(&initial_inspired_n2_pressure,
&nitrogen_rate,
&high_bound,
&nitrogen_time_constant[*i - 1],
&initial_nitrogen_pressure[*i - 1]);
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.) {
printf("\nERROR! ROOT IS NOT WITHIN BRACKETS");
pause();
}
/* =============================================================================== */
/* 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.) {
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 * .5;
mid_range_time = time + differential_change;
mid_range_ambient_pressure = *starting_ambient_pressure + *rate * mid_range_time;
mid_range_helium_pressure =
schreiner_equation__(&initial_inspired_he_pressure,
&helium_rate,
&mid_range_time,
&helium_time_constant[*i - 1],
&initial_helium_pressure[*i - 1]);
mid_range_nitrogen_pressure =
schreiner_equation__(&initial_inspired_n2_pressure,
&nitrogen_rate,
&mid_range_time,
&nitrogen_time_constant[*i - 1],
&initial_nitrogen_pressure[*i - 1]);
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.) {
time = mid_range_time;
}
if (fabs(differential_change) < .001 ||
function_at_mid_range == 0.) {
goto L100;
}
}
printf("\nERROR! ROOT SEARCH EXCEEDED MAXIMUM ITERATIONS");
pause();
/* =============================================================================== */
/* 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 - 1] = mid_range_ambient_pressure;
gas_tension_onset_of_imperm[*i - 1] = gas_tension_at_mid_range;
return 0;
} /* onset_of_impermeability */
