view Discovery/Src/vpm.c @ 542:6960df7ddb09

Added new t3 view showing first deco stop and tts in parallel: TTS and deco stop may already be shown individually. Now also available in one view where TTS is limited to 99minutes (switch to hours) to avoid string overlap.
author Ideenmodellierer
date Sat, 10 Oct 2020 21:02:05 +0200
parents 305f251cc981
children b7d93ff6b3b2
line wrap: on
line source

///////////////////////////////////////////////////////////////////////////////
/// -*- coding: UTF-8 -*-
///
/// \file   Discovery/Src/vpm.c
/// \brief  critical_volume comment by hw
/// \author Heinrichs Weikamp, Erik C. Baker
/// \date   19-April-2014
///
/// \details
///
/// $Id$
///////////////////////////////////////////////////////////////////////////////
/// \par Copyright (c) 2014-2018 Heinrichs Weikamp gmbh
///
///     This program is free software: you can redistribute it and/or modify
///     it under the terms of the GNU General Public License as published by
///     the Free Software Foundation, either version 3 of the License, or
///     (at your option) any later version.
///
///     This program is distributed in the hope that it will be useful,
///     but WITHOUT ANY WARRANTY; without even the implied warranty of
///     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
///     GNU General Public License for more details.
///
///     You should have received a copy of the GNU General Public License
///     along with this program.  If not, see <http://www.gnu.org/licenses/>.
//////////////////////////////////////////////////////////////////////////////
/// \par Varying Permeability Model (VPM) Decompression Program in c (converted from FORTRAN)
///
///     Author:  Erik C. Baker
///
///     "DISTRIBUTE FREELY - CREDIT THE AUTHORS"
///
///     This program extends the 1986 VPM algorithm (Yount & Hoffman) to include
///     mixed gas, repetitive, and altitude diving.  Developments to the algorithm
///     were made by David E. Yount, Eric B. Maiken, and Erik C. Baker over a
///     period from 1999 to 2001.  This work is dedicated in remembrance of
///     Professor David E. Yount who passed away on April 27, 2000.
///
///     Notes:
///     1.  This program uses the sixteen (16) half-time compartments of the
///         Buhlmann ZH-L16 model.  The optional Compartment 1b is used here with
///         half-times of 1.88 minutes for helium and 5.0 minutes for nitrogen.
///
///     2.  This program uses various DEC, IBM, and Microsoft extensions which
///         may not be supported by all FORTRAN compilers.  Comments are made with
///         a capital "C" in the first column or an exclamation point "!" placed
///         in a line after code.  An asterisk "*" in column 6 is a continuation
///         of the previous line.  All code, except for line numbers, starts in
///         column 7.
///
///     3.  Comments and suggestions for improvements are welcome.  Please
///         respond by e-mail to:  EBaker@se.aeieng.com
///
///     Acknowledgment:  Thanks to Kurt Spaugh for recommendations on how to clean
///     up the code.
/// ===============================================================================
///     Converted to vpmdeco.c using f2c; R.McGinnis (CABER Swe) 5/01
/// ===============================================================================
///
/// ************************ Heirichs Weipkamp **************************************
///
/// The original Yount & Baker code has been adjusted for real life calculation.
///
/// 1) The original main function has been split in several functions
///
/// 2) When the deco zone is reached (while ascending) the gradient factors are kept fix
///     and critical volume algorithm is switched of. maxfirststopdepth is kept fix
///     to make shure Boeyls Law algorithm works correctly
///
/// 4) gas_loadings_ascent_descend heeds all gaschanges and CCR support has been added
///

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>

#include "vpm.h"
#include "decom.h"

#define GAS_N2 0
#define GAS_HE 1

static const _Bool buehlmannSafety = true;
/* Common Block Declarations */

extern const float SURFACE_TENSION_GAMMA;				//!Adj. Range: 0.015 to 0.065 N/m
extern const float SKIN_COMPRESSION_GAMMAC;				//!Adj. Range: 0.160 to 0.290 N/m
extern const float UNITS_FACTOR;
extern const float WATER_VAPOR_PRESSURE;				// (Schreiner value)  based on respiratory quotien
extern const float CRIT_VOLUME_PARAMETER_LAMBDA;			//!Adj. Range: 6500 to 8300 fsw-min
//extern const float GRADIENT_ONSET_OF_IMPERM_ATM;			//!Adj. Range: 5.0 to 10.0 atm
extern const float REGENERATION_TIME_CONSTANT;			//!Adj. Range: 10080 to 51840 min
//extern const float PRESSURE_OTHER_GASES_MMHG;				//!Constant value for PO2 up to 2 atm
extern const float CONSTANT_PRESSURE_OTHER_GASES; // PRESSURE_OTHER_GASES_MMHG / 760. * UNITS_FACTOR;

extern const float HELIUM_TIME_CONSTANT[];
extern const float NITROGEN_TIME_CONSTANT[];

static float    minimum_deco_stop_time;
static float    run_time, run_time_first_stop;
static float    segment_time;
static short mix_number;
static float    barometric_pressure;
static _Bool altitude_dive_algorithm_off;
static _Bool units_equal_fsw, units_equal_msw;

/* by hw 11.06.2015 to allow */
static float gCNS_VPM;

static float    helium_pressure[16], nitrogen_pressure[16];
static float    surface_phase_volume_time[16];
static float    regenerated_radius_he[16], regenerated_radius_n2[16];
static float    allowable_gradient_he[16], allowable_gradient_n2[16];

//_Bool deco_zone_reached;
static _Bool critical_volume_algorithm_off;
static float max_first_stop_depth;
static float max_deco_ceiling_depth;
//Boylslaw compensation
static float deco_gradient_he[16];
static float deco_gradient_n2[16];
static int vpm_calc_what;
static int count_critical_volume_iteration;
static short number_of_changes;
static float   depth_change[11];
static float  step_size_change[11];
static float    rate_change[11];
static short mix_change[11];

static const _Bool vpm_b = true;

extern const float float_buehlmann_N2_factor_expositon_20_seconds[];
extern const float float_buehlmann_He_factor_expositon_20_seconds[];
extern const float float_buehlmann_N2_factor_expositon_one_minute[];
extern const float float_buehlmann_He_factor_expositon_one_minute[];
extern const float float_buehlmann_N2_factor_expositon_five_minutes[];
extern const float float_buehlmann_He_factor_expositon_five_minutes[];
extern const float float_buehlmann_N2_factor_expositon_one_hour[];
extern const float float_buehlmann_He_factor_expositon_one_hour[];

static float    depth_start_of_deco_calc;
static float    depth_start_of_deco_zone;
static float    first_stop_depth;
static float    run_time_start_of_deco_zone;

static float r_nint(float *x);
static float r_int(float *x);
static _Bool nullzeit_unter60;
static int vpm_calc_status;
static _Bool buehlmann_wait_exceeded = false;

static SLifeData* pInput = NULL;
static SVpm* pVpm = NULL;
static SDecoinfo* pDecoInfo = NULL;
static SDiveSettings* pDiveSettings = NULL;

static float r_nint(float *x)
{
    return( (*x)>=0 ?
    floorf(*x + 0.5f) : -floorf(0.5f - *x) );
}

static float r_int(float *x)
{
    return( (*x>0) ? floorf(*x) : -floorf(- *x) );
}

/** private functions
*/
extern int radius_root_finder (float *a, float *b, float *c,float *low_bound,  float *high_bound, float *ending_radius);

static int nuclear_regeneration(float *dive_time);// clock_();
static int calc_deco_ceiling(float *deco_ceiling_depth,_Bool fallowablw);
static int critical_volume(float *deco_phase_volume_time);	    ;
static int calc_start_of_deco_zone(float *starting_depth, float *rate, float *depth_start_of_deco_zone);
static int calc_initial_allowable_gradient(void);
static void decompression_stop(float *deco_stop_depth, float *step_size, _Bool final_deco_calculation);
static int gas_loadings_ascent_descen(float* helium_pressure, float* nitrogen_pressure, float starting_depth,float ending_depth, float rate,_Bool check_gas_change);
static int calc_surface_phase_volume_time(void);
static int calc_max_actual_gradient(float *deco_stop_depth);
static int projected_ascent(float *starting_depth, float *rate, float *deco_stop_depth, float *step_size);
static void vpm_calc_deco(void);
static int vpm_calc_critcal_volume(_Bool begin,_Bool calc_nulltime);
static int vpm_check_converged(_Bool calc_nulltime);
static int vpm_calc_final_deco(_Bool begin);
static void BOYLES_LAW_COMPENSATION (float* First_Stop_Depth,float * Deco_Stop_Depth,float* Step_Size);
static int vpm_calc_ndl(void);
static void  vpm_init_1(void);
static void vpm_calc_deco_ceiling(void);

static void vpm_init_1(void)
{
     units_equal_msw = true;
     units_equal_fsw = false;
     altitude_dive_algorithm_off= true;			//!Options: ON or OFF
     minimum_deco_stop_time=1.0;					//!Options: float positive number
     critical_volume_algorithm_off= false;		//!Options: ON or OFF
     run_time = 0.;
    //barometric_pressure = dive_data.surface * 10;

    //mix_number = dive_data.selected_gas + 1;

    max_first_stop_depth = 0;
    max_deco_ceiling_depth = 0;
    //deco_zone_reached = false;
    depth_start_of_deco_calc = 0;
    depth_start_of_deco_zone = 0;
    first_stop_depth = 0;
    run_time_start_of_deco_zone = 0;

    gCNS_VPM = 0;
}

float vpm_get_CNS(void)
{
    return gCNS_VPM;
}

int  vpm_calc(SLifeData* pINPUT,
              SDiveSettings* pSettings,
              SVpm* pVPM,
              SDecoinfo*
              pDECOINFO,
              int calc_what)
{
    vpm_init_1();
    //decom_CreateGasChangeList(pSettings, pINPUT);
    vpm_calc_what = calc_what;
    /**clear decoInfo*/
    pDECOINFO->output_time_to_surface_seconds = 0;
    pDECOINFO->output_ndl_seconds = 0;
    pDECOINFO->output_ceiling_meter = 0;
    pDECOINFO->super_saturation = 0;
    uint8_t tmp_calc_status;
    for(int i=0;i<DECOINFO_STRUCT_MAX_STOPS;i++)
    {
        pDECOINFO->output_stop_length_seconds[i] = 0;
    }

    if(pINPUT->dive_time_seconds_without_surface_time < 60)
    {
        vpm_calc_status = CALC_NDL;
        return vpm_calc_status;
    }
    pVpm = pVPM;
    pInput = pINPUT;
    pDecoInfo = pDECOINFO;
    pDiveSettings = pSettings;

    if(vpm_calc_status == CALC_NDL)
    {
        tmp_calc_status = vpm_calc_ndl();
    }
    else
    {
        tmp_calc_status = CALC_BEGIN;
    }
    //Normal Deco calculation
    if(tmp_calc_status != CALC_NDL)
    {
        max_first_stop_depth =  pVpm->max_first_stop_depth_save;
        run_time_start_of_deco_zone = pVpm->run_time_start_of_deco_zone_save;
        depth_start_of_deco_zone = pVpm->depth_start_of_deco_zone_save;
        for (int i = 0; i < 16; ++i) {
            helium_pressure[i] = pInput->tissue_helium_bar[i] * 10;
            nitrogen_pressure[i] = pInput->tissue_nitrogen_bar[i] * 10;
        }
        vpm_calc_deco();
        tmp_calc_status = vpm_calc_critcal_volume(true,false);
        if(vpm_calc_what == DECOSTOPS)
        {
            pVpm->max_first_stop_depth_save = max_first_stop_depth;
            pVpm->run_time_start_of_deco_zone_save = run_time_start_of_deco_zone;
            pVpm->depth_start_of_deco_zone_save = depth_start_of_deco_zone;
        }
    }

    //Only Decostops not futute stops
    if(vpm_calc_what == DECOSTOPS)
        vpm_calc_status = tmp_calc_status;
    return vpm_calc_status;
}

void  vpm_saturation_after_ascent(SLifeData* input)
{
    int i = 0;
    for (i = 0; i < 16; ++i) {
        pInput->tissue_helium_bar[i] = helium_pressure[i] / 10;
        pInput->tissue_nitrogen_bar[i] = nitrogen_pressure[i] / 10;
    }
    pInput->pressure_ambient_bar = pInput->pressure_surface_bar;
}
/* =============================================================================== */
/*     NOTE ABOUT PRESSURE UNITS USED IN CALCULATIONS: */
/*     It is the convention in decompression calculations to compute all gas */
/*     loadings, absolute pressures, partial pressures, etc., in the units of */
/*     depth pressure that you are diving - either feet of seawater (fsw) or */
/*     meters of seawater (msw).  This program follows that convention with the */
/*     the exception that all VPM calculations are performed in SI units (by */
/*     necessity).  Accordingly, there are several conversions back and forth */
/*     between the diving pressure units and the SI units. */
/* =============================================================================== */
/* =============================================================================== */
/*     FUNCTION SUBPROGRAM FOR GAS LOADING CALCULATIONS - ASCENT AND DESCENT */
/* =============================================================================== */



