'*****************************************************************************
'                      - PROGRAM VPM.BAS (V 09132000)-
'
'                        Eric Maiken, (C)1989-2000
'
'                  Distribute freely. Credit the authors.
'
'Code for the algorithm of D.E. Yount and D.C. Hoffman,
'using simplifications discussed in the appendix to the paper:
'"On the Use of a Bubble Formation Model to Calculate Diving Tables"
'Aviation, Space, and Environmental Medicine, February, 1986.
'
'*****************************************************************************
'             NOTE: For output file, remove comments in SUB "putin" 
'*****************************************************************************
'
'
'
'*****************************************************************************
COMMON SHARED ceiling(), convflag, degast(), decotable(), dt(), halfT(), kN2()
COMMON SHARED ped, pet(), pnw(), pressure(), stagetension(), tension()
COMMON SHARED dfirst, g, gc, lambda, maxod, num, N2, nitroxStops, ncompartment, nstops, max
COMMON SHARED O2, oxygenStops, outfile$, outflag, P, pdeep, pin, psmin
COMMON SHARED t, tdeco, too
'
DECLARE FUNCTION gasExponential (timeconst, pfinal, pinitial, endtime, startime)
'
DECLARE SUB calc (k())
DECLARE SUB divedata ()
DECLARE SUB first ()
DECLARE SUB mindy ()
DECLARE SUB profile ()
DECLARE SUB pnew ()
DECLARE SUB putin ()
DECLARE SUB runtime ()
DECLARE SUB sort (matrix(), length)
DECLARE SUB stops (k())
'
'----------------------------------------------------------------------------
'
CLS
CALL putin                      'get input
CALL divedata                   'input depths and run-times, set constants
CLS
'
'----------------------------------------------------------------------------
' Initialize some arrays
'----------------------------------------------------------------------------
'
DIM tension(ncompartment)       'matrix to store compartment tensions at completion of all dive levels
DIM stagetension(ncompartment)  'matrix to store compartment tensions at completion of an ascent stage
DIM ceiling(ncompartment)       'matrix of ascent ceilings
DIM degast(ncompartment)        'matrix of degas times for each compartment at each stop
'
FOR i = 1 TO ncompartment       'zero all elements
 tension(i) = 0
 stagetension(i) = 0
 ceiling(i) = 0
 degast(i) = 0
NEXT i
'
DIM decotable(2, 40)            'dim array to hold time and depth of stops
FOR i = 1 TO 2                  'array is over dimensioned to 40 stops
 FOR j = 1 TO 40
  decotable(i, j) = 0
 NEXT j
NEXT i
'
'
'----------------------------------------------------------------------------
' Now do the iterative calculation of ascent schedule
'----------------------------------------------------------------------------
'
tdeco = 0                       'clear
convflag = 0                    'ascent schedule is not converged
iteration1flag = 0              'clear
'
CALL calc(kN2())                'calc ppN2 of compartments at end of multilevel dive
'                               'store in array tension()
'
CALL mindy                      'list of ascent ceilings from list of compartment tensions
'                               'using conservative compartment supersaturation
'       
CALL first                      'find depth of first stop
'
iteration1flag = 1              'first iteration
CALL stops(kN2())               'in or out gas compartments at series of stops
CALL profile                    'output
'
LOCATE 22, 1
PRINT "THIS IS THE FIRST, CONSERVATIVE ITERATION OF A DECOMPRESSION PROFILE"
PRINT "Press any key to continue"
DO: LOOP WHILE INKEY$ = ""
CLS
'
IF tdeco > 0 THEN                    'if this is a deco dive continue with calculation
200
         iteration1flag = 0
         tee = tdeco                 'total decompression time from last iteration
  
         CLS
         CALL pnew                   'update alowable compartment supersaturations to more liberal values

         FOR i = 1 TO ncompartment   'clear matricies from last iteration
          ceiling(i) = 0
         NEXT i
       
         CALL mindy                  'make a new list of ascent ceilings from list of compartment tensions
                                     'using relaxed compartment supersaturation
         CALL first
         CALL stops(kN2())
      
         diff = (tee - tdeco)        'difference in decompression times for last iteration
          
        IF diff >= 1 THEN

                GOTO 200             'loop through again until total stop time converges
        ELSE                         'go through one more iteration to store converged ascent schedule
         
                convflag = 1         'ascent schedule is converged
                CALL stops(kN2())
                CALL profile
        END IF                       'dif < 1, final iteration

