Thanks again Robert, On 30 October 2015 at 08:44, Robert C. Helling <[email protected]> wrote:
> > On 29 Oct 2015, at 12:53, Rick Walsh <[email protected]> wrote: > > Comments added to commit > > > What I already said on IRC: > > It should probably be > > if (prefs.deco_mode == VPMB && !in_planner() && entry->ceiling >= > first_ceiling && first_iteration == true) { > > > > with a good initialisation for first_ceiling so this is true when it is > first called (where no ceiling so far has been computed). Otherwise, at > least for my real dives, the condition is never true (a breakpoint on the > next statement was never reached). > Thanks for picking this up. It's somewhat surprising it nearly worked before. Actually, the comparison should be (entry - 1)->ceiling >= first_ceiling, because we want to run it after the first_ceiling has been set. To make things a bit simpler, I have shifted it down to after first_ceiling has been set. While we *could* recalculate first_ceiling at this point based on the updated gradients, but there's no need as it will be picked up in the next CVA iteration. > Then, have a look at this file: > > > > This is the same dive planned once with VPM-B and once with Buehlmann > 30/80. For the VPM-B the ceiling looks ok. The VPM-B version has a runtime > of 55 as compared to that of Buehlmann which has 45. So I would expect to > have in the end a TTS of about 10min. But it shows 37min (and a ceiling > violation much earlier). Maybe that is what it is. But I am a bit surprised > by thus behaviour. > > The problem was that the function calculate_ndl_tts keeps adding to entry->tts_calc. That's fine for stepping through the dive the first time, but with iterating for the CVA, it needs to be reset to zero. The attached patch does that. Cheers, Rick
From 294135384d7785edeb358bb9fe8dd4a73a2a0aff Mon Sep 17 00:00:00 2001 From: Rick Walsh <[email protected]> Date: Sun, 25 Oct 2015 11:26:15 +1100 Subject: [PATCH 1/2] Calculate VPM-B ceiling outside of planner Doing VPM-B calculations for dives outside of the planner has not been possible because real dive data does not record either: - reference pressure for the Boyle's law compensation (i.e. first_stop_pressure), or - deco_time for the vpmb_next_gradient function used to do the CVA calculations However, we can infer these values to be: - first_stop_pressure is the deepest ceiling in the dive - deco_time is dive time from the deepest ceiling until the ceiling clears (or would have cleared if the diver finished their deco obligations) With these assumptions, the CVA converges rapidly. Signed-off-by: Rick Walsh <[email protected]> --- deco.c | 4 +- profile.c | 183 ++++++++++++++++++++++++++++++++++++++++++-------------------- 2 files changed, 128 insertions(+), 59 deletions(-) diff --git a/deco.c b/deco.c index 86acc03..506d8b5 100644 --- a/deco.c +++ b/deco.c @@ -244,7 +244,7 @@ double tissue_tolerance_calc(const struct dive *dive, double pressure) double lowest_ceiling = 0.0; double tissue_lowest_ceiling[16]; - if (prefs.deco_mode != VPMB || !in_planner()) { + if (prefs.deco_mode != VPMB) { for (ci = 0; ci < 16; ci++) { tissue_inertgas_saturation[ci] = tissue_n2_sat[ci] + tissue_he_sat[ci]; buehlmann_inertgas_a[ci] = ((buehlmann_N2_a[ci] * tissue_n2_sat[ci]) + (buehlmann_He_a[ci] * tissue_he_sat[ci])) / tissue_inertgas_saturation[ci]; @@ -509,7 +509,7 @@ void add_segment(double pressure, const struct gasmix *gasmix, int period_in_sec tissue_n2_sat[ci] += n2_satmult * pn2_oversat * n2_f; tissue_he_sat[ci] += he_satmult * phe_oversat * he_f; } - if(prefs.deco_mode == VPMB && in_planner()) + if(prefs.deco_mode == VPMB) calc_crushing_pressure(pressure); return; } diff --git a/profile.c b/profile.c index d39133c..edba51e 100644 --- a/profile.c +++ b/profile.c @@ -28,6 +28,9 @@ unsigned int dc_number = 0; 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 */ static void dump_pi(struct plot_info *pi) @@ -866,6 +869,8 @@ static void calculate_ndl_tts(struct plot_data *entry, struct dive *dive, double int ascent_depth = entry->depth; /* at what time should we give up and say that we got enuff NDL? */ int cylinderindex = entry->cylinderindex; + /* If iterating through a dive, entry->tts_calc needs to be reset */ + entry->tts_calc = 0; /* If we don't have a ceiling yet, calculate ndl. Don't try to calculate * a ndl for lower values than 3m it would take forever */ @@ -927,69 +932,133 @@ static void calculate_ndl_tts(struct plot_data *entry, struct dive *dive, double */ void calculate_deco_information(struct dive *dive, struct divecomputer *dc, struct plot_info *pi, bool print_mode) { - int i; + int i, count_iteration = 0; double surface_pressure = (dc->surface_pressure.mbar ? dc->surface_pressure.mbar : get_surface_pressure_in_mbar(dive, true)) / 1000.0; int last_ndl_tts_calc_time = 0; - for (i = 1; i < pi->nr; i++) { - struct plot_data *entry = pi->entry + i; - int j, t0 = (entry - 1)->sec, t1 = entry->sec; - int time_stepsize = 20; - - entry->ambpressure = depth_to_bar(entry->depth, dive); - entry->gfline = MAX((double)prefs.gflow, (entry->ambpressure - surface_pressure) / (gf_low_pressure_this_dive - surface_pressure) * - (prefs.gflow - prefs.gfhigh) + - prefs.gfhigh) * - (100.0 - AMB_PERCENTAGE) / 100.0 + AMB_PERCENTAGE; - if (t0 > t1) { - fprintf(stderr, "non-monotonous dive stamps %d %d\n", t0, t1); - int xchg = t1; - t1 = t0; - t0 = xchg; - } - if (t0 != t1 && t1 - t0 < time_stepsize) - time_stepsize = t1 - t0; - for (j = t0 + time_stepsize; j <= t1; j += time_stepsize) { - int depth = interpolate(entry[-1].depth, entry[0].depth, j - t0, t1 - t0); - add_segment(depth_to_bar(depth, dive), - &dive->cylinder[entry->cylinderindex].gasmix, time_stepsize, entry->o2pressure.mbar, dive, entry->sac); - if ((t1 - j < time_stepsize) && (j < t1)) - time_stepsize = t1 - j; - } - if (t0 == t1) - entry->ceiling = (entry - 1)->ceiling; - else - entry->ceiling = deco_allowed_depth(tissue_tolerance_calc(dive, depth_to_bar(entry->depth, dive)), surface_pressure, dive, !prefs.calcceiling3m); - for (j = 0; j < 16; j++) { - double m_value = buehlmann_inertgas_a[j] + entry->ambpressure / buehlmann_inertgas_b[j]; - entry->ceilings[j] = deco_allowed_depth(tolerated_by_tissue[j], surface_pressure, dive, 1); - entry->percentages[j] = tissue_inertgas_saturation[j] < entry->ambpressure ? - tissue_inertgas_saturation[j] / entry->ambpressure * AMB_PERCENTAGE : - AMB_PERCENTAGE + (tissue_inertgas_saturation[j] - entry->ambpressure) / (m_value - entry->ambpressure) * (100.0 - AMB_PERCENTAGE); - } + int first_ceiling = 0; + bool first_iteration = true; + int final_tts = 0 , time_clear_ceiling = 0, time_deep_ceiling = 0, deco_time = 0, prev_deco_time = 10000000; + char *cache_data_initial = NULL; + /* For VPM-B outside the planner, cache the initial deco state for CVA iterations */ + if (prefs.deco_mode == VPMB && !in_planner()) + cache_deco_state(&cache_data_initial); + /* For VPM-B outside the planner, iterate until deco time converges (usually one or two iterations after the initial) + * Set maximum number of iterations to 10 just in case */ + while ((abs(prev_deco_time - deco_time) >= 30) && (count_iteration < 10)) { + for (i = 1; i < pi->nr; i++) { + struct plot_data *entry = pi->entry + i; + int j, t0 = (entry - 1)->sec, t1 = entry->sec; + int time_stepsize = 20; + + entry->ambpressure = depth_to_bar(entry->depth, dive); + entry->gfline = MAX((double)prefs.gflow, (entry->ambpressure - surface_pressure) / (gf_low_pressure_this_dive - surface_pressure) * + (prefs.gflow - prefs.gfhigh) + + prefs.gfhigh) * + (100.0 - AMB_PERCENTAGE) / 100.0 + AMB_PERCENTAGE; + if (t0 > t1) { + fprintf(stderr, "non-monotonous dive stamps %d %d\n", t0, t1); + int xchg = t1; + t1 = t0; + t0 = xchg; + } + if (t0 != t1 && t1 - t0 < time_stepsize) + time_stepsize = t1 - t0; + for (j = t0 + time_stepsize; j <= t1; j += time_stepsize) { + int depth = interpolate(entry[-1].depth, entry[0].depth, j - t0, t1 - t0); + add_segment(depth_to_bar(depth, dive), + &dive->cylinder[entry->cylinderindex].gasmix, time_stepsize, entry->o2pressure.mbar, dive, entry->sac); + if ((t1 - j < time_stepsize) && (j < t1)) + time_stepsize = t1 - j; + } + if (t0 == t1) { + entry->ceiling = (entry - 1)->ceiling; + } else { + /* Keep updating the VPM-B gradients until the start of the ascent phase of the dive. */ + if (prefs.deco_mode == VPMB && !in_planner() && (entry - 1)->ceiling >= first_ceiling && first_iteration == true) { + nuclear_regeneration(t1); + vpmb_start_gradient(); + /* For CVA calculations, start by guessing deco time = dive time remaining */ + deco_time = pi->maxtime - t1; + vpmb_next_gradient(deco_time, surface_pressure / 1000.0); + } + entry->ceiling = deco_allowed_depth(tissue_tolerance_calc(dive, depth_to_bar(entry->depth, dive)), surface_pressure, dive, !prefs.calcceiling3m); + /* If using VPM-B outside the planner, take first_ceiling_pressure as the deepest ceiling */ + if (prefs.deco_mode == VPMB && !in_planner()) { + if (entry->ceiling >= first_ceiling) { + time_deep_ceiling = t1; + first_ceiling = entry->ceiling; + first_ceiling_pressure.mbar = depth_to_mbar(first_ceiling, dive); + if (first_iteration) { + nuclear_regeneration(t1); + vpmb_start_gradient(); + /* For CVA calculations, start by guessing deco time = dive time remaining */ + deco_time = pi->maxtime - t1; + vpmb_next_gradient(deco_time, surface_pressure / 1000.0); + } + } + // Use the point where the ceiling clears as the end of deco phase for CVA calculations + if (entry->ceiling > 0) + time_clear_ceiling = 0; + else if (time_clear_ceiling == 0) + time_clear_ceiling = t1; + } + } + for (j = 0; j < 16; j++) { + double m_value = buehlmann_inertgas_a[j] + entry->ambpressure / buehlmann_inertgas_b[j]; + entry->ceilings[j] = deco_allowed_depth(tolerated_by_tissue[j], surface_pressure, dive, 1); + entry->percentages[j] = tissue_inertgas_saturation[j] < entry->ambpressure ? + tissue_inertgas_saturation[j] / entry->ambpressure * AMB_PERCENTAGE : + AMB_PERCENTAGE + (tissue_inertgas_saturation[j] - entry->ambpressure) / (m_value - entry->ambpressure) * (100.0 - AMB_PERCENTAGE); + } - /* should we do more calculations? - * We don't for print-mode because this info doesn't show up there */ - if (prefs.calcndltts && !print_mode) { - /* only calculate ndl/tts on every 30 seconds */ - if ((entry->sec - last_ndl_tts_calc_time) < 30) { - struct plot_data *prev_entry = (entry - 1); - entry->stoptime_calc = prev_entry->stoptime_calc; - entry->stopdepth_calc = prev_entry->stopdepth_calc; - entry->tts_calc = prev_entry->tts_calc; - entry->ndl_calc = prev_entry->ndl_calc; - continue; + /* should we do more calculations? + * We don't for print-mode because this info doesn't show up there + * If the ceiling hasn't cleared by the last data point, we need tts for VPM-B CVA calculation + * It is not necessary to do these calculation on the first VPMB iteration, except for the last data point */ + if ((prefs.calcndltts && !print_mode && (prefs.deco_mode != VPMB || in_planner() || !first_iteration)) || + (prefs.deco_mode == VPMB && !in_planner() && i == pi->nr - 1)) { + /* only calculate ndl/tts on every 30 seconds */ + if ((entry->sec - last_ndl_tts_calc_time) < 30 && i != pi->nr - 1) { + struct plot_data *prev_entry = (entry - 1); + entry->stoptime_calc = prev_entry->stoptime_calc; + entry->stopdepth_calc = prev_entry->stopdepth_calc; + entry->tts_calc = prev_entry->tts_calc; + entry->ndl_calc = prev_entry->ndl_calc; + continue; + } + last_ndl_tts_calc_time = entry->sec; + + /* We are going to mess up deco state, so store it for later restore */ + char *cache_data = NULL; + cache_deco_state(&cache_data); + calculate_ndl_tts(entry, dive, surface_pressure); + if (prefs.deco_mode == VPMB && !in_planner() && i == pi->nr - 1) + final_tts = entry->tts_calc; + /* Restore "real" deco state for next real time step */ + restore_deco_state(cache_data); + free(cache_data); } - last_ndl_tts_calc_time = entry->sec; - - /* We are going to mess up deco state, so store it for later restore */ - char *cache_data = NULL; - cache_deco_state(&cache_data); - calculate_ndl_tts(entry, dive, surface_pressure); - /* Restore "real" deco state for next real time step */ - restore_deco_state(cache_data); - free(cache_data); + } + if (prefs.deco_mode == VPMB && !in_planner()) { + prev_deco_time = deco_time; + // Do we need to update deco_time? + if (final_tts > 0) + deco_time = pi->maxtime + final_tts - time_deep_ceiling; + else if (time_clear_ceiling > 0) + deco_time = time_clear_ceiling - time_deep_ceiling; + vpmb_next_gradient(deco_time, surface_pressure / 1000.0); + final_tts = 0; + last_ndl_tts_calc_time = 0; + first_ceiling = 0; + first_iteration = false; + count_iteration ++; + restore_deco_state(cache_data_initial); + } else { + // With Buhlmann, or not in planner, iterating isn't needed. This makes the while condition false. + prev_deco_time = deco_time = 0; } } + free(cache_data_initial); #if DECO_CALC_DEBUG & 1 dump_tissues(); #endif -- 2.4.3
From 93cc7f73fa2671864a670cb1f7224cf17c7a25b9 Mon Sep 17 00:00:00 2001 From: Rick Walsh <[email protected]> Date: Sat, 24 Oct 2015 20:03:46 +1100 Subject: [PATCH 2/2] Profile: Display VPM-B rather than GF when in VPM-B mode If we are planning a dive using VPM-B with +x conservatism, we want to print "VPM-B +x" at the top of the profile, instead of "GF xx/yy". Accordingly, the variable gradientFactor in profilewidget2.cpp is renamed decoModelParameters to reflect what it represents. Signed-off-by: Rick Walsh <[email protected]> --- qt-ui/profile/profilewidget2.cpp | 32 +++++++++++++++++++------------- qt-ui/profile/profilewidget2.h | 2 +- 2 files changed, 20 insertions(+), 14 deletions(-) diff --git a/qt-ui/profile/profilewidget2.cpp b/qt-ui/profile/profilewidget2.cpp index 3ccd1bb..d474670 100644 --- a/qt-ui/profile/profilewidget2.cpp +++ b/qt-ui/profile/profilewidget2.cpp @@ -87,7 +87,7 @@ ProfileWidget2::ProfileWidget2(QWidget *parent) : QGraphicsView(parent), gasPressureItem(new DiveGasPressureItem()), diveComputerText(new DiveTextItem()), diveCeiling(new DiveCalculatedCeiling()), - gradientFactor(new DiveTextItem()), + decoModelParameters(new DiveTextItem()), reportedCeiling(new DiveReportedCeiling()), pn2GasItem(new PartialPressureGasItem()), pheGasItem(new PartialPressureGasItem()), @@ -204,7 +204,7 @@ void ProfileWidget2::addItemsToScene() diveComputerText->setData(SUBSURFACE_OBJ_DATA, SUBSURFACE_OBJ_DC_TEXT); scene()->addItem(diveComputerText); scene()->addItem(diveCeiling); - scene()->addItem(gradientFactor); + scene()->addItem(decoModelParameters); scene()->addItem(reportedCeiling); scene()->addItem(pn2GasItem); scene()->addItem(pheGasItem); @@ -288,11 +288,11 @@ void ProfileWidget2::setupItemOnScene() rulerItem->setAxis(timeAxis, profileYAxis); tankItem->setHorizontalAxis(timeAxis); - // show the gradient factor at the top in the center - gradientFactor->setY(0); - gradientFactor->setX(50); - gradientFactor->setBrush(getColor(PRESSURE_TEXT)); - gradientFactor->setAlignment(Qt::AlignHCenter | Qt::AlignBottom); + // show the deco model parameters at the top in the center + decoModelParameters->setY(0); + decoModelParameters->setX(50); + decoModelParameters->setBrush(getColor(PRESSURE_TEXT)); + decoModelParameters->setAlignment(Qt::AlignHCenter | Qt::AlignBottom); setupItem(reportedCeiling, timeAxis, profileYAxis, dataModel, DivePlotDataModel::CEILING, DivePlotDataModel::TIME, 1); setupItem(diveCeiling, timeAxis, profileYAxis, dataModel, DivePlotDataModel::CEILING, DivePlotDataModel::TIME, 1); @@ -509,7 +509,10 @@ void ProfileWidget2::plotDive(struct dive *d, bool force) // this copies the dive and makes copies of all the relevant additional data copy_dive(d, &displayed_dive); - gradientFactor->setText(QString("GF %1/%2").arg(prefs.gflow).arg(prefs.gfhigh)); + if (prefs.deco_mode == VPMB) + decoModelParameters->setText(QString("VPM-B +%1").arg(prefs.conservatism_level)); + else + decoModelParameters->setText(QString("GF %1/%2").arg(prefs.gflow).arg(prefs.gfhigh)); } else { DivePlannerPointsModel *plannerModel = DivePlannerPointsModel::instance(); plannerModel->createTemporaryPlan(); @@ -518,7 +521,10 @@ void ProfileWidget2::plotDive(struct dive *d, bool force) plannerModel->deleteTemporaryPlan(); return; } - gradientFactor->setText(QString("GF %1/%2").arg(diveplan.gflow).arg(diveplan.gfhigh)); + if (prefs.deco_mode == VPMB) + decoModelParameters->setText(QString("VPM-B +%1").arg(prefs.conservatism_level)); + else + decoModelParameters->setText(QString("GF %1/%2").arg(prefs.gflow).arg(prefs.gfhigh)); } // special handling for the first time we display things @@ -926,7 +932,7 @@ void ProfileWidget2::setEmptyState() toolTipItem->setVisible(false); diveComputerText->setVisible(false); diveCeiling->setVisible(false); - gradientFactor->setVisible(false); + decoModelParameters->setVisible(false); reportedCeiling->setVisible(false); rulerItem->setVisible(false); tankItem->setVisible(false); @@ -1053,7 +1059,7 @@ void ProfileWidget2::setProfileState() diveComputerText->setPos(itemPos.dcLabel.on); diveCeiling->setVisible(prefs.calcceiling); - gradientFactor->setVisible(prefs.calcceiling); + decoModelParameters->setVisible(prefs.calcceiling); reportedCeiling->setVisible(prefs.dcceiling); if (prefs.calcalltissues) { @@ -1131,7 +1137,7 @@ void ProfileWidget2::setAddState() /* show the same stuff that the profile shows. */ currentState = ADD; /* enable the add state. */ diveCeiling->setVisible(true); - gradientFactor->setVisible(true); + decoModelParameters->setVisible(true); setBackgroundBrush(QColor("#A7DCFF")); } @@ -1165,7 +1171,7 @@ void ProfileWidget2::setPlanState() /* show the same stuff that the profile shows. */ currentState = PLAN; /* enable the add state. */ diveCeiling->setVisible(true); - gradientFactor->setVisible(true); + decoModelParameters->setVisible(true); setBackgroundBrush(QColor("#D7E3EF")); } diff --git a/qt-ui/profile/profilewidget2.h b/qt-ui/profile/profilewidget2.h index 2d1a7bf..23e9398 100644 --- a/qt-ui/profile/profilewidget2.h +++ b/qt-ui/profile/profilewidget2.h @@ -171,7 +171,7 @@ private: QList<DiveEventItem *> eventItems; DiveTextItem *diveComputerText; DiveCalculatedCeiling *diveCeiling; - DiveTextItem *gradientFactor; + DiveTextItem *decoModelParameters; QList<DiveCalculatedTissue *> allTissues; DiveReportedCeiling *reportedCeiling; PartialPressureGasItem *pn2GasItem; -- 2.4.3
_______________________________________________ subsurface mailing list [email protected] http://lists.subsurface-divelog.org/cgi-bin/mailman/listinfo/subsurface
