Hallo Robert,

Am 17.09.2017 um 21:44 schrieb Robert Helling:
> Stefan,
>
> could you please try if applying this patch removes the spikes? My suspicion 
> is that those are not properly reset between dives.
>
Just one note and a question about your patch for resetting some
additional VPM state variables (I attached it again):
We meanwhile both understood that this didn't solve the strange spikes
at the beginning of the deco ceiling.
But at least for one of my logged dives it removed a kind of
"randomness" in the deco calculation. Also the patch seems not to change
deco calculation in a negative way means creating different or wrong
results - at least not for the few dives I tested.
I expect that tests also still pass?

So if you are also sure that it's the correct approach from theoretical
point of view, the patch should be submitted to master.

Best regards
Stefan


-- 

Stefan Fuchs
E-Mail: [email protected] <mailto:[email protected]>

>From 92f648a9c79e65dbef11d18c8e07d7cce2a1b1b4 Mon Sep 17 00:00:00 2001
From: "Robert C. Helling" <[email protected]>
Date: Sun, 17 Sep 2017 21:06:44 +0200
Subject: [PATCH] More VPMB state in special structure
To: [email protected]

... and reset deco information in profile ceiling computation.

Signed-off-by: Robert C. Helling <[email protected]>
---
 core/deco.c        | 11 +++++++----
 core/dive.h        |  2 ++
 core/planner.c     | 14 ++++++--------
 core/profile.c     |  3 +--
 tests/testplan.cpp | 24 +++++++++++-------------
 5 files changed, 27 insertions(+), 27 deletions(-)

diff --git a/core/deco.c b/core/deco.c
index 469679684..4f4811368 100644
--- a/core/deco.c
+++ b/core/deco.c
@@ -34,8 +34,6 @@
 extern bool in_planner();
 extern int plot_depth;
 