/* =============================================================================== */
/*     SUBROUTINE GAS_LOADINGS_ASCENT_DESCENT */
/*     Purpose: This subprogram applies the Schreiner equation to update the */
/*     gas loadings (partial pressures of helium and nitrogen) in the half-time */
/*     compartments due to a linear ascent or descent segment at a constant rate. */
/* =============================================================================== */

static int gas_loadings_ascent_descen(float* helium_pressure,
                               float* nitrogen_pressure,
                               float starting_depth,
                               float ending_depth,
                               float rate,_Bool check_gas_change)
{
    short i;
    float initial_inspired_n2_pressure,
    initial_inspired_he_pressure, nitrogen_rate,
    last_run_time,
    starting_ambient_pressure,
    ending_ambient_pressure;
    float initial_helium_pressure[16];
    float initial_nitrogen_pressure[16];
    float helium_rate;
    float fraction_helium_begin;
    float fraction_helium_end;
    float fraction_nitrogen_begin;
    float fraction_nitrogen_end;
    float ending_depth_tmp = ending_depth;
    float segment_time_tmp = 0;
    /* loop */
    /* =============================================================================== */
    /*     CALCULATIONS */
    /* =============================================================================== */
    segment_time = (ending_depth_tmp - starting_depth) / rate;
    last_run_time = run_time;
    run_time = last_run_time + segment_time;
    do {
        ending_depth_tmp = ending_depth;
        if (starting_depth > ending_depth && check_gas_change && number_of_changes > 1)
        {
            for (i = 1; i < number_of_changes; ++i)
            {
                if (depth_change[i] < starting_depth && depth_change[i] > ending_depth)
                {
                    ending_depth_tmp = depth_change[i];
                    break;
                }
            }
            for (i = 1; i < number_of_changes; ++i)
            {
                if (depth_change[i] >= starting_depth)
                {
                    mix_number = mix_change[i];
                }
            }
        }
        segment_time_tmp = (ending_depth_tmp - starting_depth) / rate;
        ending_ambient_pressure = ending_depth_tmp + barometric_pressure;
        starting_ambient_pressure = starting_depth + barometric_pressure;
        decom_get_inert_gases( starting_ambient_pressure / 10,  (&pDiveSettings->decogaslist[mix_number]), &fraction_nitrogen_begin, &fraction_helium_begin );
        decom_get_inert_gases( ending_ambient_pressure   / 10,  (&pDiveSettings->decogaslist[mix_number]), &fraction_nitrogen_end, &fraction_helium_end );

        initial_inspired_he_pressure =	(starting_ambient_pressure - WATER_VAPOR_PRESSURE) * fraction_helium_begin;
        initial_inspired_n2_pressure =	(starting_ambient_pressure - WATER_VAPOR_PRESSURE) * fraction_nitrogen_begin;
        //helium_rate = *rate * fraction_helium[mix_number - 1];
        helium_rate = ((ending_ambient_pressure  - WATER_VAPOR_PRESSURE)* fraction_helium_end - initial_inspired_he_pressure)/segment_time_tmp;
        //nitrogen_rate2 = *rate * fraction_nitrogen[mix_number - 1];
        nitrogen_rate = ((ending_ambient_pressure  - WATER_VAPOR_PRESSURE)* fraction_nitrogen_end - initial_inspired_n2_pressure)/segment_time_tmp;


        decom_oxygen_calculate_cns_stage_SchreinerStyle(segment_time_tmp,&pDiveSettings->decogaslist[mix_number],starting_ambient_pressure/10,ending_ambient_pressure/10,&gCNS_VPM);
        //if(fabs(nitrogen_rate - nitrogen_rate2) > 0.000001)
            //return -2;
        for (i = 1; i <= 16; ++i)
        {
            initial_helium_pressure[i - 1] = helium_pressure[i - 1];
            initial_nitrogen_pressure[i - 1] = nitrogen_pressure[i - 1];
            helium_pressure[i - 1] =
            schreiner_equation__2(&initial_inspired_he_pressure,
            &helium_rate,
            &segment_time,
            &HELIUM_TIME_CONSTANT[i - 1],
            &initial_helium_pressure[i - 1]);
            nitrogen_pressure[i - 1] =
            schreiner_equation__2(&initial_inspired_n2_pressure,
            &nitrogen_rate,
            &segment_time,
            &NITROGEN_TIME_CONSTANT[i - 1],
            &initial_nitrogen_pressure[i - 1]);

            //nextround???

        }
        starting_depth = ending_depth_tmp;
    } while(ending_depth_tmp > ending_depth);

    return 0;
} /* gas_loadings_ascent_descen */

static float last_phase_volume_time[16];
static float n2_pressure_start_of_deco_zone[16];
static float he_pressure_start_of_deco_zone[16];
static float phase_volume_time[16];
static float n2_pressure_start_of_ascent[16];
static float  he_pressure_start_of_ascent[16];
static float    run_time_start_of_deco_calc;
static float starting_depth;
static float last_run_time;
static float deco_phase_volume_time;
static float run_time_start_of_ascent;
static float rate;
static float    step_size;
static _Bool vpm_violates_buehlmann;

static  void  vpm_calc_deco(void)
{
    /* System generated locals */

    //float   deepest_possible_stop_depth;
//	altitude_of_dive,
    short i;
    int j = 0;

    // float    rounding_operation;

    /* =============================================================================== */
    /*     INPUT PARAMETERS TO BE USED FOR STAGED DECOMPRESSION AND SAVE IN ARRAYS. */
    /*     ASSIGN INITAL PARAMETERS TO BE USED AT START OF ASCENT */
    /*     The user has the ability to change mix, ascent rate, and step size in any */
    /*     combination at any depth during the ascent. */
    /* =============================================================================== */

    run_time = ((float)pInput->dive_time_seconds )/ 60;
    count_critical_volume_iteration = 0;
    number_of_changes = 1;

    barometric_pressure = pInput->pressure_surface_bar * 10;
    depth_change[0] =(pInput->pressure_ambient_bar - pInput->pressure_surface_bar)* 10;
    mix_change[0] = 0;
    rate_change[0 ] = -10;// neu 160215 hw, zuvor: -12;
    step_size_change[0] = 3;
    vpm_violates_buehlmann = false;

    for (i = 1; i < BUEHLMANN_STRUCT_MAX_GASES; i++)
    {
            depth_change[i] = 0;
            mix_change[i] = 0;
    }
    j = 0;

    for (i = 1; i < BUEHLMANN_STRUCT_MAX_GASES; i++)
    {
        if(pDiveSettings->decogaslist[i].change_during_ascent_depth_meter_otherwise_zero  >= depth_change[0] + 1)
            continue;

        if(pDiveSettings->decogaslist[i].change_during_ascent_depth_meter_otherwise_zero <= 0)
            break;

        j++;
        number_of_changes ++;
        depth_change[j] = pDiveSettings->decogaslist[i].change_during_ascent_depth_meter_otherwise_zero ;
        mix_change[j] = i;
        rate_change[j] = -10;// neu 160215 hw, zuvor: -12;
        step_size_change[j] = 3;
    }

    starting_depth = depth_change[0] ;
    mix_number = mix_change[0] ;
    rate = rate_change[0];
    step_size = step_size_change[0];

    for (i = 0; i < 16; ++i) {
        he_pressure_start_of_ascent[i ] = 	helium_pressure[i];
        n2_pressure_start_of_ascent[i] = 	nitrogen_pressure[i];
    }
    run_time_start_of_ascent = run_time;
    if(starting_depth <= depth_start_of_deco_zone && vpm_calc_what == DECOSTOPS)
    {
        pVpm->deco_zone_reached = true;
        depth_start_of_deco_calc = starting_depth;
        critical_volume_algorithm_off = true;
    }
    else
    {
        //if(deco_zone_reached)
        //{
            pVpm->deco_zone_reached = false;
            critical_volume_algorithm_off = false;
            //max_first_stop_depth = 0;
            //max_first_stop_depth_save = 0;
        //}
        /* =============================================================================== */
        /*     BEGIN PROCESS OF ASCENT AND DECOMPRESSION */
        /*     First, calculate the regeneration of critical radii that takes place over */
        /*     the dive time.  The regeneration time constant has a time scale of weeks */
        /*     so this will have very little impact on dives of normal length, but will */
        /*     have major impact for saturation dives. */
        /* =============================================================================== */

        nuclear_regeneration(&run_time);

        /* =============================================================================== */
        /*     CALCULATE INITIAL ALLOWABLE GRADIENTS FOR ASCENT */
        /*     This is based on the maximum effective crushing pressure on critical radii */
        /*     in each compartment achieved during the dive profile. */
        /* =============================================================================== */

        calc_initial_allowable_gradient();

        /* =============================================================================== */
        /*     SAVE VARIABLES AT START OF ASCENT (END OF BOTTOM TIME) SINCE THESE WILL */
        /*     BE USED LATER TO COMPUTE THE FINAL ASCENT PROFILE THAT IS WRITTEN TO THE */
        /*     OUTPUT FILE. */
        /*     The VPM uses an iterative process to compute decompression schedules so */
        /*     there will be more than one pass through the decompression loop. */
        /* =============================================================================== */

        /* =============================================================================== */
        /*     CALCULATE THE DEPTH WHERE THE DECOMPRESSION ZONE BEGINS FOR THIS PROFILE */
        /*     BASED ON THE INITIAL ASCENT PARAMETERS AND WRITE THE DEEPEST POSSIBLE */
        /*     DECOMPRESSION STOP DEPTH TO THE OUTPUT FILE */
        /*     Knowing where the decompression zone starts is very important.  Below */
        /*     that depth there is no possibility for bubble formation because there */
        /*     will be no supersaturation gradients.  Deco stops should never start */
        /*     below the deco zone.  The deepest possible stop deco stop depth is */
        /*     defined as the next "standard" stop depth above the point where the */
        /*     leading compartment enters the deco zone.  Thus, the program will not */
        /*     base this calculation on step sizes larger than 10 fsw or 3 msw.  The */
        /*     deepest possible stop depth is not used in the program, per se, rather */
        /*     it is information to tell the diver where to start putting on the brakes */
        /*     during ascent.  This should be prominently displayed by any deco program. */
        /* =============================================================================== */

        calc_start_of_deco_zone(&starting_depth, &rate, &depth_start_of_deco_zone);
        /* =============================================================================== */
        /*     TEMPORARILY ASCEND PROFILE TO THE START OF THE DECOMPRESSION ZONE, SAVE */
        /*     VARIABLES AT THIS POINT, AND INITIALIZE VARIABLES FOR CRITICAL VOLUME LOOP */
        /*     The iterative process of the VPM Critical Volume Algorithm will operate */
        /*     only in the decompression zone since it deals with excess gas volume */
        /*     released as a result of supersaturation gradients (not possible below the */
        /*     decompression zone). */
        /* =============================================================================== */
        gas_loadings_ascent_descen(helium_pressure,nitrogen_pressure, starting_depth, depth_start_of_deco_zone, rate, true);

        run_time_start_of_deco_zone = run_time;
        depth_start_of_deco_calc = depth_start_of_deco_zone;

        for (i = 0; i < 16; ++i)
        {
            pVpm->max_actual_gradient[i] = 0.;
        }
    }

    for (i = 0; i < 16; ++i)
    {
        surface_phase_volume_time[i] = 0.;
        last_phase_volume_time[i] = 0.;
        he_pressure_start_of_deco_zone[i] =	helium_pressure[i];
        n2_pressure_start_of_deco_zone[i] =	nitrogen_pressure[i];
        //pVpm->max_actual_gradient[i] = 0.;
    }
    run_time_start_of_deco_calc = run_time;
}
/* =============================================================================== */
/*     START OF CRITICAL VOLUME LOOP */
/*     This loop operates between Lines 50 and 100.  If the Critical Volume */
/*     Algorithm is toggled "off" in the program settings, there will only be */
/*     one pass through this loop.  Otherwise, there will be two or more passes */
/*     through this loop until the deco schedule is "converged" - that is when a */
/*     comparison between the phase volume time of the present iteration and the */
/*     last iteration is less than or equal to one minute.  This implies that */
/*     the volume of released gas in the most recent iteration differs from the */
/*     "critical" volume limit by an acceptably small amount.  The critical */
/*     volume limit is set by the Critical Volume Parameter Lambda in the program */
/*     settings (default setting is 7500 fsw-min with adjustability range from */
/*     from 6500 to 8300 fsw-min according to Bruce Wienke). */
/* =============================================================================== */
/* L50: */

