/*
 * Copyright (C) 2013-2015, 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 tissue compartments by continuous
 * functions was never implemented, but the name stuck anyway. Early versions
 * (numbered 0 to 9) were written in PARI/GP; standard C is used from 1.0 on.
 * 
 * Units are bars, liters, meters and minutes.
 *
 * WARNING: Any fitness of this software for actual dive planning is entirely
 * coincidental. This has never been the author's goal and he has taken no care
 * in this respect. You must be fully aware that diving safely requires proper
 * training and a thorough understanding of the limitations of your gear and
 * tools, including this piece of code.
 *
 * 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 99

/* 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%. Any gas with ox==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 */
} 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.01;

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

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

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

/* saltwater density */
static const double density = 1.03;

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

/* 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.6;
static const double gf2c = 0.9;
static const double gf1p = 4;
static const double gf2p = 1;

/* rough maximum operating depth */
static double mod (gas_t gas) {
	if (gas.sp) return INFTY;
	return 16 / gas.ox - 10;
}

/* 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.66) 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,
	};
	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;
}

/* holder for current gas */
static gas_t cur;

/* holder for current gas label */
static char str[MAX];

/* set gas label from gas_t */
static void gas_str (gas_t gas) {
	if (gas_same(cur,gas)) {
		*str = '\0';
		return;
	}
	cur = gas;
	int cc = 0;
	if (gas.sp) {
		sprintf(str, "SP=%04.2f,", gas.sp);
		cc = 8;
	}
	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 nitrogen in gas at given absolute pressure */
static double nitrogen (gas_t gas, double pressure) {
	if (gas.sp) return fmax(0, pressure - gas.sp) * (gas.ni / (gas.ni + gas.he));
	return gas.ni * pressure;
}

/* partial pressure of helium in gas at given absolute pressure */
static double helium (gas_t gas, double pressure) {
	if (gas.sp) return fmax(0, pressure - gas.sp) * (gas.he / (gas.ni + gas.he));
	return gas.he * pressure;
}

/* partial pressure of oxygen in gas at given absolute pressure */
static double oxygen (gas_t gas, double pressure) {
	if (gas.sp) return fmin(gas.sp, pressure);
	return gas.ox * pressure;
}


/* ****************************************************************************
 *
 * 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, atm);
		diver->he[i] = helium(air, atm);
	}
	for (i=0;i<MAX;i++) {
		diver->gas[i].ox = 0;
		diver->vol[i] = 0;
	}
}

/* update compartments after time breathing gas at pressure */
static void tension (diver_t *diver, double pressure, double time, gas_t gas) {
	int i;
	double ni = nitrogen(gas, pressure);
	double he = helium(gas, pressure);
	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) {
	int i;
	double p = pressure(depth);
	double ox = oxygen(gas, p);
	tension(diver, p, time, gas);
	diver->run += time;
	if (ox>0.5) {
		diver->cns += time / cns_limit(ox);
		diver->otu += time * pow(2 * ox - 1, 5 / 6);
	}
	if (!gas.sp) {
		for (i=0;i<MAX;i++) if (!diver->gas[i].ox || gas_same(gas,diver->gas[i])) break;
		diver->gas[i] = gas;
		diver->vol[i] += time * sac * p;
	}
}

/* gradient factor at given absolute ressure */
static double gf (double pressure) {
	double r = (gf1c * (gf2p - pressure) - gf2c * (gf1p - pressure)) / (gf2p - gf1p);
	return r<gf1c ? gf1c : r;
}

/* minimum tolerated absolute pressure */
static double ceiling (diver_t *diver) {
	int i;
	double r = -INFTY;
	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;
		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) {
	int i;
	double min = 0;
	double time = 0;
	double max = INFTY;
	while (max-min>step) {
		diver_t pig = *diver;
		tension(&pig, pressure(from), time, gas);
		if (ceiling(&pig)<pressure(to)) max = time; else min = 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);
}

/* decompression schedule segment from depth "from" to "to" */
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 = ceil(time2asc(diver, depth, next, gas));
		if (time>0) {
			expose(diver, depth, time, gas);
			printf("Stop at %2.0fm for %2.0f' (%3.0f)\n",depth,time,diver->run);
		}
		while (depth>next) {
			expose(diver, depth, step, gas);
			depth -= step * rate;
		}
		depth = next;
	}
}

/* full decompression schedule from depth using gases */
static void deco (diver_t *diver, double depth, gas_t *bailout) {
	do {
		double next = bailout[1].ox ? floor(mod(bailout[1])/3)*3 : 0;
		if (next>=depth) continue;
		gas_str(bailout[0]);
		printf("Up from %2.0fm         (%3.0f) %s\n", depth, diver->run, str);
		if (mod(bailout[0])<depth) printf("==> PPO2 WARNING\n");
		segment(diver, depth, next, bailout[0]);
		depth = next;
	}
	while ((++bailout)->ox);
}

/* print global dive volumes and toxicity levels */
static void summary (diver_t *diver) {
	gas_t *g = diver->gas;
	double *v = diver->vol;
	printf("--------------------------------\n");
	cur.ox = 0;
	while (g->ox) {
		gas_str(*g);
		printf("Used %4.0fL == %3.0fb of 12L %s\n", ceil(*v), ceil(*v / 12), str);
		g++;
		v++;
	}
	printf("CNS: %4.0f/100\n", diver->cns * 100);
	printf("OTU: %4.0f/300\n", diver->otu);
}

/* full decompression plan for bottom using bailout */
static void plan (segment_t *bottom, gas_t *bailout) {
	diver_t diver;
	afresh(&diver);
	double depth;
	while (bottom->gas.ox) {
		depth = bottom->depth;
		double time = bottom->time;
		gas_t gas = bottom->gas;
		expose(&diver, depth, time, gas);
		gas_str(gas);
		printf("Stay at %2.0fm for %2.0f' (%3.0f) %s\n", depth, time, diver.run, str);
		bottom++;
		if (bottom->gas.ox) segment(&diver, depth, bottom->depth, gas);
	}
	deco(&diver, depth, 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-2015, Gaetan Bisson. All rights reserved.\n\
Version 1.6; compiled "__DATE__".\n\
\n\
Usage: %s [gas] [depth:time] ... [bailout] ...\n\
"

int main (int argc, char *argv[]) {
	segment_t bottom[MAX];
	segment_t *s = bottom;
	gas_t bailout[MAX];
	gas_t *g = bailout;
	int b = 0;
	if (argc<3) {
		printf(ABOUT, argv[0]);
		return 1;
	}
	while ((++argv)[0]) {
		if (argv[0][0]>64) {
			*g = gas_parse(argv[0]);
			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;
}