ELSE                                 'tdeco = 0

 CALL profile

END IF

END

SUB calc (k())
'
'************************************************************************
'                            
'This sub tracks the ingassing of compartments on a multi-level dive
'
'************************************************************************
'
'
FOR j = 1 TO ncompartment    'do calc for different compartment half-times
  
  kay = kN2(j)
  pin = .79 * (1 - 102 / 760)  'begin with ppN2 =(Atmospheric fN2)*(1 ata -ppH2O - ppCO2)
                             'ppN2 updated @ end of j loop
  '
  FOR q = 1 TO num           'bootstrap gasexchange for sequence of multilevel depths
'
    too = pet(2, q)                        'tzero for ea depth step
    ped = N2 * (pet(1, q + 1) - 102 / 760) 'asymptotic ppN2 at ea level
   
    FOR t = pet(2, q) TO pet(2, q + 1) STEP (pet(2, q + 1) - pet(2, q))
      PN2 = gasExponential(kay, ped, pin, t, too)   'call exponential function for in or out gas
    NEXT t

    pin = PN2                'update: last multi level dive depth end pressure is new level start pressure
'
  NEXT q                     'next multi level dive depth
'
  tension(j) = PN2           'Save each compartment's final ppN2 in matrix for ascent calculation
'
NEXT j                       'next compartment
'
END SUB

SUB divedata
'**********************************************************************
'                         
'This sub stores constants for in and off gas calculations.
'
'Breathing gas %02 is hard-coded.
'
'Dive (depth/run-time) data is also input into the array dt() and converted
'to a pressure/run-time array pet().
'
'*****************************************************************************
'
'O2 partial pressures for dive and two ascent gases: a nitrox and O2 stage
'Think in terms of oxygen % when changing. Confirm output of %O2. You trak OTUs
'O2 <-> Dive % O2
'nitroxStops <-> stops %O2 limited by otus and mod
'oxygenStops <-> stops %O2 limited by otus and mod<=20 ffw
'
'----------------------------------------------------------------------------

O2 = .21                '%O2 For all Dive multi levels (num)
N2 = (1 - O2)           '%N2 during multilevel dive
'                       'and for all ascent depths > mod of nitroxStops ascent gas
'
'                    
nitroxStops = .21       '%02 for 20 fsw < stop depths <= maxod in SUB stops
maxod = 33 * (1.6 / nitroxStops - 1) 'max ops depth of nitrox (careful: MOD is a BASIC comand!)
'                       'Automatically sets nitrox stages in SUB stops
'
oxygenStops = 1         '%O2 for: 0 < stop depths <= 20 fsw in SUB stops
'                       'use: 1 = %100 for O2. use 0.21 = %21 for air deco
'
'----------------------------------------------------------------------------
'
'Here are constants to limit ascent. The supersaturation matrix is built below
'with the 1/2-time array.
'
'There is some variation in published numbers. This set of constants will give
'decompression times intermediate between the old USN and Buhlmann tables. No-stop times
'are similar to Buhlmann.
'----------------------------------------------------------------------------
'
g = 17.9              'nucleus skin tension (2-Dimensional) [dyne/cm]
gc = 257.0              'nucleus crushing pressure (2-Dimensional) [dyne/cm]
lambda = (7180 / 33)  '[ata*min] (related to critical volume). Yount and Hoffmann's lambda = 7500 was set in ffw*min
ro = 1.0                'critical radius at 1 ata [microns]
'
'
tugr = (2 * g / ro) * (1000 / 101325) 'Laplacian skin tension. To convert to [atm], use: dyne/(cm*micron) = 1000 Pa.
ggc = (g / gc)
'
'----------------------------------------------------------------------------
'Read in depth/run-time sequence
'----------------------------------------------------------------------------
'
DIM dt(2, num + 1)               'depth/run-time array
FOR l = 1 TO 2                   'initialize
  FOR i = 1 TO (num + 1)         'dt(1,i)'s are dive depths (levels)
    dt(l, i) = 0                 'dt(2,i)'s are run times at end of levels
  NEXT i