-extern pressure_t first_ceiling_pressure;
-
 //! Option structure for Buehlmann decompression.
 struct buehlmann_config {
        double satmult;                 //! safety at inert gas accumulation as 
percentage of effect (more than 100).
@@ -199,7 +197,7 @@ double solve_cubic2(double B, double C)
 
 double update_gradient(double next_stop_pressure, double first_gradient)
 {
-       double B = cube(first_gradient) / (first_ceiling_pressure.mbar / 1000.0 
+ first_gradient);
+       double B = cube(first_gradient) / 
(deco_state->first_ceiling_pressure.mbar / 1000.0 + first_gradient);
        double C = next_stop_pressure * B;
 
        double new_gradient = solve_cubic2(B, C);
@@ -213,7 +211,7 @@ double vpmb_tolerated_ambient_pressure(double 
reference_pressure, int ci)
 {
        double n2_gradient, he_gradient, total_gradient;
 
-       if (reference_pressure >= first_ceiling_pressure.mbar / 1000.0 || 
!first_ceiling_pressure.mbar) {
+       if (reference_pressure >= deco_state->first_ceiling_pressure.mbar / 
1000.0 || !deco_state->first_ceiling_pressure.mbar) {
                n2_gradient = deco_state->bottom_n2_gradient[ci];
                he_gradient = deco_state->bottom_he_gradient[ci];
        } else {
@@ -298,6 +296,7 @@ double tissue_tolerance_calc(const struct dive *dive, 
double pressure)
                // We are doing ok if the gradient was computed within ten 
centimeters of the ceiling.
                } while (fabs(ret_tolerance_limit_ambient_pressure - 
reference_pressure) > 0.01);
 
+               printf("initial n2 graditent 0: %f\n", 
deco_state->initial_n2_gradient[0]);
                if (plot_depth) {
                        ++sum1;
                        sumx += plot_depth;
@@ -537,6 +536,8 @@ void clear_vpmb_state() {
                deco_state->max_he_crushing_pressure[ci] = 0.0;
        }
        deco_state->max_ambient_pressure = 0;
+       deco_state->first_ceiling_pressure.mbar = 0;
+       deco_state->max_bottom_ceiling_pressure.mbar = 0;
 }
 
 void clear_deco(double surface_pressure)
@@ -578,6 +579,8 @@ void restore_deco_state(struct deco_state *data, bool 
keep_vpmb_state)
                        data->initial_n2_gradient[ci] = 
deco_state->initial_n2_gradient[ci];
                        data->initial_he_gradient[ci] = 
deco_state->initial_he_gradient[ci];
                }
+               data->first_ceiling_pressure = 
deco_state->first_ceiling_pressure;
+               data->max_bottom_ceiling_pressure = 
deco_state->max_bottom_ceiling_pressure;
        }
        *deco_state = *data;
 
diff --git a/core/dive.h b/core/dive.h
index 43464118f..d128fc79a 100644
--- a/core/dive.h
+++ b/core/dive.h
@@ -858,6 +858,8 @@ struct deco_state {
 
        double initial_n2_gradient[16];
        double initial_he_gradient[16];
+       pressure_t first_ceiling_pressure;
+       pressure_t max_bottom_ceiling_pressure;
        int ci_pointing_to_guiding_tissue;
        double gf_low_pressure_this_dive;
 
diff --git a/core/planner.c b/core/planner.c
index da6ed7d73..c4d478a03 100644
--- a/core/planner.c
+++ b/core/planner.c
@@ -39,8 +39,6 @@ extern double regressiona();
 extern double regressionb();
 extern void reset_regression();
 
-pressure_t first_ceiling_pressure, max_bottom_ceiling_pressure = {};
-
 char *disclaimer;
 int plot_depth = 0;
 #if DEBUG_PLAN
@@ -184,8 +182,8 @@ unsigned int tissue_at_end(struct dive *dive, struct 
deco_state **cached_datap)
                                                                                
dive,
                                                                                
1),
                                                                dive);
-                       if (ceiling_pressure.mbar > 
max_bottom_ceiling_pressure.mbar)
-                               max_bottom_ceiling_pressure.mbar = 
ceiling_pressure.mbar;
+                       if (ceiling_pressure.mbar > 
deco_state->max_bottom_ceiling_pressure.mbar)
+                               deco_state->max_bottom_ceiling_pressure.mbar = 
ceiling_pressure.mbar;
                }
 
                interpolate_transition(dive, t0, t1, lastdepth, sample->depth, 
&gas, setpoint);
@@ -691,7 +689,7 @@ bool plan(struct diveplan *diveplan, struct dive *dive, int 
timestep, struct dec
                diveplan->surface_pressure = SURFACE_PRESSURE;
        dive->surface_pressure.mbar = diveplan->surface_pressure;
        clear_deco(dive->surface_pressure.mbar / 1000.0);
-       max_bottom_ceiling_pressure.mbar = first_ceiling_pressure.mbar = 0;
+       deco_state->max_bottom_ceiling_pressure.mbar = 
deco_state->first_ceiling_pressure.mbar = 0;
        create_dive_from_plan(diveplan, dive, is_planner);
 
        // Do we want deco stop array in metres or feet?
@@ -859,14 +857,14 @@ bool plan(struct diveplan *diveplan, struct dive *dive, 
int timestep, struct dec
                breaktime = -1;
                breakcylinder = 0;
                o2time = 0;
-               first_ceiling_pressure.mbar = 
depth_to_mbar(deco_allowed_depth(tissue_tolerance_calc(dive,
+               deco_state->first_ceiling_pressure.mbar = 
depth_to_mbar(deco_allowed_depth(tissue_tolerance_calc(dive,
                                                                                
                     depth_to_bar(depth, dive)),
                                                                            
diveplan->surface_pressure / 1000.0,
                                                                            
dive,
                                                                            1),
                                                         dive);
-               if (max_bottom_ceiling_pressure.mbar > 
first_ceiling_pressure.mbar)
-                       first_ceiling_pressure.mbar = 
max_bottom_ceiling_pressure.mbar;
+               if (deco_state->max_bottom_ceiling_pressure.mbar > 
deco_state->first_ceiling_pressure.mbar)
+                       deco_state->first_ceiling_pressure.mbar = 
deco_state->max_bottom_ceiling_pressure.mbar;
 
                last_ascend_rate = ascent_velocity(depth, avg_depth, 
bottom_time);
                if ((current_cylinder = get_gasidx(dive, &gas)) == -1) {
diff --git a/core/profile.c b/core/profile.c
index 4ff27fcd7..813e69fc2 100644
--- a/core/profile.c
+++ b/core/profile.c
@@ -31,7 +31,6 @@ static struct plot_data *last_pi_entry_new = NULL;
 void populate_pressure_information(struct dive *, struct divecomputer *, 
struct plot_info *, int);
 
 extern bool in_planner();
-extern pressure_t first_ceiling_pressure;
 
 #ifdef DEBUG_PI
 /* debugging tool - not normally used */
@@ -1015,7 +1014,7 @@ void calculate_deco_information(struct dive *dive, struct 
divecomputer *dc, stru
                                        if  (current_ceiling > first_ceiling) {
                                                time_deep_ceiling = t1;
                                                first_ceiling = current_ceiling;
-                                               first_ceiling_pressure.mbar = 
depth_to_mbar(first_ceiling, dive);
+                                               
deco_state->first_ceiling_pressure.mbar = depth_to_mbar(first_ceiling, dive);
                                                if (first_iteration) {
                                                        
nuclear_regeneration(t1);
                                                        vpmb_start_gradient();
diff --git a/tests/testplan.cpp b/tests/testplan.cpp
index 5dbcf89b9..0fdb51abb 100644
--- a/tests/testplan.cpp
+++ b/tests/testplan.cpp
@@ -13,8 +13,6 @@
 struct decostop stoptable[60];
 extern bool plan(struct diveplan *diveplan, struct dive *dive, int timestep, 
struct decostop *decostoptable, struct deco_state **cached_datap, bool 
is_planner, bool show_disclaimer);
 
-extern pressure_t first_ceiling_pressure;
-
 void setupPrefs()
 {
        copy_prefs(&default_prefs, &prefs);
@@ -458,7 +456,7 @@ void TestPlan::testVpmbMetric45m30minTx()
        while(!dp->minimum_gas.mbar && dp->next) dp = dp->next;
        QCOMPARE(lrint(dp->minimum_gas.mbar / 1000.0), 108l);
        // print first ceiling
-       printf("First ceiling %.1f m\n", 
(mbar_to_depth(first_ceiling_pressure.mbar, &displayed_dive) * 0.001));
+//     printf("First ceiling %.1f m\n", 
(mbar_to_depth(deco_state->first_ceiling_pressure.mbar, &displayed_dive) * 
0.001));
        // check benchmark run time of 141 minutes, and known Subsurface 
runtime of 139 minutes
        //QVERIFY(compareDecoTime(displayed_dive.dc.duration.seconds, 141u * 
60u + 20u, 139u * 60u + 20u));
 }
@@ -488,7 +486,7 @@ void TestPlan::testVpmbMetric60m10minTx()
        while(!dp->minimum_gas.mbar && dp->next) dp = dp->next;
        QCOMPARE(lrint(dp->minimum_gas.mbar / 1000.0), 162l);
        // print first ceiling
-       printf("First ceiling %.1f m\n", 
(mbar_to_depth(first_ceiling_pressure.mbar, &displayed_dive) * 0.001));
+//     printf("First ceiling %.1f m\n", 
(mbar_to_depth(deco_state->first_ceiling_pressure.mbar, &displayed_dive) * 
0.001));
        // check benchmark run time of 141 minutes, and known Subsurface 
runtime of 139 minutes
        //QVERIFY(compareDecoTime(displayed_dive.dc.duration.seconds, 141u * 
60u + 20u, 139u * 60u + 20u));
 }
@@ -518,7 +516,7 @@ void TestPlan::testVpmbMetric60m30minAir()
        while(!dp->minimum_gas.mbar && dp->next) dp = dp->next;
        QCOMPARE(lrint(dp->minimum_gas.mbar / 1000.0), 180l);
        // print first ceiling
-       printf("First ceiling %.1f m\n", 
(mbar_to_depth(first_ceiling_pressure.mbar, &displayed_dive) * 0.001));
+//     printf("First ceiling %.1f m\n", 
(mbar_to_depth(deco_state->first_ceiling_pressure.mbar, &displayed_dive) * 
0.001));
        // check benchmark run time of 141 minutes, and known Subsurface 
runtime of 139 minutes
        QVERIFY(compareDecoTime(displayed_dive.dc.duration.seconds, 141u * 60u 
+ 20u, 139u * 60u + 20u));
 }
