/*
 * Copyright (C) 2013-2022, Gaetan Bisson <gaetan@fenua.org>.
 *
 * Permission to use, copy, modify, and/or distribute this software for any
 * purpose with or without fee is hereby granted, provided that the above
 * copyright notice and this permission notice appear in all copies.
 *
 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
 * SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION
 * OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
 * CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
 */

/*
 * Decox. Decreasing continuous exponentials; toy decompression planner.
 *
 * Decox merely aims to give its author a better understanding of the behavior
 * of decreasing exponentials in the context of Haldanian decompression models.
 * The original idea to eventually model gas loading over a continuous range of
 * tissue compartments was never implemented, but the name stuck anyway. Early
 * versions (0 to 9) were written in PARI/GP; standard C is used from 1.0 on.
 *
 * Units are bars, liters, meters and minutes.
 *
 * WARNING: This software is not intended to facilitate actual dive planning.
 * Using it safely requires, as for any other dive tool, proper training and
 * full awareness of its limitations. USE AT YOUR OWN RISK.
 *
 * Compile with:
 *
 *   cc -O2 -lm -o decox decox.c
 */


/* ****************************************************************************
 *
 * HEADERS
 *
 */


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


/* ****************************************************************************
 *
 * FANCY STRUCTURES
 *
 */


/* We use static arrays of predefined lengths; increase if needed. */
#define ZHL_N 16
#define MAX 70

/* We store gases as a closed circuit setpoint (zero for open circuit mixes)
 * and fractions of oxygen, nitrogen, and helium; though this is not enforced,
 * those fractions are expected to add up to 100%. Open circuit volume use is
 * also tracked and set to zero for disabled gases.
 */
typedef struct {
	double sp; /* closed circuit setpoint */
	double ox; /* fraction of oxygen */
	double ni; /* fraction of nitrogen */
	double he; /* fraction of helium */
	double vl; /* open circuit volume used */
} gas_t;

/* We store a diver's runtime through the dive, oxygen toxicy exposure levels
 * (CNS and OTU), tensions of compartments for inert gas (nitrogen and helium),
 * and open circuit gas consumption. */
typedef struct {
	double run; /* dive runtime */
	double cns; /* CNS fraction */
	double otu; /* OTU level */
	double ni[ZHL_N]; /* compartments nitrogen tension */
	double he[ZHL_N]; /* compartments helium tension */
	gas_t gases[MAX]; /* open circuit gases used */
} diver_t;


/* ****************************************************************************
 *
 * GLOBAL PHYSICAL CONSTANTS
 *
 */


/* infinity for physicists */
double INFTY = 999999999;

/* time step for discretization */
double step = 0.003;

/* surface air consumption */
double sac = 15;

/* safe ascent rate */
double rate = 9;

/* atmospheric pressure */
double atm = 1.013;

/* saltwater density in kilograms per liter */
double density = 1.03;

/* weight of one kilogram in newtons */
double weight = 9.81;

/* ambient pressure at depth */
double pressure (double depth) {
	return atm + weight * density * depth / 100;
}

/* alveolar water vapor pressure */
double vapor = 0.0627;

/* inert gas compartments' half-time and m-value coefficients from ZH-L16B */
double zhl_ni_h[] = {    5.0,    8.0,   12.5,   18.5,   27.0,   38.3,   54.3,   77.0,  109.0,  146.0,  187.0,  239.0,  305.0,  390.0,  498.0,  635.0 };
double zhl_ni_a[] = { 1.1696, 1.0000, 0.8618, 0.7562, 0.6667, 0.5600, 0.4947, 0.4500, 0.4187, 0.3798, 0.3497, 0.3223, 0.2850, 0.2737, 0.2523, 0.2327 };
double zhl_ni_b[] = { 0.5578, 0.6514, 0.7222, 0.7825, 0.8126, 0.8434, 0.8693, 0.8910, 0.9092, 0.9222, 0.9319, 0.9403, 0.9477, 0.9544, 0.9602, 0.9653 };
double zhl_he_h[] = {   1.88,   3.02,   4.72,   6.99,  10.21,  14.48,  20.53,  29.11,  41.20,  55.19,  70.69,  90.34, 115.29, 147.42, 188.24, 240.03 };
double zhl_he_a[] = { 1.6189, 1.3830, 1.1919, 1.0458, 0.9220, 0.8205, 0.7305, 0.6502, 0.5950, 0.5545, 0.5333, 0.5189, 0.5181, 0.5176, 0.5172, 0.5119 };
double zhl_he_b[] = { 0.4770, 0.5747, 0.6527, 0.7223, 0.7582, 0.7957, 0.8279, 0.8553, 0.8757, 0.8903, 0.8997, 0.9073, 0.9122, 0.9171, 0.9217, 0.9267 };

