Merge branch 'new-vpm' of https://github.com/Slagvi/subsurface

10 files changed, 387 insertions, 137 deletions

diff --git a/deco.c b/deco.c
index 9d0ec9d53..f1bc8039f 100644
--- a/deco.c
+++ b/deco.c
@@ -32,6 +32,19 @@ struct buehlmann_config {
 };
 struct buehlmann_config buehlmann_config = { 1.0, 1.01, 0, 0.75, 0.35, 1.0, false };
 
+//! Option structure for VPM-B decompression.
+struct vpmb_config {
+	double crit_radius_N2;            //! Critical radius of N2 nucleon (microns).
+	double crit_radius_He;            //! Critical radius of He nucleon (microns).
+	double crit_volume_lambda;        //! Constant corresponding to critical gas volume.
+	double gradient_of_imperm;        //! Gradient after which bubbles become impermeable.
+	double surface_tension_gamma;     //! Nucleons surface tension constant.
+	double skin_compression_gammaC;   //!
+	double regeneration_time;         //! Time needed for the bubble to regenerate to the start radius.
+	double other_gases_pressure;      //! Always present pressure of other gasses in tissues.
+};
+struct vpmb_config vpmb_config = { 0.6, 0.5, 250.0, 8.2, 0.179, 2.57, 20160, 0.1359888 };
+
 const double buehlmann_N2_a[] = { 1.1696, 1.0, 0.8618, 0.7562,
 				  0.62, 0.5043, 0.441, 0.4,
 				  0.375, 0.35, 0.3295, 0.3065,
@@ -89,6 +102,19 @@ double tolerated_by_tissue[16];
 double tissue_inertgas_saturation[16];
 double buehlmann_inertgas_a[16], buehlmann_inertgas_b[16];
 
+double max_n2_crushing_pressure[16];
+double max_he_crushing_pressure[16];
+
+double crushing_onset_tension[16];            // total inert gas tension in the t* moment
+double n2_regen_radius[16];                   // rs
+double he_regen_radius[16];
+double max_ambient_pressure;                  // last moment we were descending
+
+double allowable_n2_gradient[16];
+double allowable_he_gradient[16];
+double total_gradient[16];
+
+
 static double tissue_tolerance_calc(const struct dive *dive)
 {
 	int ci = -1;
@@ -184,6 +210,135 @@ double he_factor(int period_in_seconds, int ci)
 	return cache[ci].last_factor;
 }
 
+bool is_vpmb_ok(double pressure)
+{
+	int ci;
+	double gradient;
+	double gas_tension;
+
+	for (ci = 0; ci < 16; ++ci) {
+		gas_tension = tissue_n2_sat[ci] + tissue_he_sat[ci] + vpmb_config.other_gases_pressure;
+		gradient = gas_tension - pressure;
+		if (gradient > total_gradient[ci])
+			return false;
+	}
+	return true;
+}
+
+void vpmb_start_gradient()
+{
+	int ci;
+	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]);
+
+		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]);
+	}
+}
+
+void vpmb_next_gradient(double deco_time)
+{
+	int ci;
+	double gradient_n2, gradient_he;
+	double n2_b, n2_c;
+	double he_b, he_c;
+	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_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)));
+		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)));
+
+		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));
+
+		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]);
+	}
+}
+
+void nuclear_regeneration(double time)
+{
+	time /= 60.0;
+	int ci;
+	double crushing_radius_N2, crushing_radius_He;
+	for (ci = 0; ci < 16; ++ci) {
+		//rm
+		crushing_radius_N2 = 1.0 / (max_n2_crushing_pressure[ci] / (2.0 * (vpmb_config.skin_compression_gammaC - vpmb_config.surface_tension_gamma)) + 1.0 / vpmb_config.crit_radius_N2);
+		crushing_radius_He = 1.0 / (max_he_crushing_pressure[ci] / (2.0 * (vpmb_config.skin_compression_gammaC - vpmb_config.surface_tension_gamma)) + 1.0 / vpmb_config.crit_radius_He);
+		//rs
+		n2_regen_radius[ci] = crushing_radius_N2 + (vpmb_config.crit_radius_N2 - crushing_radius_N2) * (1.0 - exp (-time / vpmb_config.regeneration_time));
+		he_regen_radius[ci] = crushing_radius_He + (vpmb_config.crit_radius_He - crushing_radius_He) * (1.0 - exp (-time / vpmb_config.regeneration_time));
+	}
+}
+
+// Calculates the nucleons inner pressure during the impermeable period
+double calc_inner_pressure(double crit_radius, double onset_tension, double current_ambient_pressure)
+{
+	double onset_radius = 1.0 / (vpmb_config.gradient_of_imperm / (2.0 * (vpmb_config.skin_compression_gammaC - vpmb_config.surface_tension_gamma)) + 1.0 / crit_radius);
+
+	// A*r^3 + B*r^2 + C == 0
+	// Solved with the help of mathematica
+
+	double A = current_ambient_pressure - vpmb_config.gradient_of_imperm + (2.0 * (vpmb_config.skin_compression_gammaC - vpmb_config.surface_tension_gamma)) / onset_radius;
+	double B = 2.0 * (vpmb_config.skin_compression_gammaC - vpmb_config.surface_tension_gamma);
+	double C = onset_tension * pow(onset_radius, 3);
+
+	double BA = B/A;
+	double CA = C/A;
+
+	double discriminant = CA * (4 * BA * BA * BA + 27 * CA);
+
+	// Let's make sure we have a real solution:
+	if (discriminant < 0.0) {
+		// This should better not happen
+		report_error("Complex solution for inner pressure encountered!\n A=%f\tB=%f\tC=%f\n", A, B, C);
+		return 0.0;
+	}
+	double denominator = pow(BA * BA * BA + 1.5 * (9 * CA + sqrt(3.0) * sqrt(discriminant)), 1/3.0);
+	double current_radius = (BA + BA * BA / denominator + denominator) / 3.0;
+
+	return onset_tension * onset_radius * onset_radius * onset_radius / (current_radius * current_radius * current_radius);
+}
+
+// Calculates the crushing pressure in the given moment. Updates crushing_onset_tension and critical radius if needed
+void calc_crushing_pressure(double pressure)
+{
+	int ci;
+	double gradient;
+	double gas_tension;
+	double n2_crushing_pressure, he_crushing_pressure;
+	double n2_inner_pressure, he_inner_pressure;
+
+	for (ci = 0; ci < 16; ++ci) {
+		gas_tension = tissue_n2_sat[ci] + tissue_he_sat[ci] + vpmb_config.other_gases_pressure;
+		gradient = pressure - gas_tension;
+
+		if (gradient <= vpmb_config.gradient_of_imperm) {	// permeable situation
+			n2_crushing_pressure = he_crushing_pressure = gradient;
+			crushing_onset_tension[ci] = gas_tension;
+		}
+		else {	// impermeable
+			if (max_ambient_pressure >= pressure)
+				return;
+
+			n2_inner_pressure = calc_inner_pressure(vpmb_config.crit_radius_N2, crushing_onset_tension[ci], pressure);
+			he_inner_pressure = calc_inner_pressure(vpmb_config.crit_radius_He, crushing_onset_tension[ci], pressure);
+
+			n2_crushing_pressure = pressure - n2_inner_pressure;
+			he_crushing_pressure = pressure - he_inner_pressure;
+		}
+		max_n2_crushing_pressure[ci] = MAX(max_n2_crushing_pressure[ci], n2_crushing_pressure);
+		max_he_crushing_pressure[ci] = MAX(max_he_crushing_pressure[ci], he_crushing_pressure);
+	}
+	max_ambient_pressure = MAX(pressure, max_ambient_pressure);
+}
+
 /* add period_in_seconds at the given pressure and gas to the deco calculation */
 double add_segment(double pressure, const struct gasmix *gasmix, int period_in_seconds, int ccpo2, const struct dive *dive, int sac)
 {
@@ -206,6 +361,7 @@ double add_segment(double pressure, const struct gasmix *gasmix, int period_in_s
 		tissue_n2_sat[ci] += n2_satmult * pn2_oversat * n2_f;
 		tissue_he_sat[ci] += he_satmult * phe_oversat * he_f;
 	}
+	calc_crushing_pressure(pressure);
 	return tissue_tolerance_calc(dive);
 }
 
