/*
 * Copyright (C) 2013-2018, Gaetan Bisson <bisson@archlinux.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%. The maximum operating depth
 * is also specified and any gas with md==0 is considered invalid. */
typedef struct {
	double sp; /* closed circuit setpoint */
	double ox; /* fraction of oxygen */
	double ni; /* fraction of nitrogen */
	double he; /* fraction of helium */
	double md; /* maximum operating depth */
} 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 gas[MAX]; /* open circuit gases used */
	double vol[MAX]; /* corresponding volumes */
} diver_t;

/* A dive segment exposes the diver to a gas at a depth for a time. */
typedef struct {
	double depth;
	double time;
	gas_t gas;
} segment_t;


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


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

/* time step for discretization */
static const double step = 0.1;

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

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

/* atmospheric pressure */
static const double atm = 1.013;

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

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

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

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

/* inert gas compartments' half-time and m-value coefficients from ZH-L16B */
static const 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 };
static const 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 };
static const 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 };
static const 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 };
static const 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 };
static const 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 */
static const double gf1c = 0.30;
static const double gf2c = 0.85;
static const double gf1p = 7;
static const double gf2p = 1;

/* CNS toxicity exposure limit */
static 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 */
static const gas_t air = {
	.sp = 0,
	.ox = 0.21,
	.ni = 0.79,
	.he = 0,
};


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


/* test if two gases are the same */
static 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 */
static gas_t gas_parse (char *str) {
	gas_t gas = {
		.sp = 0,
		.ox = 1,
		.ni = 1,
		.he = 1,
		.md = 1,
	};
	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;
	}
	gas.md = gas.sp ? INFTY : 16 / gas.ox - 10;
	return gas;
}

/* set gas label from gas_t */
static 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 */
static 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 */
static 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 */
static 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 */
static void afresh (diver_t *diver) {
	int i;
	diver->run = 0;
	diver->cns = 0;
	diver->otu = 0;
	for (i=0;i<ZHL_N;i++) {
		diver->ni[i] = nitrogen(air, 0);
		diver->he[i] = helium(air, 0);
	}
	for (i=0;i<MAX;i++) {
		diver->gas[i].md = 0;
		diver->vol[i] = 0;
	}
}

/* update compartments after time breathing gas at depth */
static void tension (diver_t *diver, double depth, double time, gas_t gas) {
	int i;
	double ni = nitrogen(gas, depth);
	double he = helium(gas, depth);
	for (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 */
static 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) {
		int i;
		for (i=0;i<MAX;i++) if (!diver->gas[i].md || gas_same(gas,diver->gas[i])) break;
		diver->gas[i] = gas;
		diver->vol[i] += time * sac * pressure(depth);
	}
}

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

/* minimum tolerated ambiant pressure */
static double ceiling (diver_t *diver) {
	int i, j;
	double r = 0;
	for (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 (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 */
static 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 */
static double nofly (diver_t *diver) {
	return time2asc(diver, 0, -2.5, air);
}

/* print decompression schedule entry */
static void schedule (diver_t *diver, double depth, double time, gas_t gas) {
	int i;
	char buf[MAX+1];
	buf[MAX] = '\0';
	memset(buf, ' ', MAX);
	expose(diver, depth, time, gas);
	if (!gas_same(diver->gas[MAX-1],gas)) {
		sprintf(buf, "Up from %3.0fm            at %3.0f' on", depth, diver->run);
		gas_print(buf + 35, gas);
		diver->gas[MAX-1] = gas;
	} else if (!time) return;
	if (time) {
		int sec = ceil(time * 6);
		int min = floor(sec / 6);
		sec = (sec - min * 6) * 10;
		sprintf(buf, "Stay at %3.0fm for %3i:%02i to %3.0f'", depth, min, sec, diver->run);
	} 
	double ox = oxygen(gas, depth);
	if (depth>gas.md || ox>1.5 || ox<.18) sprintf(buf + 50, "PO2=%04.2f", ox);
	if (depth>gas.md || ox>1.7 || ox<.16) sprintf(buf + 59, "WARNING");
	for (i=0;i<MAX;i++) if (!buf[i]) buf[i]=' ';
	puts(buf);
}

/* decompression schedule segment from depth "from" to "to" on gas */
static void segment (diver_t *diver, double from, double to, gas_t gas) {
	double depth = from;
	while (depth>to) {
		double next = fmax(floor((depth - 1) / 3) * 3, to);
		double time = time2asc(diver, depth, next, gas);
		if (time) time = ceil(time + diver->run) - diver->run;
		schedule(diver, depth, time, gas);
		while (depth>next) {
			depth -= step * rate;
			expose(diver, depth, step, gas);
		}
		depth = next;
	}
}

/* print global dive volumes and toxicity levels */
static void summary (diver_t *diver) {
	gas_t *g = diver->gas;
	double *v = diver->vol;
	char str[MAX];
	printf("------------------------------------------\n");
	while (g->md) {
		gas_print(str, *g);
		printf("Consumed %4.0fL == %3.0fb from 12L of %-7s\n", ceil(*v), ceil(*v / 12), str);
		g++;
		v++;
	}
	printf("Load CNS: %3.0f/100\n", diver->cns * 100);
	printf("Load OTU: %3.0f/300\n", diver->otu);
}

/* full decompression schedule for bottom using bailout */
static void plan (segment_t *bottom, gas_t *bailout) {
	double depth;
	diver_t diver;
	afresh(&diver);
	while (bottom->gas.md) {
		double time = bottom->time;
		gas_t gas = bottom->gas;
		depth = bottom->depth;
		schedule(&diver, depth, time, gas);
		bottom++;
		if (bottom->gas.md) segment(&diver, depth, bottom->depth, gas);
	}
	while (bailout->md) {
		double next = floor(bailout[1].md / 3) * 3;
		if (depth>next) {
			segment(&diver, depth, next, bailout[0]);
			depth = next;
		}
		bailout++;
	}
	summary(&diver);
}


/* ****************************************************************************
 *
 * MAIN
 *
 * Parse arguments into bottom/gases arrays; compute deco plan.
 */


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

int main (int argc, char *argv[]) {
	segment_t bottom[MAX];
	gas_t bailout[MAX];
	memset(bottom, 0, sizeof bottom);
	memset(bailout, 0, sizeof bailout);
	segment_t *s = bottom;
	gas_t *g = bailout;
	int b = 0;
	if (argc<3) {
		printf(ABOUT, argv[0]);
		return 1;
	}
	while ((++argv)[0]) {
		if (argv[0][0]>64) {
			char *t = argv[0]; while (*t && *t!='@') t++;
			double m = *t ? atof(t+1) : 0; *t = '\0';
			*g = gas_parse(argv[0]);
			if (m) g->md = m;
			if (!argv[1] || argv[1][0]>64) b = 1;
			if (b) g++;
		} else {
			char *t = argv[0]; while (*t!=':') t++; *t = '\0';
			s->depth = atof(argv[0]);
			s->time = atof(t+1);
			s->gas = *g;
			s++;
		}
	}
	plan(bottom, bailout);
	return 0;
}