static float  deco_stop_depth;
static int vpm_calc_critcal_volume(_Bool begin,
                            _Bool calc_nulltime)
{   /* loop will run continuous there is an exit stateme */

    short i;

    float rounding_operation2;
    //float    ending_depth;
    float deco_ceiling_depth;

    //float deco_time;
    int count = 0;
    _Bool first_stop;
    int dp = 0;
    float tissue_He_saturation[16];
    float tissue_N2_saturation[16];
    float vpm_buehlmann_safety_gradient = 1.0f - (((float)pDiveSettings->vpm_conservatism) / 40);
    /* =============================================================================== */
    /*     CALCULATE CURRENT DECO CEILING BASED ON ALLOWABLE SUPERSATURATION */
    /*     GRADIENTS AND SET FIRST DECO STOP.  CHECK TO MAKE SURE THAT SELECTED STEP */
    /*     SIZE WILL NOT ROUND UP FIRST STOP TO A DEPTH THAT IS BELOW THE DECO ZONE. */
    /* =============================================================================== */
    if(begin)
    {
        if(depth_start_of_deco_calc < max_first_stop_depth )
        {
            if(vpm_b)
            {
                BOYLES_LAW_COMPENSATION(&max_first_stop_depth, &depth_start_of_deco_calc, &step_size);
            }
            calc_deco_ceiling(&deco_ceiling_depth, false);
        }
        else
            calc_deco_ceiling(&deco_ceiling_depth, true);


        if (deco_ceiling_depth <= 0.0f) {
            deco_stop_depth = 0.0f;
        } else {
            rounding_operation2 = deco_ceiling_depth / step_size + ( float)0.5f;
            deco_stop_depth = r_nint(&rounding_operation2) * step_size;
        }

        // buehlmann safety
        if(buehlmannSafety)
        {
            for (i = 0; i < 16; i++)
            {
                tissue_He_saturation[i] = helium_pressure[i] / 10;
                tissue_N2_saturation[i] = nitrogen_pressure[i] / 10;
            }

            if(!decom_tissue_test_tolerance(tissue_N2_saturation, tissue_He_saturation, vpm_buehlmann_safety_gradient, (deco_stop_depth / 10.0f) + pInput->pressure_surface_bar))
            {

                vpm_violates_buehlmann = true;
                do {
                    deco_stop_depth += 3;
                } while (!decom_tissue_test_tolerance(tissue_N2_saturation, tissue_He_saturation, vpm_buehlmann_safety_gradient, (deco_stop_depth / 10.0f) + pInput->pressure_surface_bar));
            }
        }

        /* =============================================================================== */
        /*     PERFORM A SEPARATE "PROJECTED ASCENT" OUTSIDE OF THE MAIN PROGRAM TO MAKE */
        /*     SURE THAT AN INCREASE IN GAS LOADINGS DURING ASCENT TO THE FIRST STOP WILL */
        /*     NOT CAUSE A VIOLATION OF THE DECO CEILING.  IF SO, ADJUST THE FIRST STOP */
        /*     DEEPER BASED ON STEP SIZE UNTIL A SAFE ASCENT CAN BE MADE. */
        /*     Note: this situation is a possibility when ascending from extremely deep */
        /*     dives or due to an unusual gas mix selection. */
        /*     CHECK AGAIN TO MAKE SURE THAT ADJUSTED FIRST STOP WILL NOT BE BELOW THE */
        /*     DECO ZONE. */
        /* =============================================================================== */
        if (deco_stop_depth < depth_start_of_deco_calc)
        {
            projected_ascent(&depth_start_of_deco_calc, &rate, &deco_stop_depth, &step_size);
        }

        /*if (deco_stop_depth > depth_start_of_deco_zone) {
            printf("\t\n");
            printf(fmt_905);
            printf(fmt_900);
            printf("\nPROGRAM TERMINATED\n");
            exit(1);
        }*/

        /* =============================================================================== */
        /*     HANDLE THE SPECIAL CASE WHEN NO DECO STOPS ARE REQUIRED - ASCENT CAN BE */
        /*     MADE DIRECTLY TO THE SURFACE */
        /*     Write ascent data to output file and exit the Critical Volume Loop. */
        /* =============================================================================== */

        if (deco_stop_depth == 0.0f)
        {
            if(calc_nulltime)
            {
                return CALC_END;
            }
            if(pVpm->deco_zone_reached)
            {
                for(dp = 0;dp < DECOINFO_STRUCT_MAX_STOPS;dp++)
                {
                    pDecoInfo->output_stop_length_seconds[dp] = 0;
                }
                pDecoInfo->output_ndl_seconds = 0;
            }

            return CALC_NDL;
                    /* exit the critical volume l */
        }

    /* =============================================================================== */
    /*     ASSIGN VARIABLES FOR ASCENT FROM START OF DECO ZONE TO FIRST STOP.  SAVE */
    /*     FIRST STOP DEPTH FOR LATER USE WHEN COMPUTING THE FINAL ASCENT PROFILE */
    /* =============================================================================== */
        deco_stop_depth = fmaxf(deco_stop_depth,(float)pDiveSettings->last_stop_depth_bar * 10);
        starting_depth = depth_start_of_deco_calc;
        first_stop_depth = deco_stop_depth;
        first_stop = true;
    }
    /* =============================================================================== */
    /*     DECO STOP LOOP BLOCK WITHIN CRITICAL VOLUME LOOP */
    /*     This loop computes a decompression schedule to the surface during each */
    /*     iteration of the critical volume loop.  No output is written from this */
    /*     loop, rather it computes a schedule from which the in-water portion of the */
    /*     total phase volume time (Deco_Phase_Volume_Time) can be extracted.  Also, */
    /*     the gas loadings computed at the end of this loop are used the subroutine */
    /*     which computes the out-of-water portion of the total phase volume time */
    /*     (Surface_Phase_Volume_Time) for that schedule. */

    /*     Note that exit is made from the loop after last ascent is made to a deco */
    /*     stop depth that is less than or equal to zero.  A final deco stop less */
    /*     than zero can happen when the user makes an odd step size change during */
    /*     ascent - such as specifying a 5 msw step size change at the 3 msw stop! */
    /* =============================================================================== */

    while(true)           /* loop will run continuous there is an break statement */
    {
        if(starting_depth > deco_stop_depth )
            gas_loadings_ascent_descen(helium_pressure, nitrogen_pressure, starting_depth, deco_stop_depth, rate,first_stop);

        first_stop = false;
        if (deco_stop_depth <= 0.0f)
        {
            break;
        }
        if (number_of_changes > 1)
        {
            int i1 = number_of_changes;
            for (i = 2; i <= i1; ++i) {
                if (depth_change[i - 1] >= deco_stop_depth)
                {
                    mix_number = mix_change[i - 1];
                    rate = rate_change[i - 1];
                    step_size = step_size_change[i - 1];
                }
            }
        }
        if(vpm_b)
        {
            float fist_stop_depth2 =  fmaxf(first_stop_depth,max_first_stop_depth);
            BOYLES_LAW_COMPENSATION(&fist_stop_depth2, &deco_stop_depth, &step_size);
        }
        decompression_stop(&deco_stop_depth, &step_size, false);
        starting_depth = deco_stop_depth;

        if(deco_stop_depth == (float)pDiveSettings->last_stop_depth_bar * 10)
            deco_stop_depth = 0;
        else
        {
            deco_stop_depth = deco_stop_depth - step_size;
            deco_stop_depth = fmaxf(deco_stop_depth,(float)pDiveSettings->last_stop_depth_bar * 10);
        }

        count++;
        //if(count > 14)
            //return CALC_CRITICAL2;
        /* L60: */
    }

        return vpm_check_converged(calc_nulltime);
}
/* =============================================================================== */
/*     COMPUTE TOTAL PHASE VOLUME TIME AND MAKE CRITICAL VOLUME COMPARISON */
/*     The deco phase volume time is computed from the run time.  The surface */
/*     phase volume time is computed in a subroutine based on the surfacing gas */
/*     loadings from previous deco loop block.  Next the total phase volume time */
/*     (in-water + surface) for each compartment is compared against the previous */
/*     total phase volume time.  The schedule is converged when the difference is */
/*     less than or equal to 1 minute in any one of the 16 compartments. */

/*     Note:  the "phase volume time" is somewhat of a mathematical concept. */
/*     It is the time divided out of a total integration of supersaturation */
/*     gradient x time (in-water and surface).  This integration is multiplied */
/*     by the excess bubble number to represent the amount of free-gas released */
/*     as a result of allowing a certain number of excess bubbles to form. */
/* =============================================================================== */
/* end of deco stop loop */

static int vpm_check_converged(_Bool calc_nulltime)
{

    short i;
    float    critical_volume_comparison;
    float    r1;
    _Bool schedule_converged = false;


    deco_phase_volume_time = run_time - run_time_start_of_deco_zone;
    calc_surface_phase_volume_time();

    for (i = 1; i <= 16; ++i)
    {
        phase_volume_time[i - 1] =
        deco_phase_volume_time + surface_phase_volume_time[i - 1];
        critical_volume_comparison = (r1 = phase_volume_time[i - 1] - last_phase_volume_time[i - 1], fabs(r1));

        if (critical_volume_comparison <= 1.0f)
        {
            schedule_converged = true;
        }
    }

/* =============================================================================== */
/*     CRITICAL VOLUME DECISION TREE BETWEEN LINES 70 AND 99 */
/*     There are two options here.  If the Critical Volume Agorithm setting is */
/*     "on" and the schedule is converged, or the Critical Volume Algorithm */
/*     setting was "off" in the first place, the program will re-assign variables */
/*     to their values at the start of ascent (end of bottom time) and process */
/*     a complete decompression schedule once again using all the same ascent */
/*     parameters and first stop depth.  This decompression schedule will match */
/*     the last iteration of the Critical Volume Loop and the program will write */
/*     the final deco schedule to the output file. */

/*     Note: if the Critical Volume Agorithm setting was "off", the final deco */
/*     schedule will be based on "Initial Allowable Supersaturation Gradients." */
/*     If it was "on", the final schedule will be based on "Adjusted Allowable */
/*     Supersaturation Gradients" (gradients that are "relaxed" as a result of */
/*     the Critical Volume Algorithm). */

/*     If the Critical Volume Agorithm setting is "on" and the schedule is not */
/*     converged, the program will re-assign variables to their values at the */
/*     start of the deco zone and process another trial decompression schedule. */
/* =============================================================================== */
/* L70: */
    //Not more than 4 iteration allowed
    count_critical_volume_iteration++;
    if(count_critical_volume_iteration > 4)
    {
        //return CALC_FINAL_DECO;
        if(calc_nulltime)
            return CALC_FINAL_DECO;
        else
            return vpm_calc_final_deco(true);
    }
    if (schedule_converged || critical_volume_algorithm_off)
    {

        //return CALC_FINAL_DECO;
        if(calc_nulltime)
            return CALC_FINAL_DECO;
        else
            return vpm_calc_final_deco(true);
/* final deco schedule */
/* exit critical volume l */

/* =============================================================================== */
/*     IF SCHEDULE NOT CONVERGED, COMPUTE RELAXED ALLOWABLE SUPERSATURATION */
/*     GRADIENTS WITH VPM CRITICAL VOLUME ALGORITHM AND PROCESS ANOTHER */
/*     ITERATION OF THE CRITICAL VOLUME LOOP */
/* =============================================================================== */

    } else {
        critical_volume(&deco_phase_volume_time);
        deco_phase_volume_time = 0.;
        run_time = run_time_start_of_deco_calc;
        starting_depth = depth_start_of_deco_calc;
        mix_number = mix_change[0];
        rate = rate_change[0];
        step_size = step_size_change[0];
        for (i = 1; i <= 16; ++i)
        {
                last_phase_volume_time[i - 1] = phase_volume_time[i - 1];
                helium_pressure[i - 1] = he_pressure_start_of_deco_zone[i - 1];
                nitrogen_pressure[i - 1] = n2_pressure_start_of_deco_zone[i - 1];
        }
        if(calc_nulltime)
            return CALC_CRITICAL;
        else
            return vpm_calc_critcal_volume(true, false);
    }
/* end of critical volume decision */
/* L100: */
    //  }/* end of critical vol loop */
}