NEXT l
'
FOR i = 2 TO (num + 1)           'start time = 0 at depth = 0 (sea level)
50  PRINT USING "Enter Dive Level ##"; (i - 1)
    PRINT "Depth[in feet] , Run-time[in min]"
    INPUT ; dt(1, i), dt(2, i)
 
    IF dt(2, i) <= dt(2, i - 1) THEN 'trap: enter run times not stage times
       PRINT
       PRINT "Error in run-time. Redo last entry"
       PRINT
       GOTO 50
    ELSE
    END IF

    PRINT
NEXT i
'
'
DIM pet(2, num + 1)              'ambient pressure/run-time array
FOR l = 1 TO 2                   'initialize array
  FOR i = 1 TO (num + 1)         'pet(1,i)'s are pressures at each levels
    pet(l, i) = 0                'pet(2,i)'s are run-times at end of levels
  NEXT i
NEXT l
'
FOR i = 1 TO (num + 1)            'fill pressure/run-time array
    pet(1, i) = 1 + dt(1, i) / 33 'units are: atas (33 -> fsw)
    pet(2, i) = dt(2, i)
NEXT i
'
pdeep = 0                    'clear
'
DIM pressure(num)
FOR i = 2 TO (num + 1)       'starts at 2 to throw out surface value
 pressure(i - 1) = pet(1, i) 'build multilevel pressure array
NEXT i
CALL sort(pressure(), num)   'picks out largest pressure (deepest multi level dive depth)
                             'to set allowable supersaturation. max is returned by SUB Sort.
'
pdeep = max                  'pdeep is max dive pressure in [ata]
'
'----------------------------------------------------------------------------
'Compartment half-times and Supersaturation Gradients
'----------------------------------------------------------------------------
'
ncompartment = 16
'
psmin = (tugr * (1 - ggc) + ggc * (pdeep - 1))  '1st iteration allowable supersaturation [ata] to give Nsafe
'                                               'Yount and Hoffmann's eqs. 1 and 3

DIM pnw(ncompartment)         'matrix of allowed compartment supersaturations
FOR i = 1 TO ncompartment     'initialize supersaturation pressure allowed for ea compartment
 pnw(i) = psmin               'pnw is updated from the initial value after each iteration of the calculation
NEXT i
'
'
DIM halfT(ncompartment)         
halfT(1) = 4                   'Buhlmann's ZHL 16 half-times
halfT(2) = 8
halfT(3) = 12.5
halfT(4) = 18.5
halfT(5) = 27
halfT(6) = 38.3
halfT(7) = 64.3
halfT(8) = 77
halfT(9) = 109
halfT(10) = 146
halfT(11) = 187
halfT(12) = 239
halfT(13) = 305
halfT(14) = 390
halfT(15) = 498
halfT(16) = 635
'
DIM kN2(ncompartment)           'matrix of N2 time constants
'
FOR i = 1 TO ncompartment
 kN2(i) = (LOG(2) / halfT(i))
NEXT i
'
'----------------------------------------------------------------------------
' Dive level ppO2 reality check
'----------------------------------------------------------------------------
'
o2flag = 0
FOR i = 1 TO (num + 1)
 IF O2 * pet(1, i) >= 1.61 THEN o2flag = 1