/* gradient factor coefficients and pressure points */
double gf1c = 0.30;
double gf2c = 0.85;
double gf1p = 7;
double gf2p = 1;

/* CNS toxicity exposure limit */
double cns_limit (double ox) {
	if (ox<0.5) return INFTY;
	if (ox<0.6) return 1800 - 1800 * ox;
	if (ox<0.7) return 1620 - 1500 * ox;
	if (ox<0.8) return 1410 - 1200 * ox;
	if (ox<0.9) return 1170 - 900 * ox;
	if (ox<1.1) return 900 - 600 * ox;
	if (ox<1.5) return 570 - 300 * ox;
	if (ox<1.65) return 1245 - 750 * ox;
	return 1;
}

/* run-of-the-mill air */
gas_t air = {
	.sp = 0,
	.ox = 0.21,
	.ni = 0.79,
	.he = 0,
	.vl = 1,
};

/* invalid gas as default */
gas_t invalid = {
	.sp = 0,
	.ox = 0.01,
	.ni = 0.99,
	.he = 0.99,
	.vl = 1,
};

/* minimum segment duration for gas switch */
double gstime = 0.3;

/* available tank sizes and working pressure */
double tanks[] = { 6, 12, 15, 0 };
double wpress = 210;


/* ****************************************************************************
 *
 * WORKING WITH GASES
 *
 */


/* disable all gases from array */
void gas_disable (gas_t *gases) {
	for (int i=0;i<MAX;i++) gases[i].vl = 0;
}

/* test if two gases are the same */
int gas_same (gas_t a, gas_t b) {
	return a.sp==b.sp && a.ox==b.ox && a.ni==b.ni && a.he==b.he;
}

/* convert from gas label to gas_t */
gas_t gas_parse (char *str) {
	gas_t gas = invalid;
	if (!strncasecmp(str,"SP=",3)) {
		char *t = str; while (*t!=',') t++; *t = '\0';
		gas.sp = atof(str + 3);
		str = t + 1;
	}
	if (!strncasecmp(str,"AIR",3)) {
		gas.ox = 0.21;
		gas.ni = 0.79;
		gas.he = 0;
	}
	if (!strncasecmp(str,"O2",2)) {
		gas.ox = 1;
		gas.ni = 0;
		gas.he = 0;
	}
	if (!strncasecmp(str,"Nx",2)) {
		gas.ox = atof(str + 2) / 100;
		gas.ni = 1 - gas.ox;
		gas.he = 0;
	}
	if (!strncasecmp(str,"Hx",2)) {
		gas.ox = atof(str + 2) / 100;
		gas.he = 1 - gas.ox;
		gas.ni = 0;
	}
	if (!strncasecmp(str,"Tx",2)) {
		char *t = str; while (*t!='/') t++; *t = '\0';
		gas.ox = atof(str + 2) / 100;
		gas.he = atof(t + 1) / 100;
		gas.ni = 1 - gas.ox - gas.he;
	}
	return gas;
}

/* set gas label from gas_t */
void gas_print (char *str, gas_t gas) {
	int cc = 0;
	if (gas.sp) {
		sprintf(str, "SP=%03.1f,", gas.sp);
		cc = 7;
	}
	if (gas.he) {
		if (!gas.ni) sprintf(str + cc, "Hx%02.0f", gas.ox * 100);
		else sprintf(str + cc, "Tx%02.0f/%02.0f", gas.ox * 100, gas.he * 100);
	} else {
		if (!gas.ni) sprintf(str + cc, "O2");
		else if (gas.ni==0.79) sprintf(str + cc, "AIR");
		else sprintf(str + cc, "Nx%02.0f", gas.ox * 100);
	}
}