static void vpm_calc_deco_ceiling(void)
{

    short i;
// hw 1601209	 float    r1;
// hw 1601209	 float    stop_time;
// hw 1601209	 int count = 0;
    //static int dp_max;
    //static float surfacetime;
    // _Bool first_stop = false;
    float tissue_He_saturation[16];
    float tissue_N2_saturation[16];
    float vpm_buehlmann_safety_gradient = 1.0f - (((float)pDiveSettings->vpm_conservatism) / 40);
    //max_first_stop_depth = fmaxf(first_stop_depth,max_first_stop_depth);

    /** CALC DECO Ceiling ******************************************************************/
    /** Not when Future stops */
    if(vpm_calc_what == DECOSTOPS)
    {

        for (i = 1; i <= 16; ++i)
        {
            helium_pressure[i - 1] = he_pressure_start_of_deco_zone[i - 1];
            nitrogen_pressure[i - 1] = n2_pressure_start_of_deco_zone[i - 1];
        }
        run_time = run_time_start_of_ascent;// run_time_start_of_ascent;
        starting_depth = depth_change[0];
        mix_number = mix_change[0];
        rate = rate_change[0];
        //gas_loadings_ascent_descen(helium_pressure,nitrogen_pressure, starting_depth, depth_start_of_deco_calc, rate, true);

        float deco_ceiling_depth = 0.0f;
        if(depth_start_of_deco_calc > max_deco_ceiling_depth)
        {
            calc_deco_ceiling(&deco_ceiling_depth, true);
        }
        if(buehlmannSafety)
        {
            for (i = 0; i < 16; i++)
            {
                tissue_He_saturation[i] = helium_pressure[i] / 10;
                tissue_N2_saturation[i] = nitrogen_pressure[i] / 10;
            }

            if(!decom_tissue_test_tolerance(tissue_N2_saturation, tissue_He_saturation, vpm_buehlmann_safety_gradient, (deco_ceiling_depth / 10.0f) + pInput->pressure_surface_bar))
            {

                vpm_violates_buehlmann = true;
                do {
                    deco_ceiling_depth += 0.1f;
                } while (!decom_tissue_test_tolerance(tissue_N2_saturation, tissue_He_saturation, vpm_buehlmann_safety_gradient, (deco_ceiling_depth / 10.0f) + pInput->pressure_surface_bar));
            }
        }

        if (deco_ceiling_depth < depth_start_of_deco_calc)
        {
            projected_ascent(&depth_start_of_deco_calc, &rate, &deco_ceiling_depth, &step_size);
        }

        max_deco_ceiling_depth = fmaxf(max_deco_ceiling_depth,deco_ceiling_depth);

        if(depth_start_of_deco_calc > deco_ceiling_depth)
        {
            gas_loadings_ascent_descen(helium_pressure,nitrogen_pressure, depth_start_of_deco_calc,deco_ceiling_depth, rate, true);
            //surfacetime += segment_time;
        }

        if(vpm_b)
        {
            BOYLES_LAW_COMPENSATION(&max_deco_ceiling_depth, &deco_ceiling_depth, &step_size);
        }
        calc_deco_ceiling(&deco_ceiling_depth, false);

        // buehlmann safety
        if(vpm_violates_buehlmann)
        {
            for (i = 0; i < 16; i++)
            {
                tissue_He_saturation[i] = helium_pressure[i] / 10;
                tissue_N2_saturation[i] = nitrogen_pressure[i] / 10;
            }

            if(!decom_tissue_test_tolerance(tissue_N2_saturation, tissue_He_saturation, vpm_buehlmann_safety_gradient, (deco_ceiling_depth / 10.0f) + pInput->pressure_surface_bar))
            {

                vpm_violates_buehlmann = true;
                do {
                    deco_ceiling_depth += 0.1f;
                } while (!decom_tissue_test_tolerance(tissue_N2_saturation, tissue_He_saturation, vpm_buehlmann_safety_gradient, (deco_ceiling_depth / 10.0f) + pInput->pressure_surface_bar));
            }
        }
        // output_ceiling_meter
        if(deco_ceiling_depth > first_stop_depth)
            deco_ceiling_depth = first_stop_depth;
        pDecoInfo->output_ceiling_meter = deco_ceiling_depth ;
    }
    else
    {
         pDecoInfo->output_ceiling_meter = 0;
    }

    // fix hw 160627
    if(pDecoInfo->output_ceiling_meter < 0)
        pDecoInfo->output_ceiling_meter = 0;

    /*** End CALC ceiling    ***************************************************/
}


/* =============================================================================== */
/*     DECO STOP LOOP BLOCK FOR FINAL DECOMPRESSION SCHEDULE */
/* =============================================================================== */

static int vpm_calc_final_deco(_Bool begin)
{
    short i;
    float    r1;
    float    stop_time;
    int count = 0;
    static int dp_max;
    static float surfacetime;
    _Bool first_stop = false;
    max_first_stop_depth = fmaxf(first_stop_depth,max_first_stop_depth);
    if(begin)
    {
        gCNS_VPM = 0;
        dp_max = 0;
        for (i = 1; i <= 16; ++i)
        {
            helium_pressure[i - 1] =
            he_pressure_start_of_ascent[i - 1];
            nitrogen_pressure[i - 1] =
            n2_pressure_start_of_ascent[i - 1];
        }
        run_time = run_time_start_of_ascent;// run_time_start_of_ascent;
        starting_depth = depth_change[0];
        mix_number = mix_change[0];
        rate = rate_change[0];
        step_size = step_size_change[0];
        deco_stop_depth = first_stop_depth;
        max_first_stop_depth = fmaxf(first_stop_depth,max_first_stop_depth);
        last_run_time = 0.;



        /* =============================================================================== */
        /*     DECO STOP LOOP BLOCK FOR FINAL DECOMPRESSION SCHEDULE */
        /* =============================================================================== */
        surfacetime = 0;
        first_stop = true;
    }

    while(true)       /* loop will run continuous until there is an break statement */
    {
        if(starting_depth > deco_stop_depth)
        {
            gas_loadings_ascent_descen(helium_pressure,nitrogen_pressure, starting_depth,deco_stop_depth, rate, first_stop);
            surfacetime += segment_time;
        }

        /* =============================================================================== */
        /*     DURING FINAL DECOMPRESSION SCHEDULE PROCESS, COMPUTE MAXIMUM ACTUAL */
        /*     SUPERSATURATION GRADIENT RESULTING IN EACH COMPARTMENT */
        /*     If there is a repetitive dive, this will be used later in the VPM */
        /*     Repetitive Algorithm to adjust the values for critical radii. */
        /* =============================================================================== */
        if(vpm_calc_what == DECOSTOPS)
            calc_max_actual_gradient(&deco_stop_depth);

        if (deco_stop_depth <= 0.0f) {
            break;
        }
        if (number_of_changes > 1)
        {
            int i1 = number_of_changes;
            for (i = 2; i <= i1; ++i)
            {
                if (depth_change[i - 1] >= deco_stop_depth)
                {
                    mix_number = mix_change[i - 1];
                    rate = rate_change[i - 1];
                    step_size = step_size_change[i - 1];
                }
            }
        }

        if(first_stop)
        {
            run_time_first_stop = run_time;
            first_stop = false;
        }
        if(vpm_b)
        {
            BOYLES_LAW_COMPENSATION(&max_first_stop_depth, &deco_stop_depth, &step_size);
        }
        decompression_stop(&deco_stop_depth, &step_size, true);

        /* =============================================================================== */
        /*     This next bit justs rounds up the stop time at the first stop to be in */
        /*     whole increments of the minimum stop time (to make for a nice deco table). */
        /* =============================================================================== */

        if (last_run_time == 0.0f)
        {
            r1 = segment_time / minimum_deco_stop_time + 0.5f;
            stop_time = r_int(&r1) * minimum_deco_stop_time;
        } else {
            stop_time = run_time - last_run_time;
        }
        stop_time = segment_time;
        surfacetime += stop_time;
        if((vpm_calc_what == DECOSTOPS) || (vpm_calc_what == BAILOUTSTOPS))
        {
            int dp = 0;
            if(deco_stop_depth == (float)pDiveSettings->last_stop_depth_bar * 10)
            {
                dp = 0;
            }
            else
            {
                    dp = 1 + (int)((deco_stop_depth - (pDiveSettings->input_second_to_last_stop_depth_bar * 10)) / step_size);
            }
            dp_max = (int)fmaxf(dp_max,dp);
            if(dp < DECOINFO_STRUCT_MAX_STOPS)
            {
                int stop_time_seconds = fminf((999 * 60), (int)(stop_time *60));
                //

                //if(vpm_calc_what == DECOSTOPS)
                    pDecoInfo->output_stop_length_seconds[dp] = (unsigned short)stop_time_seconds;
                //else
                    //decostop_bailout[dp] = (unsigned short)stop_time_seconds;
            }
        }


        /* =============================================================================== */
        /*     DURING FINAL DECOMPRESSION SCHEDULE, IF MINIMUM STOP TIME PARAMETER IS A */
        /*     WHOLE NUMBER (i.e. 1 minute) THEN WRITE DECO SCHEDULE USING short */
        /*     NUMBERS (looks nicer).  OTHERWISE, USE DECIMAL NUMBERS. */
        /*     Note: per the request of a noted exploration diver(!), program now allows */
        /*     a minimum stop time of less than one minute so that total ascent time can */
        /*     be minimized on very long dives.  In fact, with step size set at 1 fsw or */
        /*     0.2 msw and minimum stop time set at 0.1 minute (6 seconds), a near */
        /*     continuous decompression schedule can be computed. */
        /* =============================================================================== */

        starting_depth = deco_stop_depth;
        if(deco_stop_depth == (float)pDiveSettings->last_stop_depth_bar * 10)
            deco_stop_depth = 0;
        else
        {
            deco_stop_depth = deco_stop_depth - step_size;
            deco_stop_depth = fmaxf(deco_stop_depth,(float)pDiveSettings->last_stop_depth_bar * 10);
        }

        last_run_time = run_time;
        count++;
        //if(count > 14)
            //return CALC_FINAL_DECO2;
        /* L80: */
    } /* for final deco sche */

    if( (vpm_calc_what == DECOSTOPS)  || (vpm_calc_what == BAILOUTSTOPS))
    {
        for(int dp = dp_max +1;dp < DECOINFO_STRUCT_MAX_STOPS;dp++)
        {
            //if(vpm_calc_what == DECOSTOPS)
                pDecoInfo->output_stop_length_seconds[dp] = 0;
            //else
                //decostop_bailout[dp] = 0;
        }
    }
    pDecoInfo->output_time_to_surface_seconds = (int)(surfacetime * 60);
    pDecoInfo->output_ndl_seconds = 0;

    vpm_calc_deco_ceiling();
    /* end of deco stop lo */
    return CALC_END;
}

/* =============================================================================== */
/*     SUBROUTINE NUCLEAR_REGENERATION */
/*     Purpose: This subprogram calculates the regeneration of VPM critical */
/*     radii that takes place over the dive time.  The regeneration time constant */
/*     has a time scale of weeks so this will have very little impact on dives of */
/*     normal length, but will have a major impact for saturation dives. */
/* =============================================================================== */