NEXT i
'
IF o2flag = 1 THEN
 LOCATE 20, 20
 PRINT "WARNING: EXCESSIVE ppO2"
 LOCATE 21, 20
 PRINT "Press any key to continue"
 DO: LOOP WHILE INKEY$ = ""
ELSE
END IF
'
END SUB

SUB first
'***************************************************************
'                        
'This SUB selects the deepest stop (defining the first controlling
'compartment). The stop is rounded up to a multiple of 10 ft.
'
'********************************************************************
'                         'called by main
'
CALL sort(ceiling(), ncompartment)  'sort compartment ceilings
'

IF (max - INT(max)) > 0 THEN max = INT(max) + 1 'decimal max rounded up to a padded integer
dmax = max                          'max returned from SUB sort

N = dmax MOD 10
'
IF N = 0 THEN                       'depth of first stop used in SUB stops
    dfirst = dmax
ELSE
    dfirst = 10 * (INT(dmax / 10) + 1)
END IF
'
END SUB

FUNCTION gasExponential (timeconst, pfinal, pinitial, endtime, startime)
'*****************************************************************************
'                     
'Symmetric exponential in and of gassing
'
'*****************************************************************************
'
gasExponential = pfinal + (pinitial - pfinal) * (EXP(-timeconst * (endtime - startime)))
'
END FUNCTION

SUB mindy
'*****************************************************************************
'                           
'This sub calculates and stores the ascent ceilings for each compartment
'
'*****************************************************************************
'
'
FOR i = 1 TO ncompartment                     'supersaturation gradient pnw
 ceiling(i) = 33 * (tension(i) - pnw(i) - 1) 'set in each loop thru calc
NEXT i
'       
END SUB

SUB pnew
' **********************************************************************
'                           
'This sub sets new supersaturation gradients for all compartments on each
'iteration through the calculation of a decompression profile.
'(Yount and Hoffmanns's eqs. 4a,b,c solve for G in Wienke's equation D-12)
'
'Note that "fast compartments" allow large Gs because bubbles will not have much
'time to grow before perfusion\diffusion washes inert gas out of compartments, and allows
'the bubbles to begin shrinking by offgasing into the compartments
'
'Note that the larger the tdeco calculated on the previous iteration,
'the smaller the up-dated Gs are across compartments due to the need to limit
'the volume of gas evolved in all stages of deco.
'
'***********************************************************************
'
'                 'matrix pnw() is initially defined in dive data for
'                 'use in first pass through calc
DIM b(ncompartment)
DIM C(ncompartment)
'
FOR i = 1 TO ncompartment
   b(i) = 0
   C(i) = 0
NEXT i
'
FOR i = 1 TO ncompartment     'construct new values for supersaturation for ea compartment
'
   b(i) = psmin + (lambda * g) / (gc * (tdeco + halfT(i) / .693))               'Yount and Hoffmann's eq. 4b
   C(i) = ((g / gc) ^ 2) * (lambda * (pdeep - 1) / (tdeco + halfT(i) / .693))  'eq. 4c
   pnw(i) = (b(i) + SQR(b(i) ^ 2 - 4 * C(i))) / 2                          'eq. 4a
'
NEXT i
'
END SUB

SUB profile

'*****************************************************************************
' This sub outputs the dive and ascent profile
'*****************************************************************************
'
'
    stopNumber = 0
      FOR j = 1 TO 40
       IF decotable(2, j) > 0 THEN stopNumber = stopNumber + 1
      NEXT j
    
      DIM ascent(2, stopNumber)
      FOR i = 1 TO 2
       FOR j = 1 TO stopNumber
        ascent(i, j) = decotable(i, j)
       NEXT j
      NEXT i