/* partial pressure of oxygen in gas at depth */
double oxygen (gas_t gas, double depth) {
	double p = pressure(depth) - vapor;
	return fmin(fmax(gas.sp, gas.ox * p), p);
}

/* partial pressure of nitrogen in gas at depth */
double nitrogen (gas_t gas, double depth) {
	if (!gas.ni) return 0;
	double p = pressure(depth) - vapor;
	double q = p - oxygen(gas, depth);
	return q * gas.ni / (gas.ni + gas.he);
}

/* partial pressure of helium in gas at depth */
double helium (gas_t gas, double depth) {
	if (!gas.he) return 0;
	double p = pressure(depth) - vapor;
	double q = p - oxygen(gas, depth);
	return q * gas.he / (gas.ni + gas.he);
}


/* ****************************************************************************
 *
 * DECOMPRESSION FUNCTIONS AND CONSTANTS
 *
 */


/* freshly rested diver */
void afresh (diver_t *diver) {
	diver->run = 0;
	diver->cns = 0;
	diver->otu = 0;
	for (int i=0;i<ZHL_N;i++) {
		diver->ni[i] = nitrogen(air, 0);
		diver->he[i] = helium(air, 0);
	}
	gas_disable(diver->gases);
}

/* update compartments after time breathing gas at depth */
void tension (diver_t *diver, double depth, double time, gas_t gas) {
	double ni = nitrogen(gas, depth);
	double he = helium(gas, depth);
	for (int i=0;i<ZHL_N;i++) {
		diver->ni[i] = ni + (diver->ni[i] - ni) * pow(0.5, time / zhl_ni_h[i]);
		diver->he[i] = he + (diver->he[i] - he) * pow(0.5, time / zhl_he_h[i]);
	}
}

/* expose diver to depth for time breahting gas */
void expose (diver_t *diver, double depth, double time, gas_t gas) {
	diver->run += time;
	tension(diver, depth, time, gas);
	double ox = oxygen(gas, depth);
	if (ox>0.5) {
		diver->cns += time / cns_limit(ox);
		diver->otu += time * pow(2 * ox - 1, 0.8333);
	}
	if (!gas.sp) {
		gas_t *g = diver->gases;
		while (g->vl && !gas_same(gas,*g)) g++;
		if (!g->vl) *g = gas;
		g->vl += time * sac * pressure(depth);
	}
}

/* gradient factor at given ambiant ressure */
double gf (double pressure) {
	double r = (gf1c * (gf2p - pressure) - gf2c * (gf1p - pressure)) / (gf2p - gf1p);
	return fmin(fmax(r, gf1c), gf2c);
}

/* minimum tolerated ambiant pressure */
double ceiling (diver_t *diver) {
	double r = 0;
	for (int i=0;i<ZHL_N;i++) {
		double c = diver->ni[i] + diver->he[i];
		double a = (diver->ni[i] * zhl_ni_a[i] + diver->he[i] * zhl_he_a[i]) / c;
		double b = (diver->ni[i] * zhl_ni_b[i] + diver->he[i] * zhl_he_b[i]) / c;
		double t = (c - a) * b;
		for (int j=0;j<3;j++) {
			double g = gf(t);
			t = (c - g * a) / (g / b - g + 1);
		}
		if (t>r) r = t;
	}
	return r;
}

/* time at depth "from" until direct ascent to "to" allowable */
double time2asc (diver_t *diver, double from, double to, gas_t gas) {
	double min = 0;
	double time = 0;
	double max = INFTY;
	while (max-min>step) {
		double depth = from;
		diver_t pig = *diver;
		tension(&pig, depth, time, gas);
		while (depth>to) {
			depth -= step * rate;
			if (ceiling(&pig)>pressure(depth)) break;
			tension(&pig, depth, step, gas);
		}
		if (depth>to) min = time; else max = time;
		time = (min + max) / 2;
	}
	return time;
}