@@ -548,7 +546,7 @@ void TestPlan::testVpmbMetric60m30minEan50()
        while(!dp->minimum_gas.mbar && dp->next) dp = dp->next;
        QCOMPARE(lrint(dp->minimum_gas.mbar / 1000.0), 155l);
        // print first ceiling
-       printf("First ceiling %.1f m\n", 
(mbar_to_depth(first_ceiling_pressure.mbar, &displayed_dive) * 0.001));
+//     printf("First ceiling %.1f m\n", 
(mbar_to_depth(deco_state->first_ceiling_pressure.mbar, &displayed_dive) * 
0.001));
        // check first gas change to EAN50 at 21m
        struct event *ev = displayed_dive.dc.events;
        QVERIFY(ev != NULL);
@@ -584,7 +582,7 @@ void TestPlan::testVpmbMetric60m30minTx()
        while(!dp->minimum_gas.mbar && dp->next) dp = dp->next;
        QCOMPARE(lrint(dp->minimum_gas.mbar / 1000.0), 159l);
        // print first ceiling
-       printf("First ceiling %.1f m\n", 
(mbar_to_depth(first_ceiling_pressure.mbar, &displayed_dive) * 0.001));
+//     printf("First ceiling %.1f m\n", 
(mbar_to_depth(deco_state->first_ceiling_pressure.mbar, &displayed_dive) * 
0.001));
        // check first gas change to EAN50 at 21m
        struct event *ev = displayed_dive.dc.events;
        QVERIFY(ev != NULL);
@@ -620,7 +618,7 @@ void TestPlan::testVpmbMetric100m60min()
        while(!dp->minimum_gas.mbar && dp->next) dp = dp->next;
        QCOMPARE(lrint(dp->minimum_gas.mbar / 1000.0), 157l);
        // print first ceiling
