/*
 * 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.
 *
 * WARNING: This software is unfit for actual dive planning! Diving safely
 * requires proper training and a thorough understanding of the limitations of
 * your gear and tools, including this piece of code.
 * 
 * Units are bars, liters, meters and minutes.
 *
 * Compile with:
 *
 *   cc -O2 -lm -o decox decox.c
 */


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


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


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


/* infinity for programmers */
#define MAX 99

/* 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;
}

/* nitrogen compartments' half-time and m-value from ZH-L16B */
static const double zh_ht[] = {5, 8, 12.5, 18.5, 27, 38.3, 54.3, 77, 109, 146, 187, 239, 305, 390, 498, 635};
static const double zh_m0[] = {29.6, 25.4, 22.5, 20.3, 19.0, 17.5, 16.5, 15.7, 15.2, 14.6, 14.2, 13.9, 13.4, 13.2, 12.9, 12.7};
static const double zh_md[] = {1.7928, 1.5352, 1.3847, 1.2780, 1.2306, 1.1857, 1.1504, 1.1223, 1.0999, 1.0844, 1.0731, 1.0635, 1.0552, 1.0478, 1.0414, 1.0359};
#define ZH_N 16

/* gradient factor coefficients and pressure points */
static const double gf1c = 0.9;
static const double gf2c = 0.3;
static const double gf1p = 1;
static const double gf2p = 7;

/* gradient factor at ambiant ressure */
static double gf (double pressure) {
	return (gf1c * (gf2p - pressure) - gf2c*(gf1p - pressure)) / (gf2p - gf1p);
}

/* m-value from ambiant pressure */
static double mval (double pressure, int i) {
	return pressure + gf(pressure) * (zh_m0[i] / 10 + zh_md[i] * (pressure - 1) - pressure);
}

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


/* ****************************************************************************
 *
 * BREATHABLE GASES
 *
 */

#define AIR 128
#define EAN 256
#define CCR 512

/* convert from gas label to internal representation */
static double parse_gas (char *gas) {
	if (!strncmp(gas,"SP=",3)) return CCR + atof(gas + 3);
	if (!strncmp(gas,"EAN",3)) return EAN + atof(gas + 3);
	if (!strncmp(gas,"AIR",3)) return AIR;
	return 0;
}

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

/* set gas label from internal representation */
static void str_gas (double gas) {
	if (gas>=CCR) {
		double sp = gas - CCR;
		sprintf(str, "SP=%04.2f", sp);
	} else
	if (gas>=EAN) {
		double o2 = gas - EAN;
		sprintf(str, "EAN%02.0f", o2);
	} else
	if (gas>=AIR) {
		sprintf(str, "AIR");
	}
}

/* partial pressure of nitrogen in gas at given absolute pressure */
static double nitrogen (double gas, double pressure) {
	if (gas>=CCR) {
		double sp = gas - CCR;
		return fmax(0, pressure - sp);
	}
	if (gas>=EAN) {
		double o2 = gas - EAN;
		return (1 - o2 / 100) * pressure;
	}
	if (gas>=AIR) {
		return 0.79 * pressure;
	}
	return INFTY;
}

/* partial pressure of oxygen in gas at given absolute pressure */
static double oxygen (double gas, double pressure) {
	if (gas>=CCR) {
		double sp = gas - CCR;
		return fmin(sp, pressure);
	}
	if (gas>=EAN) {
		double o2 = gas - EAN;
		return (o2 / 100) * pressure;
	}
	if (gas>=AIR) {
		return 0.21 * pressure;
	}
	return INFTY;
}

/* maximum operation depth of gas */
static double mod (double gas) {
	if (gas>=CCR) {
		return INFTY;
	}
	if (gas>=EAN) {
		double o2 = gas - EAN;
		return 1600 / o2 - 10;
	}
	if (gas>=AIR) {
		return 66;
	}
	return 0;
}

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


typedef struct {
	double run;
	double comp[ZH_N];
	double vols[MAX];
	double cns;
	double otu;
} diver_t;

/* freshly rested diver */
static void afresh (diver_t *diver) {
	int i;
	diver->cns = 0;
	diver->otu = 0;
	diver->run = 0;
	for (i=0;i<ZH_N;i++) diver->comp[i] = nitrogen(AIR, atm);
	for (i=0;i<MAX;i++) diver->vols[i] = 0;
}

/* compartments' tension after time at pn2 */
static void tension (double *comp, double time, double pn2) {
	int i;
	for (i=0;i<ZH_N;i++) comp[i] = pn2 + (comp[i] - pn2) * pow(0.5, time / zh_ht[i]);
}

/* expose diver to depth for time breahting gas */
static void expose (diver_t *diver, double depth, double time, double gas) {
	double p = pressure(depth);
	double pn2 = nitrogen(gas, p);
	double po2 = oxygen(gas, p);
	diver->run += time;
	tension(diver->comp, time, pn2);
	if (po2>0.5) {
		diver->cns += time / cns_limit(po2);
		diver->otu += time * pow(2 * po2 - 1, 5 / 6);
	}
	int i;
	if (gas<CCR) {
		for (i=0;i<MAX;i+=2) if (!diver->vols[i] || diver->vols[i]==gas) break;
		diver->vols[i] = gas;
		diver->vols[i+1] += time * sac * pressure(depth);
	}
}