/* time at sea level before 0.75 atm allowable */
double nofly (diver_t *diver) {
	return time2asc(diver, 0, -2.5, air);
}

/* print runtime headers */
void banner() {
	printf("Depth  Gas             Stop   Run   PO2\n");
	printf("-----  --------------  -----  ----  ----\n");
}

#define LASTGAS diver->gases[MAX-1]

/* print dive segment */
void segment (diver_t *diver, double depth, double time, gas_t gas) {
	char str[MAX];
	if (!gas_same(gas, LASTGAS)) {
		gas_print(str, gas);
		LASTGAS = gas;
	} else {
		strcpy(str,"--");
	}
	expose(diver, depth, time, gas);

	int sec = ceil(time * 6);
	int min = floor(sec / 6);
	sec = (sec - min * 6) * 10;
	printf("%4.0fm  %-14s  %2i:%02i  %3.0f'", depth, str, min, sec, diver->run);

	double ox = oxygen(gas, depth);
	if (ox>1.4 || ox<.18) printf("  %04.2f", ox);
	printf("\n");
}

/* decompression schedule from depth "from" to "to" using gases */
void schedule (diver_t *diver, double from, double to, gas_t *gas) {
	double depth = from;
	while (depth>to) {
		while (gas[1].vl && oxygen(gas[1],depth)<1.6) gas++;
		double next = fmax(floor((depth - 1) / 3) * 3, to);
		double time = time2asc(diver, depth, next, *gas);
		if (!gas_same(*gas, LASTGAS)) time = fmax(time, gstime);
		if (time) {
			/* round up so runtime matches printed value */
			time = ceil(diver->run + time) - diver->run;
			segment(diver, depth, time, *gas);
		}
		while (depth>next) {
			expose(diver, depth, step, *gas);
			depth -= step * rate;
		}
		/* depth overshoot may decrease effective rate */
		depth = next;
	}
}

/* print gas info and toxicity levels */
void summary (diver_t *diver) {
	printf("\n");
	printf("OC Gas   MOD   Volume  Tank  Fill\n");
	printf("-------  ----  ------  ----  ----\n");
	gas_t *gas = diver->gases-1;
	while ((++gas)->vl) {
		char str[MAX];
		gas_print(str, *gas);

		double mod = floor(16 / gas->ox - 10);
		double volume = floor(gas->vl);

		/* find fitting tank */
		double *size = tanks;
		while (*size && volume > *size*wpress) size++;
		if (!*size) size--; /* use last even if too small */
		double fill = ceil(volume / *size);

		printf("%-7s  %3.0fm  %5.0fL  %3.0fL  %3.0fb\n", str, mod, volume, *size, fill);
	}

	printf("\n");
	printf("CNS: %3.0f/100\n", diver->cns * 100);
	printf("OTU: %3.0f/300\n", diver->otu);
}


/* ****************************************************************************
 *
 * MAIN
 *
 * Parse arguments and compute deco plan.
 */


#define ABOUT "\
Decox. Decreasing continuous exponentials; toy decompression planner.\n\
Copyright (C) 2013-2022, Gaetan Bisson. All rights reserved.\n\
Version 2.6; compiled "__DATE__".\n\
\n\
Usage: %s gas depth:time ...\n\
"

void main (int argc, char *argv[]) {
	if (argc<3) {
		printf(ABOUT, argv[0]);
		return;
	}

	gas_t buf[MAX];
	gas_disable(buf);
	buf[0] = invalid;
	gas_t *gases, *gas = buf+1;
	double last = -INFTY;

	diver_t diver;
	afresh(&diver);
	banner();

	while (*(++argv)) {
		if (**argv>64) *(gas++) = gas_parse(*argv);
		if (**argv<64) {
			char *t = *argv; while (*t!=':') t++; *t = '\0';
			double depth = atof(*argv);
			double time = atof(t+1);
			if (last>depth) schedule(&diver, last, depth, gases);
			gases = gas-1;
			last = depth;
			segment(&diver, depth, time, *gases);
		}
	}
	schedule(&diver, last, 0, gases);
	summary(&diver);
}