-       printf("First ceiling %.1f m\n", 
(mbar_to_depth(first_ceiling_pressure.mbar, &displayed_dive) * 0.001));
+//     printf("First ceiling %.1f m\n", 
(mbar_to_depth(deco_state->first_ceiling_pressure.mbar, &displayed_dive) * 
0.001));
        // check first gas change to EAN50 at 21m
        struct event *ev = displayed_dive.dc.events;
        QVERIFY(ev != NULL);
@@ -662,7 +660,7 @@ void TestPlan::testVpmbMetricMultiLevelAir()
        while(!dp->minimum_gas.mbar && dp->next) dp = dp->next;
        QCOMPARE(lrint(dp->minimum_gas.mbar / 1000.0), 101l);
        // print first ceiling
-       printf("First ceiling %.1f m\n", 
(mbar_to_depth(first_ceiling_pressure.mbar, &displayed_dive) * 0.001));
+//     printf("First ceiling %.1f m\n", 
(mbar_to_depth(deco_state->first_ceiling_pressure.mbar, &displayed_dive) * 
0.001));
        // check benchmark run time of 167 minutes, and known Subsurface 
runtime of 169 minutes
        QVERIFY(compareDecoTime(displayed_dive.dc.duration.seconds, 167u * 60u 
+ 20u, 169u * 60u + 20u));
 }
@@ -692,7 +690,7 @@ void TestPlan::testVpmbMetric100m10min()
        while(!dp->minimum_gas.mbar && dp->next) dp = dp->next;
        QCOMPARE(lrint(dp->minimum_gas.mbar / 1000.0), 175l);
        // print first ceiling
-       printf("First ceiling %.1f m\n", 
(mbar_to_depth(first_ceiling_pressure.mbar, &displayed_dive) * 0.001));
+//     printf("First ceiling %.1f m\n", 
(mbar_to_depth(deco_state->first_ceiling_pressure.mbar, &displayed_dive) * 
0.001));
        // check first gas change to EAN50 at 21m
        struct event *ev = displayed_dive.dc.events;
        QVERIFY(ev != NULL);
@@ -738,7 +736,7 @@ void TestPlan::testVpmbMetricRepeat()
        while(!dp->minimum_gas.mbar && dp->next) dp = dp->next;
        QCOMPARE(lrint(dp->minimum_gas.mbar / 1000.0), 61l);
        // print first ceiling
-       printf("First ceiling %.1f m\n", 
(mbar_to_depth(first_ceiling_pressure.mbar, &displayed_dive) * 0.001));
+//     printf("First ceiling %.1f m\n", 
(mbar_to_depth(deco_state->first_ceiling_pressure.mbar, &displayed_dive) * 
0.001));
        // check benchmark run time of 27 minutes, and known Subsurface runtime 
of 28 minutes
        QVERIFY(compareDecoTime(displayed_dive.dc.duration.seconds, 27u * 60u + 
20u, 27u * 60u + 20u));
 
@@ -758,7 +756,7 @@ void TestPlan::testVpmbMetricRepeat()
        while(!dp->minimum_gas.mbar && dp->next) dp = dp->next;
        QCOMPARE(lrint(dp->minimum_gas.mbar / 1000.0), 80l);
        // print first ceiling
-       printf("First ceiling %.1f m\n", 
(mbar_to_depth(first_ceiling_pressure.mbar, &displayed_dive) * 0.001));
+//     printf("First ceiling %.1f m\n", 
(mbar_to_depth(deco_state->first_ceiling_pressure.mbar, &displayed_dive) * 
0.001));
        // check first gas change to 21/35 at 66m
        struct event *ev = displayed_dive.dc.events;
        QVERIFY(ev != NULL);
@@ -794,7 +792,7 @@ void TestPlan::testVpmbMetricRepeat()
        while(!dp->minimum_gas.mbar && dp->next) dp = dp->next;
        QCOMPARE(lrint(dp->minimum_gas.mbar / 1000.0), 61l);
        // print first ceiling
-       printf("First ceiling %.1f m\n", 
(mbar_to_depth(first_ceiling_pressure.mbar, &displayed_dive) * 0.001));
+//     printf("First ceiling %.1f m\n", 
(mbar_to_depth(deco_state->first_ceiling_pressure.mbar, &displayed_dive) * 
0.001));
 
        // check runtime is exactly the same as the first time
        int finalDiveRunTimeSeconds = displayed_dive.dc.duration.seconds;
-- 
2.11.0 (Apple Git-81)

_______________________________________________
subsurface mailing list
[email protected]
http://lists.subsurface-divelog.org/cgi-bin/mailman/listinfo/subsurface

Reply via email to