IF stopNumber > 0 THEN
  
   IF outflag = 1 THEN               'store ascent stage depth, run-time, and %O2
      OPEN outfile$ FOR OUTPUT AS #1 'stair-step format suitable for ascii input to plotting program
      WRITE #1, 0, 0'store dive level time, depth
   
      FOR i = 2 TO (num + 1)
       WRITE #1, dt(2, (i - 1)), -dt(1, i)'store dive level time, depth
       WRITE #1, dt(2, i), -dt(1, i)
      NEXT i
   
      CLOSE #1
   
      OPEN outfile$ FOR APPEND AS #1
      WRITE #1, dt(2, num + 1), ascent(1, 1)
   
      FOR j = 1 TO stopNumber - 1
       WRITE #1, ascent(2, j), ascent(1, j)
       WRITE #1, ascent(2, j), ascent(1, j + 1)
      NEXT j
   
      WRITE #1, ascent(2, stopNumber), ascent(1, stopNumber)
      WRITE #1, ascent(2, stopNumber), 0
      CLOSE #1
 
    ELSE   'outflag =/= 1
    END IF

ELSE 'stopNumber = 0


END IF
'
PRINT
PRINT
PRINT
PRINT USING "TOTAL DECOMPRESSION TIME: #### min"; tdeco
IF (tdeco = 0) THEN PRINT " No Staged Decompression Required"
'
CALL runtime            'Display sequence of dive levels and run-times
'
END SUB

SUB putin
'*****************************************************************************
'
' This sub reads in a multilevel dive profile
'
'*****************************************************************************
'
outflag = 0   'Remove this line, and un-comment all code below for output of ascent
'                             
'                               'pick your own path
'outfile$ = "c:\myfile.txt"    'dive levels output in last iteration
'PRINT                         'ascent levels output in SUB stops
'PRINT "Should dive profile and ascent schedule be output to ",
'PRINT outfile$;
'PRINT " ?"
'PRINT "Note: If Yes, then old file will be overwriten"
'PRINT
'INPUT "Enter '1'for Yes or '0'for no. ", outflag
'
'IF outflag = 1 THEN
'   PRINT "Confirm Output"
'ELSEIF outflag = 0 THEN
'   PRINT "No Output"
'ELSE
'   outflag = 0
'   PRINT "No Output"
'END IF
'PRINT
'
INPUT "input number of dive levels: ", num
'
END SUB

SUB runtime
'*************************************************************
'                       
'This sub prints out the sequence of depths and run times
'
'**************************************************************
LOCATE 1, 1


PRINT "__________________________DIVE LEVELS___________________________"

 FOR i = 2 TO (num + 1)
     PRINT USING "dive level ## = ####.# ft   Run-time = ####.# min    fO2 = ###%"; (i - 1); dt(1, i); dt(2, i); 100 * O2
 NEXT i

PRINT "_____________________________ASCENT_____________________________"

END SUB

SUB sort (matrix(), length)
'*****************************************************************************
'                      
'This sub picks the largest number from an array and returns value max
'
'*****************************************************************************
max = 0
'
FOR i = 1 TO length
    number = matrix(i)
    IF number < 0 THEN                         'positive values for depth and pressure only
       number = 0
    ELSEIF number > max THEN
       max = number
    END IF
NEXT i
'
END SUB

SUB stops (k())
'*************************************************************************
'                   
'This sub calculates the time required at each stop and tracks compartment
'tensions.
'
'*************************************************************************
'
'
CLS
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT "                     DECOMPRESSION PROFILE"
'
'------------------------------------------------------------------------------
tdeco = 0                  'initialize total decompression time

FOR i = 1 TO 2
        FOR j = 1 TO 40
        decotable(i, j) = 0
        NEXT j
NEXT i
nstops = dfirst / 10       'number of stops. dfirst returned by from SUB first
'            
FOR q = 0 TO (nstops - 1)  'loop over stop depths to in or out gas compartments
'
'------------------------------------------------------------------------------
'
'  N2 partial pressures for each stage
'
'------------------------------------------------------------------------------
'
  d = (dfirst - q * 10)                                'depth of stop
  Pd = ((1 + d / 33) - 102 / 760)                      'ambient pressure