@@ -229,10 +385,13 @@ void clear_deco(double surface_pressure)
 	for (ci = 0; ci < 16; ci++) {
 		tissue_n2_sat[ci] = (surface_pressure - WV_PRESSURE) * N2_IN_AIR / 1000;
 		tissue_he_sat[ci] = 0.0;
+		max_n2_crushing_pressure[ci] = 0.0;
+		max_he_crushing_pressure[ci] = 0.0;
 	}
 	gf_low_pressure_this_dive = surface_pressure;
 	if (!buehlmann_config.gf_low_at_maxdepth)
 		gf_low_pressure_this_dive += buehlmann_config.gf_low_position_min;
+	max_ambient_pressure = 0.0;
 }
 
 void cache_deco_state(double tissue_tolerance, char **cached_datap)
diff --git a/dive.h b/dive.h
index f7044fba2..042878c1b 100644
--- a/dive.h
+++ b/dive.h
@@ -783,6 +783,10 @@ extern unsigned int deco_allowed_depth(double tissues_tolerance, double surface_
 extern void set_gf(short gflow, short gfhigh, bool gf_low_at_maxdepth);
 extern void cache_deco_state(double, char **datap);
 extern double restore_deco_state(char *data);
+extern void nuclear_regeneration(double time);
+extern void vpmb_start_gradient();
+extern void vpmb_next_gradient(double deco_time);
+extern bool is_vpmb_ok(double pressure);
 
 /* this should be converted to use our types */
 struct divedatapoint {
diff --git a/planner.c b/planner.c
index 35158fe12..83e475758 100644
--- a/planner.c
+++ b/planner.c
@@ -836,11 +836,15 @@ bool trial_ascent(int trial_depth, int stoplevel, int avg_depth, int bottom_time
 		tissue_tolerance = add_segment(depth_to_mbar(trial_depth, &displayed_dive) / 1000.0,
 					       gasmix,
 					       TIMESTEP, po2, &displayed_dive, prefs.decosac);
-		if (deco_allowed_depth(tissue_tolerance, surface_pressure, &displayed_dive, 1) > trial_depth - deltad) {
+		if (prefs.deco_mode != VPMB && deco_allowed_depth(tissue_tolerance, surface_pressure, &displayed_dive, 1) > trial_depth - deltad) {
 			/* We should have stopped */
 			clear_to_ascend = false;
 			break;
 		}
+		if (prefs.deco_mode == VPMB && (!is_vpmb_ok(depth_to_mbar(trial_depth, &displayed_dive) / 1000.0))){
+			clear_to_ascend = false;
+			break;
+		}
 		trial_depth -= deltad;
 	}
 	restore_deco_state(trial_cache);
@@ -865,6 +869,14 @@ bool enough_gas(int current_cylinder)
 
 bool plan(struct diveplan *diveplan, char **cached_datap, bool is_planner, bool show_disclaimer)
 {
+	int bottom_depth;
+	int bottom_gi;
+	int bottom_stopidx;
+
+	bool is_final_plan = true;
+	int deco_time;
+	int previous_deco_time;
+	char *bottom_cache = NULL;
 	struct sample *sample;
 	int po2;
 	int transitiontime, gi;
@@ -882,7 +894,7 @@ bool plan(struct diveplan *diveplan, char **cached_datap, bool is_planner, bool
 	int avg_depth, max_depth, bottom_time = 0;
 	int last_ascend_rate;
 	int best_first_ascend_cylinder;
-	struct gasmix gas;
+	struct gasmix gas, bottom_gas;
 	int o2time = 0;
 	int breaktime = -1;
 	int breakcylinder = 0;
@@ -955,7 +967,8 @@ bool plan(struct diveplan *diveplan, char **cached_datap, bool is_planner, bool
 	/* Keep time during the ascend */
 	bottom_time = clock = previous_point_time = displayed_dive.dc.sample[displayed_dive.dc.samples - 1].time.seconds;
 	gi = gaschangenr - 1;
-	if(prefs.recreational_mode) {
+
+	if(prefs.deco_mode == RECREATIONAL) {
 		bool safety_stop = prefs.safetystop && max_depth >= 10000;
 		track_ascent_gas(depth, &displayed_dive.cylinder[current_cylinder], avg_depth, bottom_time, safety_stop);
 		// How long can we stay at the current depth and still directly ascent to the surface?
@@ -1003,7 +1016,6 @@ bool plan(struct diveplan *diveplan, char **cached_datap, bool is_planner, bool
 
 		free(stoplevels);
 		free(gaschanges);
-
 		return(false);
 	}
 
@@ -1018,147 +1030,185 @@ bool plan(struct diveplan *diveplan, char **cached_datap, bool is_planner, bool
 		       (get_o2(&gas) + 5) / 10, (get_he(&gas) + 5) / 10, gaschanges[best_first_ascend_cylinder].depth / 1000.0);
 #endif
 	}
-	while (1) {
-		/* We will break out when we hit the surface */
-		do {
-			/* Ascend to next stop depth */
-			int deltad = ascent_velocity(depth, avg_depth, bottom_time) * TIMESTEP;
-			if (ascent_velocity(depth, avg_depth, bottom_time) != last_ascend_rate) {
-				plan_add_segment(diveplan, clock - previous_point_time, depth, gas, po2, false);
+
+	// VPM-B or Buehlmann Deco
+	nuclear_regeneration(clock);
+	vpmb_start_gradient();
+	previous_deco_time = 100000000;
+	deco_time = 10000000;
+	cache_deco_state(tissue_tolerance, &bottom_cache);  // Lets us make several iterations
+	bottom_depth = depth;
+	bottom_gi = gi;
+	bottom_gas = gas;
+	bottom_stopidx = stopidx;
+	//CVA
+	do {
+		is_final_plan = (prefs.deco_mode == BUEHLMANN) || (previous_deco_time - deco_time < 10);  // CVA time converges
+		previous_deco_time = deco_time;
+		restore_deco_state(bottom_cache);
+
+		depth = bottom_depth;
+		gi = bottom_gi;
+		clock = previous_point_time = bottom_time;
+		gas = bottom_gas;
+		stopping = false;
+		decodive = false;
+		stopidx = bottom_stopidx;
+		breaktime = -1;
+		breakcylinder = 0;
+		o2time = 0;
+		last_ascend_rate = ascent_velocity(depth, avg_depth, bottom_time);
+		if ((current_cylinder = get_gasidx(&displayed_dive, &gas)) == -1) {
+			report_error(translate("gettextFromC", "Can't find gas %s"), gasname(&gas));
+			current_cylinder = 0;
+		}
+		vpmb_next_gradient(deco_time);
+
+		while (1) {
+			/* We will break out when we hit the surface */
+			do {
+				/* Ascend to next stop depth */
+				int deltad = ascent_velocity(depth, avg_depth, bottom_time) * TIMESTEP;
+				if (ascent_velocity(depth, avg_depth, bottom_time) != last_ascend_rate) {
+					if (is_final_plan)
+						plan_add_segment(diveplan, clock - previous_point_time, depth, gas, po2, false);
+					previous_point_time = clock;
+					stopping = false;
+					last_ascend_rate = ascent_velocity(depth, avg_depth, bottom_time);
+				}
+				if (depth - deltad < stoplevels[stopidx])
+					deltad = depth - stoplevels[stopidx];
+
+				tissue_tolerance = add_segment(depth_to_mbar(depth, &displayed_dive) / 1000.0,
+								&displayed_dive.cylinder[current_cylinder].gasmix,
+								TIMESTEP, po2, &displayed_dive, prefs.decosac);
+				clock += TIMESTEP;
+				depth -= deltad;
+			} while (depth > 0 && depth > stoplevels[stopidx]);
+
+			if (depth <= 0)
+				break; /* We are at the surface */
+
+			if (gi >= 0 && stoplevels[stopidx] <= gaschanges[gi].depth) {
+				/* We have reached a gas change.
+				 * Record this in the dive plan */
+				if (is_final_plan)
+					plan_add_segment(diveplan, clock - previous_point_time, depth, gas, po2, false);
 				previous_point_time = clock;
-				stopping = false;
-				last_ascend_rate = ascent_velocity(depth, avg_depth, bottom_time);
-			}
-			if (depth - deltad < stoplevels[stopidx])
-				deltad = depth - stoplevels[stopidx];
+				stopping = true;
 
-			tissue_tolerance = add_segment(depth_to_mbar(depth, &displayed_dive) / 1000.0,
-						       &displayed_dive.cylinder[current_cylinder].gasmix,
-						       TIMESTEP, po2, &displayed_dive, prefs.decosac);
-			clock += TIMESTEP;
-			depth -= deltad;
-		} while (depth > 0 && depth > stoplevels[stopidx]);
-
-		if (depth <= 0)
-			break; /* We are at the surface */
-
-		if (gi >= 0 && stoplevels[stopidx] <= gaschanges[gi].depth) {
-			/* We have reached a gas change.
-			 * Record this in the dive plan */
-			plan_add_segment(diveplan, clock - previous_point_time, depth, gas, po2, false);
-			previous_point_time = clock;
-			stopping = true;
-
-			/* Check we need to change cylinder.
-			 * We might not if the cylinder was chosen by the user
-			 * or user has selected only to switch only at required stops.
-			 * If current gas is hypoxic, we want to switch asap */
-			if (current_cylinder != gaschanges[gi].gasidx) {
-				if (!prefs.switch_at_req_stop ||
-						!trial_ascent(depth, stoplevels[stopidx - 1], avg_depth, bottom_time, tissue_tolerance,
-						&displayed_dive.cylinder[current_cylinder].gasmix, po2, diveplan->surface_pressure / 1000.0) || get_o2(&displayed_dive.cylinder[current_cylinder].gasmix) < 160) {
-					current_cylinder = gaschanges[gi].gasidx;
+				/* Check we need to change cylinder.
+				 * We might not if the cylinder was chosen by the user
+				 * or user has selected only to switch only at required stops.
+				 * If current gas is hypoxic, we want to switch asap */
+				if (current_cylinder != gaschanges[gi].gasidx) {
+					if (!prefs.switch_at_req_stop ||
+							!trial_ascent(depth, stoplevels[stopidx - 1], avg_depth, bottom_time, tissue_tolerance,
+							&displayed_dive.cylinder[current_cylinder].gasmix, po2, diveplan->surface_pressure / 1000.0) || get_o2(&displayed_dive.cylinder[current_cylinder].gasmix) < 160) {
+						current_cylinder = gaschanges[gi].gasidx;
+						gas = displayed_dive.cylinder[current_cylinder].gasmix;
+#if DEBUG_PLAN & 16
+						printf("switch to gas %d (%d/%d) @ %5.2lfm\n", gaschanges[gi].gasidx,
+							(get_o2(&gas) + 5) / 10, (get_he(&gas) + 5) / 10, gaschanges[gi].depth / 1000.0);
+#endif
+						/* Stop for the minimum duration to switch gas */
+						tissue_tolerance = add_segment(depth_to_mbar(depth, &displayed_dive) / 1000.0,
+							&displayed_dive.cylinder[current_cylinder].gasmix,
+							prefs.min_switch_duration, po2, &displayed_dive, prefs.decosac);
+						clock += prefs.min_switch_duration;
+					} else {
+						pendinggaschange = true;
+					}
+				gi--;
+				}
+			}
+			--stopidx;
+
+			/* Save the current state and try to ascend to the next stopdepth */
+			while (1) {
+				/* Check if ascending to next stop is clear, go back and wait if we hit the ceiling on the way */
+				if (trial_ascent(depth, stoplevels[stopidx], avg_depth, bottom_time, tissue_tolerance,
+						&displayed_dive.cylinder[current_cylinder].gasmix, po2, diveplan->surface_pressure / 1000.0))
+					break; /* We did not hit the ceiling */
+
+				/* Add a minute of deco time and then try again */
+				decodive = true;
+				if (!stopping) {
+					/* The last segment was an ascend segment.
+					 * Add a waypoint for start of this deco stop */
+					if (is_final_plan)
+						plan_add_segment(diveplan, clock - previous_point_time, depth, gas, po2, false);
+					previous_point_time = clock;
+					stopping = true;
+				}
+				if (pendinggaschange) {
+					current_cylinder = gaschanges[gi + 1].gasidx;
 					gas = displayed_dive.cylinder[current_cylinder].gasmix;
 #if DEBUG_PLAN & 16
-					printf("switch to gas %d (%d/%d) @ %5.2lfm\n", gaschanges[gi].gasidx,
-						(get_o2(&gas) + 5) / 10, (get_he(&gas) + 5) / 10, gaschanges[gi].depth / 1000.0);
+					printf("switch to gas %d (%d/%d) @ %5.2lfm\n", gaschanges[gi + 1].gasidx,
+						(get_o2(&gas) + 5) / 10, (get_he(&gas) + 5) / 10, gaschanges[gi + 1].depth / 1000.0);
 #endif
 					/* Stop for the minimum duration to switch gas */
 					tissue_tolerance = add_segment(depth_to_mbar(depth, &displayed_dive) / 1000.0,
 						&displayed_dive.cylinder[current_cylinder].gasmix,
 						prefs.min_switch_duration, po2, &displayed_dive, prefs.decosac);
 					clock += prefs.min_switch_duration;
-				} else {
-					pendinggaschange = true;
+					pendinggaschange = false;
 				}
-			gi--;
-			}
-		}
 
-		--stopidx;
+				/* Deco stop should end when runtime is at a whole minute */
+				int this_decotimestep;
+				this_decotimestep = DECOTIMESTEP - clock % DECOTIMESTEP;
 
-		/* Save the current state and try to ascend to the next stopdepth */
-		while (1) {
-			/* Check if ascending to next stop is clear, go back and wait if we hit the ceiling on the way */
-			if (trial_ascent(depth, stoplevels[stopidx], avg_depth, bottom_time, tissue_tolerance,
-					 &displayed_dive.cylinder[current_cylinder].gasmix, po2, diveplan->surface_pressure / 1000.0))
-				break; /* We did not hit the ceiling */
-
-			/* Add a minute of deco time and then try again */
-			decodive = true;
-			if (!stopping) {
-				/* The last segment was an ascend segment.
-				 * Add a waypoint for start of this deco stop */
-				plan_add_segment(diveplan, clock - previous_point_time, depth, gas, po2, false);
-				previous_point_time = clock;
-				stopping = true;
-			}
-
-			if (pendinggaschange) {
-				current_cylinder = gaschanges[gi + 1].gasidx;
-				gas = displayed_dive.cylinder[current_cylinder].gasmix;
-#if DEBUG_PLAN & 16
-				printf("switch to gas %d (%d/%d) @ %5.2lfm\n", gaschanges[gi + 1].gasidx,
-					(get_o2(&gas) + 5) / 10, (get_he(&gas) + 5) / 10, gaschanges[gi + 1].depth / 1000.0);
-#endif
-				/* Stop for the minimum duration to switch gas */
 				tissue_tolerance = add_segment(depth_to_mbar(depth, &displayed_dive) / 1000.0,
-					&displayed_dive.cylinder[current_cylinder].gasmix,
-					prefs.min_switch_duration, po2, &displayed_dive, prefs.decosac);
-				clock += prefs.min_switch_duration;
-				pendinggaschange = false;
-			}
-
-			/* Deco stop should end when runtime is at a whole minute */
-			int this_decotimestep;
-			this_decotimestep = DECOTIMESTEP - clock % DECOTIMESTEP;
-
-			tissue_tolerance = add_segment(depth_to_mbar(depth, &displayed_dive) / 1000.0,
-						       &displayed_dive.cylinder[current_cylinder].gasmix,
-						       this_decotimestep, po2, &displayed_dive, prefs.decosac);
-			clock += this_decotimestep;
-			/* Finish infinite deco */
-			if(clock >= 48 * 3600 && depth >= 6000) {
-				error = LONGDECO;
-				break;
-			}
-			if (prefs.doo2breaks) {
-				if (get_o2(&displayed_dive.cylinder[current_cylinder].gasmix) == 1000) {
-					o2time += DECOTIMESTEP;
-					if (o2time >= 12 * 60) {
-						breaktime = 0;
-						breakcylinder = current_cylinder;
-						plan_add_segment(diveplan, clock - previous_point_time, depth, gas, po2, false);
-						previous_point_time = clock;
-						current_cylinder = 0;
-						gas = displayed_dive.cylinder[current_cylinder].gasmix;
-					}
-				} else {
-					if (breaktime >= 0) {
-						breaktime += DECOTIMESTEP;
-						if (breaktime >= 6 * 60) {
-							o2time = 0;
-							plan_add_segment(diveplan, clock - previous_point_time, depth, gas, po2, false);
+								&displayed_dive.cylinder[current_cylinder].gasmix,
+								this_decotimestep, po2, &displayed_dive, prefs.decosac);
+				clock += this_decotimestep;
+				/* Finish infinite deco */
+				if(clock >= 48 * 3600 && depth >= 6000) {
+					error = LONGDECO;
+					break;
+				}
+				if (prefs.doo2breaks) {
+					if (get_o2(&displayed_dive.cylinder[current_cylinder].gasmix) == 1000) {
+						o2time += DECOTIMESTEP;
+						if (o2time >= 12 * 60) {
+							breaktime = 0;
+							breakcylinder = current_cylinder;
+							if (is_final_plan)
+								plan_add_segment(diveplan, clock - previous_point_time, depth, gas, po2, false);
 							previous_point_time = clock;
-							current_cylinder = breakcylinder;
+							current_cylinder = 0;
 							gas = displayed_dive.cylinder[current_cylinder].gasmix;
-							breaktime = -1;
+						}
+					} else {
+						if (breaktime >= 0) {
+							breaktime += DECOTIMESTEP;
+							if (breaktime >= 6 * 60) {
+								o2time = 0;
+								if (is_final_plan)
+									plan_add_segment(diveplan, clock - previous_point_time, depth, gas, po2, false);
+								previous_point_time = clock;
+								current_cylinder = breakcylinder;
+								gas = displayed_dive.cylinder[current_cylinder].gasmix;
+								breaktime = -1;
+							}
 						}
 					}
 				}
 			}
+			if (stopping) {
+				/* Next we will ascend again. Add a waypoint if we have spend deco time */
+				if (is_final_plan)
+					plan_add_segment(diveplan, clock - previous_point_time, depth, gas, po2, false);
+				previous_point_time = clock;
+				stopping = false;
+			}
 		}
-		if (stopping) {
-			/* Next we will ascend again. Add a waypoint if we have spend deco time */
-			plan_add_segment(diveplan, clock - previous_point_time, depth, gas, po2, false);
-			previous_point_time = clock;
-			stopping = false;
-		}
-	}
 
-	/* We made it to the surface
-	 * Create the final dive, add the plan to the notes and fixup some internal
-	 * data that we need to be there when plotting the dive */
+		deco_time = clock - bottom_time;
+	} while (!is_final_plan);
+
 	plan_add_segment(diveplan, clock - previous_point_time, 0, gas, po2, false);
 	create_dive_from_plan(diveplan, is_planner);
 	add_plan_to_notes(diveplan, &displayed_dive, show_disclaimer, error);
diff --git a/pref.h b/pref.h
index 0470d0e3b..fcb051bd9 100644
--- a/pref.h
+++ b/pref.h
@@ -32,6 +32,12 @@ typedef struct {
 	enum taxonomy_category category[3];
 } geocoding_prefs_t;
 
+enum deco_mode {
+	BUEHLMANN,
+	RECREATIONAL,
+	VPMB
+};
+
 struct preferences {
 	const char *divelist_font;
 	const char *default_filename;
@@ -89,7 +95,6 @@ struct preferences {
 	bool display_runtime;
 	bool display_duration;
 	bool display_transitions;
-	bool recreational_mode;
 	bool safetystop;
 	bool switch_at_req_stop;
 	int reserve_gas;
@@ -110,6 +115,7 @@ struct preferences {
 	short cloud_verification_status;
 	bool cloud_background_sync;
 	geocoding_prefs_t geocoding;
+	enum deco_mode deco_mode;
 };
 enum unit_system_values {
 	METRIC,
diff --git a/qt-models/diveplannermodel.cpp b/qt-models/diveplannermodel.cpp
index adbcdba9a..e28b4d146 100644
--- a/qt-models/diveplannermodel.cpp
+++ b/qt-models/diveplannermodel.cpp
@@ -449,10 +449,10 @@ void DivePlannerPointsModel::setDisplayTransitions(bool value)
 	emit dataChanged(createIndex(0, 0), createIndex(rowCount() - 1, COLUMNS - 1));
 }
 
-void DivePlannerPointsModel::setRecreationalMode(bool value)
+void DivePlannerPointsModel::setDecoMode(int mode)
 {
-	prefs.recreational_mode = value;
-	emit recreationChanged(value);
+	prefs.deco_mode = deco_mode(mode);
+	emit recreationChanged(mode == int(RECREATIONAL));
 	emit dataChanged(createIndex(0, 0), createIndex(rowCount() - 1, COLUMNS -1));
 }
 
diff --git a/qt-models/diveplannermodel.h b/qt-models/diveplannermodel.h
index 422499084..a8524393a 100644
--- a/qt-models/diveplannermodel.h
+++ b/qt-models/diveplannermodel.h
@@ -78,7 +78,7 @@ slots:
 	void setDisplayRuntime(bool value);
 	void setDisplayDuration(bool value);
 	void setDisplayTransitions(bool value);
-	void setRecreationalMode(bool value);
+	void setDecoMode(int mode);
 	void setSafetyStop(bool value);
 	void savePlan();
 	void saveDuplicatePlan();
diff --git a/qt-ui/diveplanner.cpp b/qt-ui/diveplanner.cpp
index 8b0e93cd8..684ef44aa 100644
--- a/qt-ui/diveplanner.cpp
+++ b/qt-ui/diveplanner.cpp
@@ -243,7 +243,7 @@ PlannerSettingsWidget::PlannerSettingsWidget(QWidget *parent, Qt::WindowFlags f)
 	prefs.display_duration = s.value("display_duration", prefs.display_duration).toBool();
 	prefs.display_runtime = s.value("display_runtime", prefs.display_runtime).toBool();
 	prefs.display_transitions = s.value("display_transitions", prefs.display_transitions).toBool();
-	prefs.recreational_mode = s.value("recreational_mode", prefs.recreational_mode).toBool();
+	prefs.deco_mode = deco_mode(s.value("deco_mode", prefs.deco_mode).toInt());
 	prefs.safetystop = s.value("safetystop", prefs.safetystop).toBool();
 	prefs.reserve_gas = s.value("reserve_gas", prefs.reserve_gas).toInt();
 	prefs.ascrate75 = s.value("ascrate75", prefs.ascrate75).toInt();
@@ -269,7 +269,6 @@ PlannerSettingsWidget::PlannerSettingsWidget(QWidget *parent, Qt::WindowFlags f)
 	ui.display_duration->setChecked(prefs.display_duration);
 	ui.display_runtime->setChecked(prefs.display_runtime);
 	ui.display_transitions->setChecked(prefs.display_transitions);
-	ui.recreational_mode->setChecked(prefs.recreational_mode);
 	ui.safetystop->setChecked(prefs.safetystop);
 	ui.reserve_gas->setValue(prefs.reserve_gas / 1000);
 	ui.bottompo2->setValue(prefs.bottompo2 / 1000.0);
@@ -278,17 +277,30 @@ PlannerSettingsWidget::PlannerSettingsWidget(QWidget *parent, Qt::WindowFlags f)
 	ui.drop_stone_mode->setChecked(prefs.drop_stone_mode);
 	ui.switch_at_req_stop->setChecked(prefs.switch_at_req_stop);
 	ui.min_switch_duration->setValue(prefs.min_switch_duration / 60);
+	ui.recreational_deco->setChecked(prefs.deco_mode == RECREATIONAL);
+	ui.buehlmann_deco->setChecked(prefs.deco_mode == BUEHLMANN);
+	ui.vpmb_deco->setChecked(prefs.deco_mode == VPMB);
+
 	// should be the same order as in dive_comp_type!
 	rebreater_modes << tr("Open circuit") << tr("CCR") << tr("pSCR");
 	ui.rebreathermode->insertItems(0, rebreater_modes);
 
+	modeMapper = new QSignalMapper(this);
+	connect(modeMapper, SIGNAL(mapped(int)) , plannerModel, SLOT(setDecoMode(int)));
+	modeMapper->setMapping(ui.recreational_deco, int(RECREATIONAL));
+	modeMapper->setMapping(ui.buehlmann_deco, int(BUEHLMANN));
+	modeMapper->setMapping(ui.vpmb_deco, int(VPMB));
+	
+	connect(ui.recreational_deco, SIGNAL(clicked()), modeMapper, SLOT(map()));
+	connect(ui.buehlmann_deco, SIGNAL(clicked()), modeMapper, SLOT(map()));
+	connect(ui.vpmb_deco, SIGNAL(clicked()), modeMapper, SLOT(map()));
+	
 	connect(ui.lastStop, SIGNAL(toggled(bool)), plannerModel, SLOT(setLastStop6m(bool)));
 	connect(ui.verbatim_plan, SIGNAL(toggled(bool)), plannerModel, SLOT(setVerbatim(bool)));
 	connect(ui.display_duration, SIGNAL(toggled(bool)), plannerModel, SLOT(setDisplayDuration(bool)));
 	connect(ui.display_runtime, SIGNAL(toggled(bool)), plannerModel, SLOT(setDisplayRuntime(bool)));
 	connect(ui.display_transitions, SIGNAL(toggled(bool)), plannerModel, SLOT(setDisplayTransitions(bool)));
 	connect(ui.safetystop, SIGNAL(toggled(bool)), plannerModel, SLOT(setSafetyStop(bool)));
-	connect(ui.recreational_mode, SIGNAL(toggled(bool)), plannerModel, SLOT(setRecreationalMode(bool)));
 	connect(ui.reserve_gas, SIGNAL(valueChanged(int)), plannerModel, SLOT(setReserveGas(int)));
 	connect(ui.ascRate75, SIGNAL(valueChanged(int)), this, SLOT(setAscRate75(int)));
 	connect(ui.ascRate75, SIGNAL(valueChanged(int)), plannerModel, SLOT(emitDataChanged()));
@@ -341,7 +353,6 @@ PlannerSettingsWidget::~PlannerSettingsWidget()
 	s.setValue("display_duration", prefs.display_duration);
 	s.setValue("display_runtime", prefs.display_runtime);
 	s.setValue("display_transitions", prefs.display_transitions);
-	s.setValue("recreational_mode", prefs.recreational_mode);
 	s.setValue("safetystop", prefs.safetystop);
 	s.setValue("reserve_gas", prefs.reserve_gas);
 	s.setValue("ascrate75", prefs.ascrate75);
@@ -357,6 +368,7 @@ PlannerSettingsWidget::~PlannerSettingsWidget()
 	s.setValue("min_switch_duration", prefs.min_switch_duration);
 	s.setValue("bottomsac", prefs.bottomsac);
 	s.setValue("decosac", prefs.decosac);
+	s.setValue("deco_mode", int(prefs.deco_mode));
 	s.endGroup();
 }
 
diff --git a/qt-ui/diveplanner.h b/qt-ui/diveplanner.h
index 8c7ff9c46..c4210aadf 100644
--- a/qt-ui/diveplanner.h
+++ b/qt-ui/diveplanner.h
@@ -5,6 +5,7 @@
 #include <QAbstractTableModel>
 #include <QAbstractButton>
 #include <QDateTime>
+#include <QSignalMapper>
 
 #include "dive.h"
 
@@ -84,6 +85,7 @@ slots:
 private:
 	Ui::plannerSettingsWidget ui;
 	void updateUnitsUI();
+	QSignalMapper *modeMapper;
 };
 
 #include "ui_plannerDetails.h"
diff --git a/qt-ui/plannerSettings.ui b/qt-ui/plannerSettings.ui
index 55ca83a9b..5916b2777 100644
--- a/qt-ui/plannerSettings.ui
+++ b/qt-ui/plannerSettings.ui
@@ -359,13 +359,30 @@
            </widget>
           </item>
           <item row="0" column="1">
-           <widget class="QCheckBox" name="recreational_mode">
+           <widget class="QRadioButton" name="recreational_deco">
             <property name="text">
              <string>Recreational mode</string>
             </property>
            </widget>
           </item>
           <item row="1" column="1">
+           <widget class="QRadioButton" name="buehlmann_deco">
+            <property name="text">
+             <string>Buehlmann deco</string>
+            </property>
+            <property name="checked">
+             <bool>true</bool>
+            </property>
+           </widget>
+          </item>
+          <item row="2" column="1">
+           <widget class="QRadioButton" name="vpmb_deco">
+            <property name="text">
+             <string>VPM-B deco</string>
+            </property>
+           </widget>
+          </item>
+          <item row="3" column="1">
            <widget class="QCheckBox" name="safetystop">
             <property name="text">
              <string>Safety stop</string>
@@ -375,14 +392,14 @@
             </property>
            </widget>
           </item>
-          <item row="2" column="1">
+          <item row="4" column="1">
            <widget class="QLabel" name="label_3">
             <property name="text">
              <string>Reserve gas</string>
             </property>
            </widget>
           </item>
-          <item row="2" column="2">
+          <item row="4" column="2">
            <widget class="QSpinBox" name="reserve_gas">
             <property name="suffix">
              <string>bar</string>
diff --git a/subsurfacestartup.c b/subsurfacestartup.c
index 3005e7e04..17bd43e45 100644
--- a/subsurfacestartup.c
+++ b/subsurfacestartup.c
@@ -54,7 +54,6 @@ struct preferences default_prefs = {
 	.display_runtime = true,
 	.display_duration = true,
 	.display_transitions = true,
-	.recreational_mode = false,
 	.safetystop = true,
 	.bottomsac = 20000,
 	.decosac = 17000,
@@ -75,7 +74,8 @@ struct preferences default_prefs = {
 		.parse_dive_without_gps = false,
 		.tag_existing_dives = false,
 		.category = { 0 }
-	}
+	},
+	.deco_mode = BUEHLMANN
 };
 
 int run_survey;