/* time breathing pn2 until any compartment greater than lim */
static double time2any (double pn2, double *comp, double *lim) {
	double time = INFTY;
	double q;
	int i;
	for (i=0;i<ZH_N;i++) {
		q = (pn2 - comp[i]) / (pn2 - lim[i]);
		if (q>0) {
			time = fmin(time, zh_ht[i] * log(q) / log(2));
		} else if (pn2>lim[i]) {
			time = 0;
		}
	}
	return time;
}

/* time breathing pn2 until all compartments smaller than lim */
static double time2all (double pn2, double *comp, double *lim) {
	double time = 0;
	double q;
	int i;
	for (i=0;i<ZH_N;i++) {
		q = (pn2 - comp[i]) / (pn2 - lim[i]);
		if (q>0) {
			time = fmax(time, zh_ht[i] * log(q) / log(2));
		} else if (pn2>lim[i]) {
			time = INFTY;
		}
	}
	return time;
}

/* time at sea level before 0.7 atm allowable */
static double nofly (double *comp) {
	int i;
	double cabin[ZH_N];
	for (i=0;i<ZH_N;i++) cabin[i] = mval(0.7 * atm, i);
	return time2all(nitrogen(AIR, atm), comp, cabin);
}

/* time at sea level before complete desaturation */
static double desat (double *comp) {
	int i;
	double fresh[ZH_N];
	for (i=0;i<ZH_N;i++) fresh[i] = 1.01 * nitrogen(AIR, atm);
	return time2all(nitrogen(AIR, atm), comp, fresh);
}

/* time at depth "from" until direct ascent to "to" allowable */
static double time2asc (double *comp, double from, double to, double gas) {
	int i;
	double depth = to;
	double pn2;
	double p;
	double worst[ZH_N];
	for (i=0;i<ZH_N;i++) worst[i] = INFTY;
	while (depth<from) {
		p = pressure(depth);
		pn2 = nitrogen(gas, p);
		for (i=0;i<ZH_N;i++) worst[i] = fmin(worst[i], mval(pressure(depth), i));
		tension(worst, -step, pn2);
		depth += step * rate;
	}
	p = pressure(from);
	double stop = time2all(nitrogen(gas, p), comp, worst);
	double stay = time2any(nitrogen(gas, p), comp, worst);
	if (stop>0) return stop;
	if (stay>0) return -stay;
	return -INFTY;
}

/* decompression schedule segment from depth "from" to "to" */
static void segment (diver_t *diver, double from, double to, double gas) {
	double depth = from;
	double next;
	double time;
	while (depth>to) {
		next = fmax(floor((depth - 1) / 3) * 3, to);
		time = ceil(time2asc(diver->comp, 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;
	}
}


/* no decompression limit breathing gas at depth */
static double ndl (double gas, double depth) {
	int i;
	double fresh[ZH_N];
	for (i=0;i<ZH_N;i++) fresh[i] = nitrogen(AIR, atm);
	return fmax(0, -time2asc(fresh, depth, 0, gas));
}

/* print global dive volumes and toxicity levels */
static void summary (diver_t *diver) {
	double *v = diver->vols;
	printf("---------------------------------\n");
	while (*v) {
		str_gas(v[0]);
		printf("Used %4.0fL == %3.0fb of 12L %s\n", ceil(v[1]), ceil(v[1] / 12), str);
		v += 2;
	}
	printf("CNS: %4.0f/100\n", diver->cns * 100);
	printf("OTU: %4.0f/300\n", diver->otu);
}

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

/* full decompression plan for bottom using gases */
static void plan (double *bottom, double *gases) {
	diver_t diver;
	double time, depth, gas;
	double last = 0;
	double next;
	afresh(&diver);
	while (*bottom) {
		depth = bottom[0];
		time = bottom[1];
		gas = bottom[2];
		if (last) segment(&diver, last, depth, gas);
		expose(&diver, depth, time, gas);
		str_gas(gas);
		printf("Stay at %2.0fm for %2.0f' (%3.0f) %s\n", depth, time, diver.run, str);
		last = depth;
		bottom += 3;
	}
	deco(&diver, depth, gases);
}


/* ****************************************************************************
 *
 * 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.1; compiled "__DATE__".\n\
\n\
Usage: %s [depth time gas] ... [deco_gas] ...\n\
"

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

	int i;
	double bottom[MAX];
	double gases[MAX];

	i = 0;
	while (*argv && *argv[0]<64) {
		bottom[i+0] = atof(argv[0]);
		bottom[i+1] = atof(argv[1]);
		bottom[i+2] = parse_gas(argv[2]);
		argv += 3;
		i += 3;
	}
	bottom[i] = 0;

	i = 0;
	while (*argv) {
		gases[i] = parse_gas(*argv);
		argv++;
		i++;
	}
	gases[i] = 0;

	plan(bottom, gases);
}