aboutsummaryrefslogtreecommitdiffstatshomepage
path: root/deco.c
diff options
context:
space:
mode:
Diffstat (limited to 'deco.c')
-rw-r--r--deco.c279
1 files changed, 279 insertions, 0 deletions
diff --git a/deco.c b/deco.c
new file mode 100644
index 0000000..66182d9
--- /dev/null
+++ b/deco.c
@@ -0,0 +1,279 @@
+/* SPDX-License-Identifier: MIT-0 */
+
+#include <assert.h>
+#include <math.h>
+#include <stdbool.h>
+
+#include "deco.h"
+
+enum ALGO ALGO_VER = ZHL_16C;
+
+#define PO2_MAX (1.6)
+#define END_MAX (abs_depth(msw_to_bar(30)))
+#define RND(x) (round((x) *10000) / 10000)
+
+typedef struct zhl_n2_t {
+ double t;
+ double a[3];
+ double b;
+} zhl_n2_t;
+
+typedef struct zhl_he_t {
+ double t;
+ double a;
+ double b;
+} zhl_he_t;
+
+const zhl_n2_t ZHL16N[] = {
+ {.t = 5.0, .a = {1.1696, 1.1696, 1.1696}, .b = 0.5578},
+ {.t = 8.0, .a = {1.0000, 1.0000, 1.0000}, .b = 0.6514},
+ {.t = 12.5, .a = {0.8618, 0.8618, 0.8618}, .b = 0.7222},
+ {.t = 18.5, .a = {0.7562, 0.7562, 0.7562}, .b = 0.7825},
+ {.t = 27.0, .a = {0.6667, 0.6667, 0.6200}, .b = 0.8126},
+ {.t = 38.3, .a = {0.5933, 0.5600, 0.5043}, .b = 0.8434},
+ {.t = 54.3, .a = {0.5282, 0.4947, 0.4410}, .b = 0.8693},
+ {.t = 77.0, .a = {0.4701, 0.4500, 0.4000}, .b = 0.8910},
+ {.t = 109.0, .a = {0.4187, 0.4187, 0.3750}, .b = 0.9092},
+ {.t = 146.0, .a = {0.3798, 0.3798, 0.3500}, .b = 0.9222},
+ {.t = 187.0, .a = {0.3497, 0.3497, 0.3295}, .b = 0.9319},
+ {.t = 239.0, .a = {0.3223, 0.3223, 0.3065}, .b = 0.9403},
+ {.t = 305.0, .a = {0.2971, 0.2850, 0.2835}, .b = 0.9477},
+ {.t = 390.0, .a = {0.2737, 0.2737, 0.2610}, .b = 0.9544},
+ {.t = 498.0, .a = {0.2523, 0.2523, 0.2480}, .b = 0.9602},
+ {.t = 635.0, .a = {0.2327, 0.2327, 0.2327}, .b = 0.9653},
+};
+
+const zhl_he_t ZHL16He[] = {
+ {.t = 1.88, .a = 1.6189, .b = 0.4770},
+ {.t = 3.02, .a = 1.3830, .b = 0.5747},
+ {.t = 4.72, .a = 1.1919, .b = 0.6527},
+ {.t = 6.99, .a = 1.0458, .b = 0.7223},
+ {.t = 10.21, .a = 0.9220, .b = 0.7582},
+ {.t = 14.48, .a = 0.8205, .b = 0.7957},
+ {.t = 20.53, .a = 0.7305, .b = 0.8279},
+ {.t = 29.11, .a = 0.6502, .b = 0.8553},
+ {.t = 41.20, .a = 0.5950, .b = 0.8757},
+ {.t = 55.19, .a = 0.5545, .b = 0.8903},
+ {.t = 70.69, .a = 0.5333, .b = 0.8997},
+ {.t = 90.34, .a = 0.5189, .b = 0.9073},
+ {.t = 115.29, .a = 0.5181, .b = 0.9122},
+ {.t = 147.42, .a = 0.5176, .b = 0.9171},
+ {.t = 188.24, .a = 0.5172, .b = 0.9217},
+ {.t = 240.03, .a = 0.5119, .b = 0.9267},
+};
+
+double bar_to_msw(const double bar)
+{
+ return bar * 10;
+}
+
+double msw_to_bar(const double msw)
+{
+ return msw / 10;
+}
+
+double abs_depth(const double gd)
+{
+ return gd + SURFACE_PRESSURE;
+}
+
+double gauge_depth(const double ad)
+{
+ return ad - SURFACE_PRESSURE;
+}
+
+gas_t gas_new(const unsigned char o2, const unsigned char he, double mod)
+{
+ assert(o2 + he <= 100);
+
+ if (mod == MOD_AUTO) {
+ double mod_po2 = PO2_MAX / (o2 / 100.0);
+ double mod_end = END_MAX / (1 - he / 100.0);
+
+ mod = min(mod_po2, mod_end);
+ }
+
+ return (gas_t){.o2 = o2, .he = he, .n2 = 100 - o2 - he, .mod = mod};
+}
+
+int gas_equal(const gas_t *g1, const gas_t *g2)
+{
+ return g1->o2 == g2->o2 && g1->he == g2->he;
+}
+
+unsigned char gas_o2(const gas_t *gas)
+{
+ return gas->o2;
+}
+
+unsigned char gas_he(const gas_t *gas)
+{
+ return gas->he;
+}
+
+unsigned char gas_n2(const gas_t *gas)
+{
+ return gas->n2;
+}
+
+double gas_mod(const gas_t *gas)
+{
+ return gas->mod;
+}
+
+double add_segment_ascdec(decostate_t *ds, const double dstart, const double dend, const double time, const gas_t *gas)
+{
+ assert(time > 0);
+
+ const double rate = (dend - dstart) / time;
+
+ for (int i = 0; i < 16; i++) {
+ double pio = gas_he(gas) / 100.0 * (dstart - P_WV);
+ double po = ds->phe[i];
+ double r = gas_he(gas) / 100.0 * rate;
+ double k = log(2) / ZHL16He[i].t;
+ double t = time;
+
+ ds->phe[i] = pio + r * (t - 1 / k) - (pio - po - (r / k)) * exp(-k * t);
+ }
+
+ for (int i = 0; i < 16; i++) {
+ double pio = gas_n2(gas) / 100.0 * (dstart - P_WV);
+ double po = ds->pn2[i];
+ double r = gas_n2(gas) / 100.0 * rate;
+ double k = log(2) / ZHL16N[i].t;
+ double t = time;
+
+ ds->pn2[i] = pio + r * (t - 1 / k) - (pio - po - (r / k)) * exp(-k * t);
+ }
+
+ /* TODO add CNS */
+ /* TODO add OTU */
+
+ if (dend > ds->max_depth)
+ ds->max_depth = dend;
+
+ return time;
+}
+
+double add_segment_const(decostate_t *ds, const double depth, const double time, const gas_t *gas)
+{
+ assert(time > 0);
+
+ for (int i = 0; i < 16; i++) {
+ double pio = gas_he(gas) / 100.0 * (depth - P_WV);
+ double po = ds->phe[i];
+ double k = log(2) / ZHL16He[i].t;
+ double t = time;
+
+ ds->phe[i] = po + (pio - po) * (1 - exp(-k * t));
+ }
+
+ for (int i = 0; i < 16; i++) {
+ double pio = gas_n2(gas) / 100.0 * (depth - P_WV);
+ double po = ds->pn2[i];
+ double k = log(2) / ZHL16N[i].t;
+ double t = time;
+
+ ds->pn2[i] = po + (pio - po) * (1 - exp(-k * t));
+ }
+
+ /* TODO add CNS */
+ /* TODO add OTU */
+
+ if (depth > ds->max_depth)
+ ds->max_depth = depth;
+
+ return time;
+}
+
+double get_gf(const decostate_t *ds, const double depth)
+{
+ const unsigned char lo = ds->gflo;
+ const unsigned char hi = ds->gfhi;
+
+ if (ds->firststop == -1)
+ return lo;
+
+ if (depth < SURFACE_PRESSURE)
+ return hi;
+
+ if (depth > ds->firststop)
+ return lo;
+
+ /* interpolate lo and hi between first stop and surface */
+ double next_stop = depth - ds->ceil_multiple;
+ return hi - (hi - lo) * gauge_depth(next_stop) / gauge_depth(ds->firststop);
+}
+
+double round_ceiling(const decostate_t *ds, const double c)
+{
+ assert(ds->ceil_multiple != 0);
+
+ return abs_depth(ds->ceil_multiple * ceil(RND(gauge_depth(c) / ds->ceil_multiple)));
+}
+
+double ceiling(const decostate_t *ds, double gf)
+{
+ double c = 0;
+ gf /= 100;
+
+ for (int i = 0; i < 16; i++) {
+ /* n2 a and b values */
+ double an = ZHL16N[i].a[ALGO_VER];
+ double bn = ZHL16N[i].b;
+
+ /* he a and b values */
+ double ah = ZHL16He[i].a;
+ double bh = ZHL16He[i].b;
+
+ /* scale n2 and he values for a and b proportional to their pressure */
+ double pn2 = ds->pn2[i];
+ double phe = ds->phe[i];
+
+ double a = ((an * pn2) + (ah * phe)) / (pn2 + phe);
+ double b = ((bn * pn2) + (bh * phe)) / (pn2 + phe);
+
+ /* update ceiling */
+ c = max(c, ((pn2 + phe) - (a * gf)) / (gf / b + 1 - gf));
+ }
+
+ return round_ceiling(ds, c);
+}
+
+void init_tissues(decostate_t *ds)
+{
+ const double pn2 = 0.79 * (SURFACE_PRESSURE - P_WV);
+ const double phe = 0.00 * (SURFACE_PRESSURE - P_WV);
+
+ for (int i = 0; i < 16; i++)
+ ds->pn2[i] = pn2;
+
+ for (int i = 0; i < 16; i++)
+ ds->phe[i] = phe;
+}
+
+void init_decostate(decostate_t *ds, const unsigned char gflo, const unsigned char gfhi, const double ceil_multiple)
+{
+ init_tissues(ds);
+
+ ds->gflo = gflo;
+ ds->gfhi = gfhi;
+ ds->firststop = -1;
+ ds->ceil_multiple = ceil_multiple;
+}
+
+double ppO2(double depth, const gas_t *gas)
+{
+ return gas_o2(gas) / 100.0 * depth;
+}
+
+double end(double depth, const gas_t *gas)
+{
+ return (gas_o2(gas) + gas_n2(gas)) / 100.0 * depth;
+}
+
+double ead(double depth, const gas_t *gas)
+{
+ return depth * gas_n2(gas) / 79.0;
+}
generated by cgit v1.2.3-70-g09d2 (git 2.47.0) at 2024-12-24 16:56:48 +0000