static int nuclear_regeneration(float *dive_time)
{
    /* Local variables */
    float crush_pressure_adjust_ratio_he,
    ending_radius_n2,
    ending_radius_he;
    short i;
    float crushing_pressure_pascals_n2,
    crushing_pressure_pascals_he,
    adj_crush_pressure_n2_pascals,
    adj_crush_pressure_he_pascals,
    crush_pressure_adjust_ratio_n2;

    /* loop */
    /* =============================================================================== */
    /*     CALCULATIONS */
    /*     First convert the maximum crushing pressure obtained for each compartment */
    /*     to Pascals.  Next, compute the ending radius for helium and nitrogen */
    /*     critical nuclei in each compartment. */
    /* =============================================================================== */

    for (i = 1; i <= 16; ++i)
    {
        crushing_pressure_pascals_he =
        pVpm->max_crushing_pressure_he[i - 1] / UNITS_FACTOR * 101325.0f;
        crushing_pressure_pascals_n2 =
        pVpm->max_crushing_pressure_n2[i - 1] / UNITS_FACTOR * 101325.0f;
        ending_radius_he =
        1.0f / (crushing_pressure_pascals_he /
        ((SKIN_COMPRESSION_GAMMAC - SURFACE_TENSION_GAMMA) * 2.0f) +
        1.0f / pVpm->adjusted_critical_radius_he[i - 1]);
        ending_radius_n2 =
        1.0f / (crushing_pressure_pascals_n2 /
        ((SKIN_COMPRESSION_GAMMAC - SURFACE_TENSION_GAMMA) * 2.0f) +
        1.0f / pVpm->adjusted_critical_radius_n2[i - 1]);

        /* =============================================================================== */
        /*     A "regenerated" radius for each nucleus is now calculated based on the */
        /*     regeneration time constant.  This means that after application of */
        /*     crushing pressure and reduction in radius, a nucleus will slowly grow */
        /*     back to its original initial radius over a period of time.  This */
        /*     phenomenon is probabilistic in nature and depends on absolute temperature. */
        /*     It is independent of crushing pressure. */
        /* =============================================================================== */

        regenerated_radius_he[i - 1] =
        pVpm->adjusted_critical_radius_he[i - 1] +
        (ending_radius_he - pVpm->adjusted_critical_radius_he[i - 1]) *
        expf(-(*dive_time) / REGENERATION_TIME_CONSTANT);
        regenerated_radius_n2[i - 1] =
        pVpm->adjusted_critical_radius_n2[i - 1] +
        (ending_radius_n2 - pVpm->adjusted_critical_radius_n2[i - 1]) *
        expf(-(*dive_time) / REGENERATION_TIME_CONSTANT);

        /* =============================================================================== */
        /*     In order to preserve reference back to the initial critical radii after */
        /*     regeneration, an "adjusted crushing pressure" for the nuclei in each */
        /*     compartment must be computed.  In other words, this is the value of */
        /*     crushing pressure that would have reduced the original nucleus to the */
        /*     to the present radius had regeneration not taken place.  The ratio */
        /*     for adjusting crushing pressure is obtained from algebraic manipulation */
        /*     of the standard VPM equations.  The adjusted crushing pressure, in lieu */
        /*     of the original crushing pressure, is then applied in the VPM Critical */
        /*     Volume Algorithm and the VPM Repetitive Algorithm. */
        /* =============================================================================== */

        crush_pressure_adjust_ratio_he =
        ending_radius_he * (pVpm->adjusted_critical_radius_he[i - 1] -
        regenerated_radius_he[i - 1]) /
        (regenerated_radius_he[i - 1] *
        (pVpm->adjusted_critical_radius_he[i - 1] -
        ending_radius_he));
        crush_pressure_adjust_ratio_n2 =
        ending_radius_n2 * (pVpm->adjusted_critical_radius_n2[i - 1] -
        regenerated_radius_n2[i - 1]) /
        (regenerated_radius_n2[i - 1] *
        (pVpm->adjusted_critical_radius_n2[i - 1] -
        ending_radius_n2));
        adj_crush_pressure_he_pascals =
        crushing_pressure_pascals_he * crush_pressure_adjust_ratio_he;
        adj_crush_pressure_n2_pascals =
        crushing_pressure_pascals_n2 * crush_pressure_adjust_ratio_n2;
        pVpm->adjusted_crushing_pressure_he[i - 1] =
        adj_crush_pressure_he_pascals / 101325.0f * UNITS_FACTOR;
        pVpm->adjusted_crushing_pressure_n2[i - 1] =
        adj_crush_pressure_n2_pascals / 101325.0f * UNITS_FACTOR;
    }
    return 0;
} /* nuclear_regeneration */

/* =============================================================================== */
/*     SUBROUTINE CALC_INITIAL_ALLOWABLE_GRADIENT */
/*     Purpose: This subprogram calculates the initial allowable gradients for */
/*     helium and nitrogren in each compartment.  These are the gradients that */
/*     will be used to set the deco ceiling on the first pass through the deco */
/*     loop.  If the Critical Volume Algorithm is set to "off", then these */
/*     gradients will determine the final deco schedule.  Otherwise, if the */
/*     Critical Volume Algorithm is set to "on", these gradients will be further */
/*     "relaxed" by the Critical Volume Algorithm subroutine.  The initial */
/*     allowable gradients are referred to as "PssMin" in the papers by Yount */
/*     and colleauges, i.e., the minimum supersaturation pressure gradients */
/*     that will probe bubble formation in the VPM nuclei that started with the */
/*     designated minimum initial radius (critical radius). */

/*     The initial allowable gradients are computed directly from the */
/*     "regenerated" radii after the Nuclear Regeneration subroutine.  These */
/*     gradients are tracked separately for helium and nitrogen. */
/* =============================================================================== */

static int calc_initial_allowable_gradient()
{
    float initial_allowable_grad_n2_pa,
    initial_allowable_grad_he_pa;
    short i;

    /* loop */
    /* =============================================================================== */
    /*     CALCULATIONS */
    /*     The initial allowable gradients are computed in Pascals and then converted */
    /*     to the diving pressure units.  Two different sets of arrays are used to */
    /*     save the calculations - Initial Allowable Gradients and Allowable */
    /*     Gradients.  The Allowable Gradients are assigned the values from Initial */
    /*     Allowable Gradients however the Allowable Gradients can be changed later */
    /*     by the Critical Volume subroutine.  The values for the Initial Allowable */
    /*     Gradients are saved in a global array for later use by both the Critical */
    /*     Volume subroutine and the VPM Repetitive Algorithm subroutine. */
    /* =============================================================================== */

    for (i = 1; i <= 16; ++i)
    {
        initial_allowable_grad_n2_pa =
            SURFACE_TENSION_GAMMA * 2.0f *
            (SKIN_COMPRESSION_GAMMAC - SURFACE_TENSION_GAMMA) /
            (regenerated_radius_n2[i - 1] * SKIN_COMPRESSION_GAMMAC);

        initial_allowable_grad_he_pa =
            SURFACE_TENSION_GAMMA * 2.0f *
            (SKIN_COMPRESSION_GAMMAC - SURFACE_TENSION_GAMMA) /
            (regenerated_radius_he[i - 1] * SKIN_COMPRESSION_GAMMAC);

        pVpm->initial_allowable_gradient_n2[i - 1] =
            initial_allowable_grad_n2_pa / 101325.0f * UNITS_FACTOR;

        pVpm->initial_allowable_gradient_he[i - 1] =
            initial_allowable_grad_he_pa / 101325.0f * UNITS_FACTOR;

        allowable_gradient_he[i - 1] =
            pVpm->initial_allowable_gradient_he[i - 1];

        allowable_gradient_n2[i - 1] =
            pVpm->initial_allowable_gradient_n2[i - 1];
    }
    return 0;
} /* calc_initial_allowable_gradient */

/* =============================================================================== */
/*     SUBROUTINE CALC_DECO_CEILING */
/*     Purpose: This subprogram calculates the deco ceiling (the safe ascent */
/*     depth) in each compartment, based on the allowable gradients, and then */
/*     finds the deepest deco ceiling across all compartments.  This deepest */
/*     value (Deco Ceiling Depth) is then used by the Decompression Stop */
/*     subroutine to determine the actual deco schedule. */
/* =============================================================================== */

static int calc_deco_ceiling(float *deco_ceiling_depth,_Bool fallowable)
{
    /* System generated locals */
    float r1, r2;
    /* Local variables */
    float weighted_allowable_gradient;
    short i;
    float compartment_deco_ceiling[16],
    gas_loading,
    tolerated_ambient_pressure;
    float gradient_he, gradient_n2;

    if(!vpm_b)
        fallowable = true;
    /* loop */
    /* =============================================================================== */
    /*     CALCULATIONS */
    /*     Since there are two sets of allowable gradients being tracked, one for */
    /*     helium and one for nitrogen, a "weighted allowable gradient" must be */
    /*     computed each time based on the proportions of helium and nitrogen in */
    /*     each compartment.  This proportioning follows the methodology of */
    /*     Buhlmann/Keller.  If there is no helium and nitrogen in the compartment, */
    /*     such as after extended periods of oxygen breathing, then the minimum value */
    /*     across both gases will be used.  It is important to note that if a */
    /*     compartment is empty of helium and nitrogen, then the weighted allowable */
    /*     gradient formula cannot be used since it will result in division by zero. */
    /* =============================================================================== */

    for (i = 1; i <= 16; ++i)
    {

        // abfrage raus und pointer stattdessen
        if(fallowable){
            gradient_he = allowable_gradient_he[i-1];
            gradient_n2 = allowable_gradient_n2[i-1];
        }
        else{
            gradient_he = deco_gradient_he[i-1];
            gradient_n2 = deco_gradient_n2[i-1];
        }

        gas_loading =	helium_pressure[i - 1] + nitrogen_pressure[i - 1];

        if (gas_loading > 0)
        {
            weighted_allowable_gradient =
            (gradient_he * helium_pressure[i - 1] +
            gradient_n2 * nitrogen_pressure[i - 1]) /
            (helium_pressure[i - 1] + nitrogen_pressure[i - 1]);

            tolerated_ambient_pressure =
            gas_loading +
            CONSTANT_PRESSURE_OTHER_GASES -
            weighted_allowable_gradient;
        }
        else
        {
            /* Computing MIN */
            r1 = gradient_he;
            r2 = gradient_n2;
            weighted_allowable_gradient = fminf(r1,r2);

            tolerated_ambient_pressure =
            CONSTANT_PRESSURE_OTHER_GASES - weighted_allowable_gradient;
        }

        /* =============================================================================== */
        /*     The tolerated ambient pressure cannot be less than zero absolute, i.e., */
        /*     the vacuum of outer space! */
        /* =============================================================================== */

        if (tolerated_ambient_pressure < 0) {
            tolerated_ambient_pressure = 0;
        }
        compartment_deco_ceiling[i - 1] =
        tolerated_ambient_pressure - barometric_pressure;
    }

    /* =============================================================================== */
    /*     The Deco Ceiling Depth is computed in a loop after all of the individual */
    /*     compartment deco ceilings have been calculated.  It is important that the */
    /*     Deco Ceiling Depth (max deco ceiling across all compartments) only be */
    /*     extracted from the compartment values and not be compared against some */
    /*     initialization value.  For example, if MAX(Deco_Ceiling_Depth . .) was */
    /*     compared against zero, this could cause a program lockup because sometimes */
    /*     the Deco Ceiling Depth needs to be negative (but not less than zero */
    /*     absolute ambient pressure) in order to decompress to the last stop at zero */
    /*     depth. */
    /* =============================================================================== */

    *deco_ceiling_depth = compartment_deco_ceiling[0];
    for (i = 2; i <= 16; ++i)
    {
        /* Computing MAX */
        r1 = *deco_ceiling_depth;
        r2 = compartment_deco_ceiling[i - 1];
        *deco_ceiling_depth = fmaxf(r1,r2);
    }
    return 0;
} /* calc_deco_ceiling */



/* =============================================================================== */
/*     SUBROUTINE CALC_MAX_ACTUAL_GRADIENT */
/*     Purpose: This subprogram calculates the actual supersaturation gradient */
/*     obtained in each compartment as a result of the ascent profile during */
/*     decompression.  Similar to the concept with crushing pressure, the */
/*     supersaturation gradients are not cumulative over a multi-level, staged */
/*     ascent.  Rather, it will be the maximum value obtained in any one discrete */
/*     step of the overall ascent.  Thus, the program must compute and store the */
/*     maximum actual gradient for each compartment that was obtained across all */
/*     steps of the ascent profile.  This subroutine is invoked on the last pass */
/*     through the deco stop loop block when the final deco schedule is being */
/*     generated. */
/*   */
/*     The max actual gradients are later used by the VPM Repetitive Algorithm to */
/*     determine if adjustments to the critical radii are required.  If the max */
/*     actual gradient did not exceed the initial alllowable gradient, then no */
/*     adjustment will be made.  However, if the max actual gradient did exceed */
/*     the intitial allowable gradient, such as permitted by the Critical Volume */
/*     Algorithm, then the critical radius will be adjusted (made larger) on the */
/*     repetitive dive to compensate for the bubbling that was allowed on the */
/*     previous dive.  The use of the max actual gradients is intended to prevent */
/*     the repetitive algorithm from being overly conservative. */
/* =============================================================================== */

static int calc_max_actual_gradient(float *deco_stop_depth)
{
    /* System generated locals */
    float r1;

    /* Local variables */
    short i;
    float compartment_gradient;

    /* loop */
    /* =============================================================================== */
    /*     CALCULATIONS */
    /*     Note: negative supersaturation gradients are meaningless for this */
    /*     application, so the values must be equal to or greater than zero. */
    /* =============================================================================== */

    for (i = 1; i <= 16; ++i)
    {
        compartment_gradient =
        helium_pressure[i - 1] +
        nitrogen_pressure[i - 1] +
        CONSTANT_PRESSURE_OTHER_GASES -
        (*deco_stop_depth + barometric_pressure);
        if (compartment_gradient <= 0.0f) {
            compartment_gradient = 0.0f;
        }
        /* Computing MAX */
        r1 = pVpm->max_actual_gradient[i - 1];
        pVpm->max_actual_gradient[i - 1] = fmaxf(r1, compartment_gradient);
    }
    return 0;
} /* calc_max_actual_gradient */

/* =============================================================================== */
/*     SUBROUTINE CALC_SURFACE_PHASE_VOLUME_TIME */
/*     Purpose: This subprogram computes the surface portion of the total phase */
/*     volume time.  This is the time factored out of the integration of */
/*     supersaturation gradient x time over the surface interval.  The VPM */
/*     considers the gradients that allow bubbles to form or to drive bubble */
/*     growth both in the water and on the surface after the dive. */