'                                                      'f is % N2
  IF d > maxod THEN f = N2                             'max operating depth of nitrox ascent stage
  IF d <= maxod AND d > 20 THEN f = (1 - nitroxStops)  'Ascent gases
  IF d <= 20 THEN f = (1 - oxygenStops)                'O2% set in SUB divedata
'
'------------------------------------------------------------------------------
  FOR j = 1 TO ncompartment 'loop over each compartment to in or out gas at stop depth d
'
    kay = kN2(j)
    tstop = 0               'clear stop tme
'
    IF q = 0 THEN           'If first stop, then
      PN2 = tension(j)      'read ppN2 at end of last dive level for each compartment
    ELSE                    'If 2nd stop etc, then read pN2 after previous stop for each compartment
      PN2 = stagetension(j)
    END IF
                             
    IF PN2 > (f * Pd) THEN  'If the compartment N2 tension is greater than ambient ppN2,  then outgas
                             
       arg = (pnw(j) + Pd * (1 - f) - 10 / 33) / (PN2 - f * Pd) 'pnw() is critical gradient from SUB divedata or SUB pnew
      
       IF arg > 0 AND arg < 1 THEN    'Condition argument of LOG
         tstop = -(1 / kay) * LOG(arg) 'find required time at stop for ea compartment
       ELSE
         tstop = 0          'no stop required by this compartment
       END IF

     ELSE                   'this compartment needs to continue ingasing
      GOTO 100              'so move to next compartment
     END IF

     degast(j) = tstop      'save stop time for each compartment
                            'last tstop zeroed at start of j loop
100 NEXT j                  'next compartment

  CALL sort(degast(), ncompartment)  'Find limiting compartment's degas time at stop
                        'Wienke's comment following his eq. A-12 should say: "Longest t_rem"
 
  IF (0 < max) AND (max < 1) THEN max = 1  'stops less than 1 min are padded
  IF (max - INT(max)) > 0 THEN max = INT(max) + 1 'decimal max rounded up to a padded integer
 
  stopTime = max            'limiting compartment degas time returned by SUB sort
 
  FOR j = 1 TO ncompartment  'clear stop times for next ascent stage
   degast(j) = 0
  NEXT j
'       
'------------------------------------------------------------------------------
'Now that the limiting compartment has been found,
'in or out gas all compartments for limiting compartment's outgas time
'------------------------------------------------------------------------------
 
  FOR j = 1 TO ncompartment    'each compartment either in or off gases at stop during time
   
    IF q = 0 THEN        'on first stop, find ending record from dive profile
       PN2 = tension(j)  'read ppN2 at end of last dive level for each compartment
    ELSE                 'If 2nd stop etc, then read ppN2 after previous stop for each compartment
       PN2 = stagetension(j)
    END IF               ' Now we have PN2

  kay = kN2(j)             'Here are the other arguments for the exponential in and out
  zero = 0               'gas excange start time is referenced to end of last stop
  ppAsymp = f * Pd       'Asymptotic N2 partial pressure
 
  P = gasExponential(kay, ppAsymp, PN2, stopTime, zero)'in or out gas all compartments for tmax
  stagetension(j) = P                      'update ppN2 for next stage

  NEXT j                                   'next compartment

  tdeco = tdeco + stopTime                     'running total of all stop times

  IF stopTime > 0 THEN PRINT USING "stop depth = ###.# ft   stop time = ###.# min   mix O2 = ### %"; d; stopTime; 100 * (1 - f)
 
  IF ((convflag = 1) OR (iteration1flag = 1)) AND (stopTime > 0) THEN  'if ascent schedule is converged
     decotable(1, (q + 1)) = -d            'make an array for output
     decotable(2, (q + 1)) = (tdeco + dt(2, num + 1))
  ELSEIF convflag = 0 THEN
  END IF
'
NEXT q                                     'next stop depth
'
END SUB