mifr/subsurface
1
0
Fork 0
mirror of https://github.com/subsurface/subsurface.git synced 2025-02-19 22:16:15 +00:00
Code Issues Projects Releases Packages Wiki Activity

VPM-B: Add simple Boyle's law compensation.

Browse source 
This is a very basic implementation that uses bin search for solving the cubic.
It's not called anywhere at this stage, needs to be changed to analytic solution.

Signed-off-by: Jan Darowski <jan.darowski@gmail.com>
This commit is contained in:
Jan Darowski 2015-08-15 14:03:08 +02:00
parent a3e87bf1a1
commit a828d528f9
1 changed files with 48 additions and 9 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

57
deco.c
View file

@ -114,6 +114,11 @@ double allowable_n2_gradient[16];
double allowable_he_gradient[16];
double total_gradient[16];
double bottom_n2_gradient[16];
double bottom_he_gradient[16];
double initial_n2_gradient[16];
double initial_he_gradient[16];
static double tissue_tolerance_calc(const struct dive *dive)
{
@ -231,8 +236,8 @@ void vpmb_start_gradient()
double gradient_n2, gradient_he;
for (ci = 0; ci < 16; ++ci) {
allowable_n2_gradient[ci] = 2.0 * (vpmb_config.surface_tension_gamma / vpmb_config.skin_compression_gammaC) * ((vpmb_config.skin_compression_gammaC - vpmb_config.surface_tension_gamma) / n2_regen_radius[ci]);
allowable_he_gradient[ci] = 2.0 * (vpmb_config.surface_tension_gamma / vpmb_config.skin_compression_gammaC) * ((vpmb_config.skin_compression_gammaC - vpmb_config.surface_tension_gamma) / he_regen_radius[ci]);
initial_n2_gradient[ci] = bottom_n2_gradient[ci] = allowable_n2_gradient[ci] = 2.0 * (vpmb_config.surface_tension_gamma / vpmb_config.skin_compression_gammaC) * ((vpmb_config.skin_compression_gammaC - vpmb_config.surface_tension_gamma) / n2_regen_radius[ci]);
initial_he_gradient[ci] = bottom_he_gradient[ci] = allowable_he_gradient[ci] = 2.0 * (vpmb_config.surface_tension_gamma / vpmb_config.skin_compression_gammaC) * ((vpmb_config.skin_compression_gammaC - vpmb_config.surface_tension_gamma) / he_regen_radius[ci]);
total_gradient[ci] = ((allowable_n2_gradient[ci] * tissue_n2_sat[ci]) + (allowable_he_gradient[ci] * tissue_he_sat[ci])) / (tissue_n2_sat[ci] + tissue_he_sat[ci]);
}
@ -244,19 +249,53 @@ void vpmb_next_gradient(double deco_time)
double gradient_n2, gradient_he;
double n2_b, n2_c;
double he_b, he_c;
deco_time /= 60.0 ;
double desat_time = deco_time / 60.0;
for (ci = 0; ci < 16; ++ci) {
n2_b = allowable_n2_gradient[ci] + ((vpmb_config.crit_volume_lambda * vpmb_config.surface_tension_gamma) / (vpmb_config.skin_compression_gammaC * (deco_time + buehlmann_N2_t_halflife[ci] * 60.0 / log(2.0))));
he_b = allowable_he_gradient[ci] + ((vpmb_config.crit_volume_lambda * vpmb_config.surface_tension_gamma) / (vpmb_config.skin_compression_gammaC * (deco_time + buehlmann_He_t_halflife[ci] * 60.0 / log(2.0))));
n2_b = initial_n2_gradient[ci] + (vpmb_config.crit_volume_lambda * vpmb_config.surface_tension_gamma) / (vpmb_config.skin_compression_gammaC * desat_time);
he_b = initial_he_gradient[ci] + (vpmb_config.crit_volume_lambda * vpmb_config.surface_tension_gamma) / (vpmb_config.skin_compression_gammaC * desat_time);
n2_c = vpmb_config.surface_tension_gamma * vpmb_config.surface_tension_gamma * vpmb_config.crit_volume_lambda * max_n2_crushing_pressure[ci];
n2_c = n2_c / (vpmb_config.skin_compression_gammaC * vpmb_config.skin_compression_gammaC * (deco_time + buehlmann_N2_t_halflife[ci] * 60.0 / log(2.0)));
n2_c = n2_c / (vpmb_config.skin_compression_gammaC * vpmb_config.skin_compression_gammaC * desat_time);
he_c = vpmb_config.surface_tension_gamma * vpmb_config.surface_tension_gamma * vpmb_config.crit_volume_lambda * max_he_crushing_pressure[ci];
he_c = he_c / (vpmb_config.skin_compression_gammaC * vpmb_config.skin_compression_gammaC * (deco_time + buehlmann_He_t_halflife[ci] * 60.0 / log(2.0)));
he_c = he_c / (vpmb_config.skin_compression_gammaC * vpmb_config.skin_compression_gammaC * desat_time);
allowable_n2_gradient[ci] = 0.5 * ( n2_b + sqrt(n2_b * n2_b - 4.0 * n2_c));
allowable_he_gradient[ci] = 0.5 * ( he_b + sqrt(he_b * he_b - 4.0 * he_c));
bottom_n2_gradient[ci] = allowable_n2_gradient[ci] = 0.5 * ( n2_b + sqrt(n2_b * n2_b - 4.0 * n2_c));
bottom_he_gradient[ci] = allowable_he_gradient[ci] = 0.5 * ( he_b + sqrt(he_b * he_b - 4.0 * he_c));
total_gradient[ci] = ((allowable_n2_gradient[ci] * tissue_n2_sat[ci]) + (allowable_he_gradient[ci] * tissue_he_sat[ci])) / (tissue_n2_sat[ci] + tissue_he_sat[ci]);
}
}
double update_gradient(double first_stop_pressure, double next_stop_pressure, double first_gradient)
{
double first_radius = 2.0 * vpmb_config.surface_tension_gamma / first_gradient;
double A = next_stop_pressure;
double B = -2.0 * vpmb_config.surface_tension_gamma;
double C = (first_stop_pressure + 2.0 * vpmb_config.surface_tension_gamma / first_radius) * pow(first_radius, 3.0);
double low = first_radius;
double high = first_radius * pow(first_stop_pressure / next_stop_pressure, (1.0/3.0));
double next_radius;
double value;
int ci;
for (ci = 0; ci < 100; ++ci){
next_radius = (high + low) /2.0;
value = A * pow(next_radius, 3.0) - B * next_radius * next_radius - C;
if (value < 0)
low = next_radius;
else
high = next_radius;
}
return 2.0 * vpmb_config.surface_tension_gamma / next_radius;
}
void boyles_law(double first_stop_pressure, double next_stop_pressure)
{
int ci;
for (ci = 0; ci < 16; ++ci) {
allowable_n2_gradient[ci] = update_gradient(first_stop_pressure, next_stop_pressure, bottom_n2_gradient[ci]);
allowable_he_gradient[ci] = update_gradient(first_stop_pressure, next_stop_pressure, bottom_he_gradient[ci]);
total_gradient[ci] = ((allowable_n2_gradient[ci] * tissue_n2_sat[ci]) + (allowable_he_gradient[ci] * tissue_he_sat[ci])) / (tissue_n2_sat[ci] + tissue_he_sat[ci]);
}