/*     This subroutine is a new development to the VPM algorithm in that it */
/*     computes the time course of supersaturation gradients on the surface */
/*     when both helium and nitrogen are present.  Refer to separate write-up */
/*     for a more detailed explanation of this algorithm. */
/* =============================================================================== */

static int calc_surface_phase_volume_time()
{
    /* Local variables */
    float decay_time_to_zero_gradient;
    short i;
    float integral_gradient_x_time,
    surface_inspired_n2_pressure;

    /* loop */
    /* =============================================================================== */
    /*     CALCULATIONS */
    /* =============================================================================== */

    surface_inspired_n2_pressure =
    (barometric_pressure - WATER_VAPOR_PRESSURE) * 0.79f;
    for (i = 1; i <= 16; ++i)
    {
        if (nitrogen_pressure[i - 1] > surface_inspired_n2_pressure)
        {
            surface_phase_volume_time[i - 1] =
            (helium_pressure[i - 1] / HELIUM_TIME_CONSTANT[i - 1] +
            (nitrogen_pressure[i - 1] - surface_inspired_n2_pressure) /
            NITROGEN_TIME_CONSTANT[i - 1]) /
            (helium_pressure[i - 1] + nitrogen_pressure[i - 1] -
            surface_inspired_n2_pressure);
        } else if (nitrogen_pressure[i - 1] <= surface_inspired_n2_pressure &&
        helium_pressure[i - 1] + nitrogen_pressure[i - 1] >= surface_inspired_n2_pressure)
        {
            decay_time_to_zero_gradient =
            1.0f / (NITROGEN_TIME_CONSTANT[i - 1] - HELIUM_TIME_CONSTANT[i - 1]) *
            log((surface_inspired_n2_pressure - nitrogen_pressure[i - 1]) /
            helium_pressure[i - 1]);
            integral_gradient_x_time =
            helium_pressure[i - 1] /
            HELIUM_TIME_CONSTANT[i - 1] *
            (1.0f - expf(-HELIUM_TIME_CONSTANT[i - 1] *
            decay_time_to_zero_gradient)) +
            (nitrogen_pressure[i - 1] - surface_inspired_n2_pressure) /
            NITROGEN_TIME_CONSTANT[i - 1] *
            (1.0f - expf(-NITROGEN_TIME_CONSTANT[i - 1] *
            decay_time_to_zero_gradient));
            surface_phase_volume_time[i - 1] =
            integral_gradient_x_time /
            (helium_pressure[i - 1] +
            nitrogen_pressure[i - 1] -
            surface_inspired_n2_pressure);
        } else {
            surface_phase_volume_time[i - 1] = 0.0f;
        }
    }
    return 0;
} /* calc_surface_phase_volume_time */

/* =============================================================================== */
/*     SUBROUTINE CRITICAL_VOLUME */
/*     Purpose: This subprogram applies the VPM Critical Volume Algorithm.  This */
/*     algorithm will compute "relaxed" gradients for helium and nitrogen based */
/*     on the setting of the Critical Volume Parameter Lambda. */
/* =============================================================================== */

static int critical_volume(float *deco_phase_volume_time)
{
    /* System generated locals */
    float r1;

    /* Local variables */
    float initial_allowable_grad_n2_pa,
    initial_allowable_grad_he_pa,
    parameter_lambda_pascals, b,
    c;
    short i;
    float new_allowable_grad_n2_pascals,
    phase_volume_time[16],
    new_allowable_grad_he_pascals,
    adj_crush_pressure_n2_pascals,
    adj_crush_pressure_he_pascals;

    /* loop */
    /* =============================================================================== */
    /*     CALCULATIONS */
    /*     Note:  Since the Critical Volume Parameter Lambda was defined in units of */
    /*     fsw-min in the original papers by Yount and colleauges, the same */
    /*     convention is retained here.  Although Lambda is adjustable only in units */
    /*     of fsw-min in the program settings (range from 6500 to 8300 with default */
    /*     7500), it will convert to the proper value in Pascals-min in this */
    /*     subroutine regardless of which diving pressure units are being used in */
    /*     the main program - feet of seawater (fsw) or meters of seawater (msw). */
    /*     The allowable gradient is computed using the quadratic formula (refer to */
    /*     separate write-up posted on the Deco List web site). */
    /* =============================================================================== */

/**
    ******************************************************************************
    * @brief		critical_volume comment by hw
    * @version 	V0.0.1
    * @date   	19-April-2014
    * @retval 	global: allowable_gradient_he[i], allowable_gradient_n2[i]
    ******************************************************************************
    */

    parameter_lambda_pascals =
    CRIT_VOLUME_PARAMETER_LAMBDA / 33.0f * 101325.0f;
    for (i = 1; i <= 16; ++i)
    {
        phase_volume_time[i - 1] =
        *deco_phase_volume_time + surface_phase_volume_time[i - 1];
    }
    for (i = 1; i <= 16; ++i)
    {

        adj_crush_pressure_he_pascals =
        pVpm->adjusted_crushing_pressure_he[i - 1] / UNITS_FACTOR * 101325.0f;

        initial_allowable_grad_he_pa =
        pVpm->initial_allowable_gradient_he[i - 1] / UNITS_FACTOR * 101325.0f;

        b = initial_allowable_grad_he_pa + parameter_lambda_pascals *
        SURFACE_TENSION_GAMMA / (
        SKIN_COMPRESSION_GAMMAC * phase_volume_time[i - 1]);

        c = SURFACE_TENSION_GAMMA * (
        SURFACE_TENSION_GAMMA * (
        parameter_lambda_pascals * adj_crush_pressure_he_pascals)) /
        (SKIN_COMPRESSION_GAMMAC *
        (SKIN_COMPRESSION_GAMMAC * phase_volume_time[i - 1]));
        /* Computing 2nd power */

        r1 = b;

        new_allowable_grad_he_pascals =
        (b + sqrtf(r1 * r1 - c * 4.0f)) / 2.0f;

        /* modify global variable */
        allowable_gradient_he[i - 1] =
        new_allowable_grad_he_pascals / 101325.0f * UNITS_FACTOR;
    }

    for (i = 1; i <= 16; ++i)
    {
        adj_crush_pressure_n2_pascals =
        pVpm->adjusted_crushing_pressure_n2[i - 1] / UNITS_FACTOR * 101325.0f;

        initial_allowable_grad_n2_pa =
        pVpm->initial_allowable_gradient_n2[i - 1] / UNITS_FACTOR * 101325.0f;

        b = initial_allowable_grad_n2_pa + parameter_lambda_pascals *
        SURFACE_TENSION_GAMMA / (
        SKIN_COMPRESSION_GAMMAC * phase_volume_time[i - 1]);

        c = SURFACE_TENSION_GAMMA *
        (SURFACE_TENSION_GAMMA *
        (parameter_lambda_pascals * adj_crush_pressure_n2_pascals)) /
        (SKIN_COMPRESSION_GAMMAC *
        (SKIN_COMPRESSION_GAMMAC * phase_volume_time[i - 1]));
        /* Computing 2nd power */

        r1 = b;

        new_allowable_grad_n2_pascals =
        (b + sqrtf(r1 * r1 - c * 4.0f)) / 2.0f;

        /* modify global variable */
        allowable_gradient_n2[i - 1] =
        new_allowable_grad_n2_pascals / 101325.0f * UNITS_FACTOR;
    }
    return 0;
} /* critical_volume */

/* =============================================================================== */
/*     SUBROUTINE CALC_START_OF_DECO_ZONE */
/*     Purpose: This subroutine uses the Bisection Method to find the depth at */
/*     which the leading compartment just enters the decompression zone. */
/*     Source:  "Numerical Recipes in Fortran 77", Cambridge University Press, */
/*     1992. */
/* =============================================================================== */

static int calc_start_of_deco_zone(float *starting_depth,
                            float *rate,
                            float *depth_start_of_deco_zone)
{
    /* Local variables */
    float last_diff_change,
    initial_helium_pressure,
    mid_range_nitrogen_pressure;
    short i, j;
    float initial_inspired_n2_pressure,
    cpt_depth_start_of_deco_zone,
    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,
    starting_ambient_pressure,
    initial_nitrogen_pressure,
    function_at_high_bound;

    float time_to_start_of_deco_zone,
    high_bound_helium_pressure,
    helium_rate,
    differential_change;
    float fraction_helium_begin;
    float fraction_helium_end;
    float fraction_nitrogen_begin;
    float fraction_nitrogen_end;
    float ending_ambient_pressure;
    float time_test;


    /* loop */
    /* =============================================================================== */
    /*     CALCULATIONS */
    /*     First initialize some variables */
    /* =============================================================================== */

    *depth_start_of_deco_zone = 0.0f;
    starting_ambient_pressure =	*starting_depth + barometric_pressure;

    //>>>>>>>>>>>>>>>>>>>>
    //Test depth to calculate helium_rate and nitrogen_rate
    ending_ambient_pressure = starting_ambient_pressure/2;

    time_test = (ending_ambient_pressure - starting_ambient_pressure) / *rate;
    decom_get_inert_gases(starting_ambient_pressure / 10, (&pDiveSettings->decogaslist[mix_number]), &fraction_nitrogen_begin, &fraction_helium_begin );
    decom_get_inert_gases(ending_ambient_pressure   / 10, (&pDiveSettings->decogaslist[mix_number]), &fraction_nitrogen_end, &fraction_helium_end );
    initial_inspired_he_pressure =	(starting_ambient_pressure - WATER_VAPOR_PRESSURE) * fraction_helium_begin;
    initial_inspired_n2_pressure =	(starting_ambient_pressure - WATER_VAPOR_PRESSURE) * fraction_nitrogen_begin;
    helium_rate = ((ending_ambient_pressure  - WATER_VAPOR_PRESSURE)* fraction_helium_end - initial_inspired_he_pressure)/time_test;
    nitrogen_rate = ((ending_ambient_pressure  - WATER_VAPOR_PRESSURE)* fraction_nitrogen_end - initial_inspired_n2_pressure)/time_test;
    //>>>>>>>>>>>>>>>>>>>>>
    /*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];*/

    /* =============================================================================== */
    /*     ESTABLISH THE BOUNDS FOR THE ROOT SEARCH USING THE BISECTION METHOD */
    /*     AND CHECK TO MAKE SURE THAT THE ROOT WILL BE WITHIN BOUNDS.  PROCESS */
    /*     EACH COMPARTMENT INDIVIDUALLY AND FIND THE MAXIMUM DEPTH ACROSS ALL */
    /*     COMPARTMENTS (LEADING COMPARTMENT) */
    /*     In this case, we are solving for time - the time when the gas tension in */
    /*     the compartment will be equal to ambient pressure.  The low bound for time */
    /*     is set at zero and the high bound is set at the time it would take to */
    /*     ascend to zero ambient pressure (absolute).  Since the ascent rate is */
    /*     negative, a multiplier of -1.0 is used to make the time positive.  The */
    /*     desired point when gas tension equals ambient pressure is found at a time */
    /*     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. */
    /* =============================================================================== */

    low_bound = 0.;
    high_bound = starting_ambient_pressure / *rate * -1.0f;
    for (i = 1; i <= 16; ++i)
    {
        initial_helium_pressure = helium_pressure[i - 1];
        initial_nitrogen_pressure = nitrogen_pressure[i - 1];
        function_at_low_bound =
        initial_helium_pressure +
        initial_nitrogen_pressure +
        CONSTANT_PRESSURE_OTHER_GASES -
        starting_ambient_pressure;
        high_bound_helium_pressure =
        schreiner_equation__2(&initial_inspired_he_pressure,
        &helium_rate,
        &high_bound,
        &HELIUM_TIME_CONSTANT[i - 1],
        &initial_helium_pressure);
        high_bound_nitrogen_pressure =
        schreiner_equation__2(&initial_inspired_n2_pressure,
        &nitrogen_rate,
        &high_bound,
        &NITROGEN_TIME_CONSTANT[i - 1],
        &initial_nitrogen_pressure);
        function_at_high_bound = high_bound_helium_pressure +
        high_bound_nitrogen_pressure +
        CONSTANT_PRESSURE_OTHER_GASES;
        if (function_at_high_bound * function_at_low_bound >= 0.0f)
        {
            printf("\nERROR! ROOT IS NOT WITHIN BRACKETS");
        }

        /* =============================================================================== */
        /*     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.0f)
        {
            time_to_start_of_deco_zone = low_bound;
            differential_change = high_bound - low_bound;
        } else {
            time_to_start_of_deco_zone = high_bound;
            differential_change = low_bound - high_bound;
        }
        for (j = 1; j <= 100; ++j)
        {
            last_diff_change = differential_change;
            differential_change = last_diff_change * 0.5f;
            mid_range_time =
            time_to_start_of_deco_zone +
            differential_change;
            mid_range_helium_pressure =
            schreiner_equation__2(&initial_inspired_he_pressure,
            &helium_rate,
            &mid_range_time,
            &HELIUM_TIME_CONSTANT[i - 1],
            &initial_helium_pressure);
            mid_range_nitrogen_pressure =
            schreiner_equation__2(&initial_inspired_n2_pressure,
            &nitrogen_rate,
            &mid_range_time,
            &NITROGEN_TIME_CONSTANT[i - 1],
            &initial_nitrogen_pressure);
            function_at_mid_range =
            mid_range_helium_pressure +
            mid_range_nitrogen_pressure +
            CONSTANT_PRESSURE_OTHER_GASES -
            (starting_ambient_pressure + *rate * mid_range_time);
            if (function_at_mid_range <= 0.0f) {
                time_to_start_of_deco_zone = mid_range_time;
            }
            if( fabs(differential_change) < 0.001f
             || function_at_mid_range == 0.0f)
            {
                goto L170;
            }
            /* L150: */
        }
        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 170 and assigns the solution value for the individual */
        /*     compartment. */
        /* =============================================================================== */

        L170:
        cpt_depth_start_of_deco_zone =
        starting_ambient_pressure +
        *rate * time_to_start_of_deco_zone -
        barometric_pressure;

        /* =============================================================================== */
        /*     The overall solution will be the compartment with the maximum depth where */
        /*     gas tension equals ambient pressure (leading compartment). */
        /* =============================================================================== */

        *depth_start_of_deco_zone =
        fmaxf(*depth_start_of_deco_zone, cpt_depth_start_of_deco_zone);
        /* L200: */
    }
    return 0;
} /* calc_start_of_deco_zone */

