Mercurial > public > ostc4
view Discovery/Src/vpm.c @ 589:8d95bab5bf6a
Minor: Update the documentation for building the code
| author | heinrichsweikamp |
|---|---|
| date | Thu, 17 Dec 2020 10:58:28 +0100 |
| 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; }