/*
utils.c

utilities
*/

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


// nearest id : p[n1]...p[n2]
int nearest(double x, int n1, int n2, const double *p)
{
	int    imin = n1;
	double dmin = fabs(x - p[n1]);
	for (int i = n1 + 1; i <= n2; i++) {
		double d = fabs(x - p[i]);
		if (d < dmin) {
			dmin = d;
			imin = i;
		}
	}

	return imin;
}


// tokenize a string
int tokenize(char *str, const char tokensep[], char *token[], size_t maxtoken)
{
	if ((str == NULL) || !maxtoken) return 0;

	char *thistoken = strtok(str, tokensep);

	int   count;
	for (count = 0; (count < maxtoken) && (thistoken != NULL); ) {
		token[count++] = thistoken;
		thistoken = strtok(NULL, tokensep);
	}

	token[count] = NULL;

	return count;
}


// p1 <= p[n1] <= p[n2] <= p2
// p : p[0] ... p[n - 1]
void getspan(double *p, int n, double p1, double p2, int *n1, int *n2, double eps)
{
	*n1 = 0;
	*n2 = n - 1;

	if (n < 1) return;

	// q1 <= q2
	const double q1 = (p1 < p2) ? p1 : p2;
	const double q2 = (p1 > p2) ? p1 : p2;

	if      ((q1 < p[0] - eps) && (q2 < p[0] - eps)) {
		// p1, p2 < p[0] -> n1 > n2
		*n1 = -1;
		*n2 = -2;
	}
	else if ((q1 > p[n - 1] + eps) && (q2 > p[n - 1] + eps)) {
		// p1, p2 > p[n -1] -> n1 > n2
		*n1 = n + 1;
		*n2 = n;
	}
	else {
		// p1 <= p[n1]
		for (int i = 0; i < n; i++) {
			if (p[i] > q1 - eps) {
				*n1 = i;
				break;
			}
		}

		// p[n2] <= p2
		for (int i = n - 1; i >= 0; i--) {
			if (p[i] < q2 + eps) {
				*n2 = i;
				break;
			}
		}
	}
}


// interpolation of the table using first column
void interpolation(double x, int ntable, int ndata, double **table, double f[])
{
	if (ntable < 2) return;
	if (ndata < 2) return;

	int    n = 0;
	double a = 0;
	if (x <= table[0][0]) {
		n = 0;
		a = 0;
	}
	else if (x >= table[ntable - 1][0]) {
		n = ntable - 2;
		a = 1;
	}
	else {
		for (int i = 0; i < ntable - 1; i++) {
			if ((table[i][0] <= x) && (x <= table[i + 1][0])) {
				n = i;
				a = (x - table[i][0]) / (table[i + 1][0] - table[i][0]);
				break;
			}
		}
	}

	for (int m = 1; m < ndata; m++) {
		f[m - 1] = ((1 - a) * table[n][m]) + (a * table[n + 1][m]);
	}
}