/* =============================================================================== */
/*     SUBROUTINE PROJECTED_ASCENT */
/*     Purpose: This subprogram performs a simulated ascent outside of the main */
/*     program to ensure that a deco ceiling will not be violated due to unusual */
/*     gas loading during ascent (on-gassing).  If the deco ceiling is violated, */
/*     the stop depth will be adjusted deeper by the step size until a safe */
/*     ascent can be made. */
/* =============================================================================== */

static int projected_ascent(float *starting_depth,
                     float *rate,
                     float *deco_stop_depth,
                     float *step_size)
{
    /* Local variables */
    float weighted_allowable_gradient,
    ending_ambient_pressure,
    temp_gas_loading[16];
    int i;
    float allowable_gas_loading[16];
    float temp_nitrogen_pressure[16];
    float temp_helium_pressure[16];
    float run_time_save = 0;

    /* loop */
    /* =============================================================================== */
    /*     CALCULATIONS */
    /* =============================================================================== */


    L665:
    ending_ambient_pressure = *deco_stop_depth + barometric_pressure;
    for (i = 1; i <= 16; ++i) {
        temp_helium_pressure[i - 1] = helium_pressure[i - 1];
        temp_nitrogen_pressure[i - 1] = nitrogen_pressure[i - 1];
    }
    run_time_save = run_time;
    gas_loadings_ascent_descen(temp_helium_pressure, temp_nitrogen_pressure, *starting_depth,*deco_stop_depth,*rate,true);
    run_time = run_time_save;

    for (i = 1; i <= 16; ++i)
    {
        temp_gas_loading[i - 1] =
        temp_helium_pressure[i - 1] +
        temp_nitrogen_pressure[i - 1];
        if (temp_gas_loading[i - 1] > 0.0f)
        {
            weighted_allowable_gradient =
            (allowable_gradient_he[i - 1] *
            temp_helium_pressure[i - 1] +
            allowable_gradient_n2[i - 1] *
            temp_nitrogen_pressure[i - 1]) / temp_gas_loading[i - 1];
        } else {
            /* Computing MIN */
            weighted_allowable_gradient = fminf(allowable_gradient_he[i - 1],allowable_gradient_n2[i - 1]);
        }
        allowable_gas_loading[i - 1] =
        ending_ambient_pressure +
        weighted_allowable_gradient -
        CONSTANT_PRESSURE_OTHER_GASES;
        /* L670: */
    }
    for (i = 1; i <= 16; ++i) {
        if (temp_gas_loading[i - 1] > allowable_gas_loading[i - 1]) {
            *deco_stop_depth += *step_size;
            goto L665;
        }
        /* L671: */
    }
    return 0;
} /* projected_ascent */

/* =============================================================================== */
/*     SUBROUTINE DECOMPRESSION_STOP */
/*     Purpose: This subprogram calculates the required time at each */
/*     decompression stop. */
/* =============================================================================== */

static void decompression_stop(float *deco_stop_depth,
                        float *step_size,
                        _Bool final_deco_calculation)
{
    /* Local variables */
    float inspired_nitrogen_pressure;
//	short last_segment_number;
//	float weighted_allowable_gradient;
    float initial_helium_pressure[16];
 /* by hw */
    float initial_CNS = gCNS_VPM;

    //static float time_counter;
    short i;
    float ambient_pressure;
    float inspired_helium_pressure,
    next_stop;
    //last_run_time,
    //temp_segment_time;

    float deco_ceiling_depth,
    initial_nitrogen_pressure[16];
    //round_up_operation;
    float fraction_helium_begin;
    float fraction_nitrogen_begin;
    int count = 0;
    _Bool buehlmann_wait = false;
    float tissue_He_saturation[16];
    float tissue_N2_saturation[16];
    float vpm_buehlmann_safety_gradient = 1.0f - (((float)pDiveSettings->vpm_conservatism) / 40);
    /* loop */
    /* =============================================================================== */
    /*     CALCULATIONS */
    /* =============================================================================== */

    segment_time = 0;
//	temp_segment_time = segment_time;
    ambient_pressure = *deco_stop_depth + barometric_pressure;
    //ending_ambient_pressure = ambient_pressure;
    decom_get_inert_gases(ambient_pressure / 10, (&pDiveSettings->decogaslist[mix_number]), &fraction_nitrogen_begin, &fraction_helium_begin );

    if(*deco_stop_depth == (float)(pDiveSettings->last_stop_depth_bar * 10))
        next_stop = 0;
    else
    {
        next_stop = *deco_stop_depth - *step_size;
        next_stop = fmaxf(next_stop,(float)pDiveSettings->last_stop_depth_bar * 10);
    }

    inspired_helium_pressure =
    (ambient_pressure - WATER_VAPOR_PRESSURE) * fraction_helium_begin;
    inspired_nitrogen_pressure =
    (ambient_pressure - WATER_VAPOR_PRESSURE) * fraction_nitrogen_begin;

    /* =============================================================================== */
    /*     Check to make sure that program won't lock up if unable to decompress */
    /*     to the next stop.  If so, write error message and terminate program. */
    /* =============================================================================== */

    //deco_ceiling_depth = next_stop +1;	//deco_ceiling_depth = next_stop + 1;
    if(!vpm_violates_buehlmann)
    {
        calc_deco_ceiling(&deco_ceiling_depth, false); //weg, weil auf jeden Fall schleife für safety und so konservativer
    }
    else
    {
        deco_ceiling_depth = next_stop + 1;
    }
        if(deco_ceiling_depth > next_stop)
        {
            while (deco_ceiling_depth > next_stop)
            {

                segment_time  +=  60;
                if(segment_time >= 999 )
                {
                    segment_time = 999 ;
                    run_time += segment_time;
                    return;
                }
                //goto L700;
                initial_CNS = gCNS_VPM;
                decom_oxygen_calculate_cns_exposure(60*60,&pDiveSettings->decogaslist[mix_number],ambient_pressure/10,&gCNS_VPM);
                for (i = 0; i < 16; i++)
                {
                    initial_helium_pressure[i] = helium_pressure[i];
                    initial_nitrogen_pressure[i] = nitrogen_pressure[i];
                    helium_pressure[i] 		+=	(inspired_helium_pressure - helium_pressure[i]) 		*	float_buehlmann_He_factor_expositon_one_hour[i];
                    nitrogen_pressure[i]	+= 	(inspired_nitrogen_pressure - nitrogen_pressure[i]) *	float_buehlmann_N2_factor_expositon_one_hour[i];
                }
                calc_deco_ceiling(&deco_ceiling_depth, false);
            }
            if(deco_ceiling_depth < next_stop)
            {
                segment_time  -=  60;
                gCNS_VPM = initial_CNS;
                for (i = 0; i < 16; i++)
                {
                    helium_pressure[i] = initial_helium_pressure[i];
                    nitrogen_pressure[i] = initial_nitrogen_pressure[i];
                }
                deco_ceiling_depth = next_stop +1;
            }
            count = 0;
            while (deco_ceiling_depth > next_stop && count < 13)
            {
                count++;
                segment_time  +=  5;
                //goto L700;
                initial_CNS = gCNS_VPM;
                decom_oxygen_calculate_cns_exposure(60*5,&pDiveSettings->decogaslist[mix_number],ambient_pressure/10,&gCNS_VPM);
                for (i = 0; i < 16; i++)
                {
                    initial_helium_pressure[i] = helium_pressure[i];
                    initial_nitrogen_pressure[i] = nitrogen_pressure[i];
                    helium_pressure[i] 		+=	(inspired_helium_pressure - helium_pressure[i]) 		*	float_buehlmann_He_factor_expositon_five_minutes[i];
                    nitrogen_pressure[i]	+= 	(inspired_nitrogen_pressure - nitrogen_pressure[i]) *	float_buehlmann_N2_factor_expositon_five_minutes[i];
                }
                calc_deco_ceiling(&deco_ceiling_depth, false);
            }
            if(deco_ceiling_depth < next_stop)
            {
                segment_time  -=  5;
                gCNS_VPM = initial_CNS;
                for (i = 0; i < 16; i++) {
                    helium_pressure[i] = initial_helium_pressure[i];
                    nitrogen_pressure[i] = initial_nitrogen_pressure[i];
                }
                deco_ceiling_depth = next_stop +1;
            }
            buehlmann_wait = false;
            while (buehlmann_wait || (deco_ceiling_depth > next_stop))
            {
                //time_counter = temp_segment_time;
                segment_time  +=  1;

                if(segment_time >= 999 )
                {
                    segment_time = 999 ;
                    run_time += segment_time;
                    return;
                }
                //goto L700;
                initial_CNS = gCNS_VPM;
                decom_oxygen_calculate_cns_exposure(60*1,&pDiveSettings->decogaslist[mix_number],ambient_pressure/10,&gCNS_VPM);
                for (i = 0; i < 16; i++)
                {
                    initial_helium_pressure[i] = helium_pressure[i];
                    initial_nitrogen_pressure[i] = nitrogen_pressure[i];
                    helium_pressure[i] 		+=	(inspired_helium_pressure - helium_pressure[i]) 		*	float_buehlmann_He_factor_expositon_one_minute[i];
                    nitrogen_pressure[i]	+= 	(inspired_nitrogen_pressure - nitrogen_pressure[i]) *	float_buehlmann_N2_factor_expositon_one_minute[i];
                }
                if(!buehlmann_wait)
                    calc_deco_ceiling(&deco_ceiling_depth, false);

                if(buehlmannSafety && final_deco_calculation && !(deco_ceiling_depth > next_stop))
                {
                    for (i = 0; i < 16; i++)
                    {
                        tissue_He_saturation[i] = helium_pressure[i] / 10;
                        tissue_N2_saturation[i] = nitrogen_pressure[i] / 10;
                    }
                    if( (fabsf(nitrogen_pressure[15] - inspired_nitrogen_pressure) < 0.00001f) && (fabsf(helium_pressure[15] - inspired_helium_pressure) < 0.00001f)
                     && (fabsf(nitrogen_pressure[0] - inspired_nitrogen_pressure) < 0.00001f) && (fabsf(helium_pressure[0] - inspired_helium_pressure) < 0.00001f))
                    {
                         buehlmann_wait_exceeded = true;
                        break;
                    }

                    if(decom_tissue_test_tolerance(tissue_N2_saturation, tissue_He_saturation, vpm_buehlmann_safety_gradient, (next_stop / 10.0f) + pInput->pressure_surface_bar))
                        break;

                    buehlmann_wait = true;
                }
            }
            if(buehlmann_wait)
            {
                vpm_violates_buehlmann = true;
            }
            if(!buehlmann_wait)
            {
                if(deco_ceiling_depth < next_stop)
                {
                    segment_time  -=  1;
                    gCNS_VPM = initial_CNS;
                    for (i = 0; i < 16; i++) {
                        helium_pressure[i] = initial_helium_pressure[i];
                        nitrogen_pressure[i] = initial_nitrogen_pressure[i];
                    }
                    deco_ceiling_depth = next_stop +1;
                }
                while (deco_ceiling_depth > next_stop)
                {
                    //time_counter = temp_segment_time;
                    segment_time  +=  (float) 1.0f / 3.0f;
                    //goto L700;
                    initial_CNS = gCNS_VPM;
                    decom_oxygen_calculate_cns_exposure(20,&pDiveSettings->decogaslist[mix_number],ambient_pressure/10,&gCNS_VPM);
                    for (i = 0; i < 16; i++)
                    {
                        helium_pressure[i] 		+=	(inspired_helium_pressure - helium_pressure[i]) 		*	float_buehlmann_He_factor_expositon_20_seconds[i];
                        nitrogen_pressure[i]	+= 	(inspired_nitrogen_pressure - nitrogen_pressure[i]) *	float_buehlmann_N2_factor_expositon_20_seconds[i];
                    }
                    calc_deco_ceiling(&deco_ceiling_depth, false);
                }
            }
    }

/*float pressure_save =dive_data.pressure;
dive_data.pressure = ambient_pressure/10;
tissues_exposure_stage(st_deco_test,(int)(segment_time * 60), &dive_data, &gaslist);
dive_data.pressure = pressure_save;*/
    run_time += segment_time;
    return;
} /* decompression_stop */

/* =============================================================================== */
// SUROUTINE BOYLES_LAW_COMPENSATION
//     Purpose: This subprogram calculates the reduction in allowable gradients
//    with decreasing ambient pressure during the decompression profile based
//    on Boyle's Law considerations.
//===============================================================================
static void BOYLES_LAW_COMPENSATION (float* First_Stop_Depth,
                              float*  Deco_Stop_Depth,
                              float* Step_Size)
{
    short i;

    float Next_Stop;
    float Ambient_Pressure_First_Stop, Ambient_Pressure_Next_Stop;
    float Amb_Press_First_Stop_Pascals, Amb_Press_Next_Stop_Pascals;
    float A, B, C, Low_Bound, High_Bound, Ending_Radius;
    float Deco_Gradient_Pascals;
    float Allow_Grad_First_Stop_He_Pa, Radius_First_Stop_He;
    float Allow_Grad_First_Stop_N2_Pa, Radius_First_Stop_N2;

    //===============================================================================
    //     LO//AL ARRAYS
    //===============================================================================
//	float Radius1_He[16], Radius2_He[16];
//	float Radius1_N2[16], Radius2_N2[16];
    float root_factor;

    //===============================================================================
    //     CALCULATIONS
    //===============================================================================
    Next_Stop = *Deco_Stop_Depth - *Step_Size;

    Ambient_Pressure_First_Stop = *First_Stop_Depth +
                                 barometric_pressure;

    Ambient_Pressure_Next_Stop = Next_Stop + barometric_pressure;

    Amb_Press_First_Stop_Pascals =  (Ambient_Pressure_First_Stop/UNITS_FACTOR) * 101325.0f;

    Amb_Press_Next_Stop_Pascals =
            (Ambient_Pressure_Next_Stop/UNITS_FACTOR) * 101325.0f;
    root_factor = powf(Amb_Press_First_Stop_Pascals/Amb_Press_Next_Stop_Pascals,1.0f / 3.0f);

    for( i = 0; i < 16;i++)
    {
            Allow_Grad_First_Stop_He_Pa =
                        (allowable_gradient_he[i]/UNITS_FACTOR) * 101325.0f;

            Radius_First_Stop_He = (2.0f * SURFACE_TENSION_GAMMA) /
                                    Allow_Grad_First_Stop_He_Pa;

//			Radius1_He[i] = Radius_First_Stop_He;
            A = Amb_Press_Next_Stop_Pascals;
            B = -2.0f * SURFACE_TENSION_GAMMA;
            C = (Amb_Press_First_Stop_Pascals + (2.0f * SURFACE_TENSION_GAMMA)/
                     Radius_First_Stop_He)* Radius_First_Stop_He*
                     (Radius_First_Stop_He*(Radius_First_Stop_He));
            Low_Bound = Radius_First_Stop_He;
            High_Bound = Radius_First_Stop_He * root_factor;
                        //*pow(Amb_Press_First_Stop_Pascals/Amb_Press_Next_Stop_Pascals,1.0/3.0);
                    //*(Amb_Press_First_Stop_Pascals/Amb_Press_Next_Stop_Pascals)**(1.0/3.0);

            radius_root_finder(&A,&B,&C, &Low_Bound, &High_Bound,
                                             &Ending_Radius);

//			Radius2_He[i] = Ending_Radius;
            Deco_Gradient_Pascals = (2.0f * SURFACE_TENSION_GAMMA) /
                                    Ending_Radius;

            deco_gradient_he[i] = (Deco_Gradient_Pascals / 101325.0f)*
                                     UNITS_FACTOR;

    }

    for( i = 0; i < 16;i++)
    {
        Allow_Grad_First_Stop_N2_Pa =
                             (allowable_gradient_n2[i]/UNITS_FACTOR) * 101325.0f;

        Radius_First_Stop_N2 = (2.0f * SURFACE_TENSION_GAMMA) /
                                                         Allow_Grad_First_Stop_N2_Pa;

//		Radius1_N2[i] = Radius_First_Stop_N2;
        A = Amb_Press_Next_Stop_Pascals;
        B = -2.0f * SURFACE_TENSION_GAMMA;
        C = (Amb_Press_First_Stop_Pascals + (2.0f * SURFACE_TENSION_GAMMA)/
                        Radius_First_Stop_N2)* Radius_First_Stop_N2*
                        (Radius_First_Stop_N2*(Radius_First_Stop_N2));
        Low_Bound = Radius_First_Stop_N2;
        High_Bound = Radius_First_Stop_N2* root_factor;//pow(Amb_Press_First_Stop_Pascals/Amb_Press_Next_Stop_Pascals,1.0/3.0);

        //High_Bound = Radius_First_Stop_N2*exp(log(Amb_Press_First_Stop_Pascals/Amb_Press_Next_Stop_Pascals)/3);
        radius_root_finder(&A,&B,&C, &Low_Bound, &High_Bound,
                                                                        &Ending_Radius);

//		Radius2_N2[i] = Ending_Radius;
        Deco_Gradient_Pascals = (2.0f * SURFACE_TENSION_GAMMA) /
                                                         Ending_Radius;

        deco_gradient_n2[i] = (Deco_Gradient_Pascals / 101325.0f)*
                                                        UNITS_FACTOR;
    }
}

/* =============================================================================== */
//	vpm_calc_ndl
//     Purpose: This function computes NDL (time where no decostops are needed)
//===============================================================================
#define MAX_NDL	240

static int vpm_calc_ndl(void)
{
    static float future_helium_pressure[16];
    static float future_nitrogen_pressure[16];
    static int temp_segment_time;
    static int mix_number;
    static float inspired_helium_pressure;
    static float inspired_nitrogen_pressure;

    float previous_helium_pressure[16];
    float previous_nitrogen_pressure[16];
    float ambient_pressure;
    float fraction_helium_begin;
    float fraction_nitrogen_begin;
    int i = 0;
    int count = 0;
    int status = CALC_END;

        for(i = 0; i < 16;i++)
        {
            future_helium_pressure[i] =  pInput->tissue_helium_bar[i] * 10;//tissue_He_saturation[st_dive][i] * 10;
            future_nitrogen_pressure[i] = pInput->tissue_nitrogen_bar[i] * 10;
        }
        temp_segment_time = 0;

        mix_number = 0;
        ambient_pressure = pInput->pressure_ambient_bar  * 10;
        decom_get_inert_gases( ambient_pressure / 10,  (&pDiveSettings->decogaslist[mix_number]) , &fraction_nitrogen_begin, &fraction_helium_begin );
        inspired_helium_pressure =(ambient_pressure - WATER_VAPOR_PRESSURE) * fraction_helium_begin;
        inspired_nitrogen_pressure =(ambient_pressure - WATER_VAPOR_PRESSURE) *fraction_nitrogen_begin;

            status = CALC_END;
            while (status == CALC_END)
            {
                count++;
                temp_segment_time  +=  60;
                if(temp_segment_time >= MAX_NDL)
                {
                    pDecoInfo->output_ndl_seconds = temp_segment_time * 60;
                    return CALC_NDL;
                }
                run_time += 60;
                //goto L700;
                for (i = 1; i <= 16; ++i) {
                    previous_helium_pressure[i-1] = future_helium_pressure[i - 1];
                    previous_nitrogen_pressure[i - 1] = future_nitrogen_pressure[i - 1];
                    future_helium_pressure[i - 1] = future_helium_pressure[i - 1] + (inspired_helium_pressure - future_helium_pressure[i - 1]) * float_buehlmann_He_factor_expositon_one_hour[i-1];
                    future_nitrogen_pressure[i - 1] = future_nitrogen_pressure[i - 1] + (inspired_nitrogen_pressure - future_nitrogen_pressure[i - 1]) * float_buehlmann_N2_factor_expositon_one_hour[i-1];
                    helium_pressure[i - 1] = future_helium_pressure[i - 1];
                    nitrogen_pressure[i - 1] = future_nitrogen_pressure[i - 1];
                }
                vpm_calc_deco();
                while((status = vpm_calc_critcal_volume(true,true)) == CALC_CRITICAL);

            }

            temp_segment_time  -=  60;
            run_time -= 60;
            for (i = 1; i <= 16; ++i)
            {
                future_helium_pressure[i - 1] = previous_helium_pressure[i-1];
                future_nitrogen_pressure[i - 1] = previous_nitrogen_pressure[i - 1];
            }

        status = CALC_END;
        if(temp_segment_time < 60)
            nullzeit_unter60 = true;

        while (status == CALC_END)
        {
            temp_segment_time  +=  5;
            if(temp_segment_time >= MAX_NDL)
            {
                pDecoInfo->output_ndl_seconds = temp_segment_time * 60;
                return CALC_NDL;
            }
            if(nullzeit_unter60 && temp_segment_time > 60)
            {
                nullzeit_unter60 = false;
                return CALC_NDL;
            }
            run_time += 5;
            //goto L700;
            for (i = 1; i <= 16; ++i) {
                previous_helium_pressure[i-1] = future_helium_pressure[i - 1];
                previous_nitrogen_pressure[i - 1] = future_nitrogen_pressure[i - 1];
                future_helium_pressure[i - 1] = future_helium_pressure[i - 1] + (inspired_helium_pressure - future_helium_pressure[i - 1]) * float_buehlmann_He_factor_expositon_five_minutes[i-1];
                future_nitrogen_pressure[i - 1] = future_nitrogen_pressure[i - 1] + (inspired_nitrogen_pressure - future_nitrogen_pressure[i - 1]) * float_buehlmann_N2_factor_expositon_five_minutes[i-1];
                helium_pressure[i - 1] = future_helium_pressure[i - 1];
                nitrogen_pressure[i - 1] = future_nitrogen_pressure[i - 1];
            }
            vpm_calc_deco();
            while((status =vpm_calc_critcal_volume(true,true)) == CALC_CRITICAL);
        }
        temp_segment_time  -=  5;
        run_time -= 5;
        for (i = 1; i <= 16; ++i) {
            future_helium_pressure[i - 1] = previous_helium_pressure[i-1];
            future_nitrogen_pressure[i - 1] = previous_nitrogen_pressure[i - 1];
        }
        status = CALC_END;

    if(temp_segment_time <= 20)
    {
        while (status == CALC_END)
        {
            temp_segment_time  +=  minimum_deco_stop_time;
            run_time += minimum_deco_stop_time;
            //goto L700;
            for (i = 1; i <= 16; ++i) {
                future_helium_pressure[i - 1] = future_helium_pressure[i - 1] + (inspired_helium_pressure - future_helium_pressure[i - 1]) *	float_buehlmann_He_factor_expositon_one_minute[i-1];
                future_nitrogen_pressure[i - 1] = future_nitrogen_pressure[i - 1] + (inspired_nitrogen_pressure - future_nitrogen_pressure[i - 1]) * float_buehlmann_N2_factor_expositon_one_minute[i-1];
                helium_pressure[i - 1] = future_helium_pressure[i - 1];
                nitrogen_pressure[i - 1] =future_nitrogen_pressure[i - 1];

            }
            vpm_calc_deco();
            while((status =vpm_calc_critcal_volume(true,true)) == CALC_CRITICAL);

        }
    }
    else
        temp_segment_time += 5;
    pDecoInfo->output_ndl_seconds = temp_segment_time * 60;
    if(temp_segment_time > 1)
        return CALC_NDL;
    else
        return CALC_BEGIN;
}