#include <../../nrnconf.h>
#include "hoc.h"
extern int hoc_is_double_arg();
extern double* getarg();
extern double* hoc_pgetarg();

#define INTERPOLATE 1

typedef struct TableArg {
	int nsize;
	double* argvec;	/* if nil use min,max */
	double min;
	double max;
#if INTERPOLATE
	double frac;	/* temporary storage */
#endif
} TableArg;

typedef struct FuncTable {
	double* table;
	TableArg* targs;	
	double value; /* if constant this is it */
} FuncTable;

static int arg_index(ta, x)
	TableArg* ta;
	double x;
{
	int j;
		
	if (ta->argvec) {
		int t0, t1;
#if INTERPOLATE
		ta->frac = 0.;
#endif
		t0 = 0;
		t1 = ta->nsize-1;
		if (x <= ta->argvec[t0]) {
			j = 0;
		}else if (x >= ta->argvec[t1]) {
			j = t1;
		}else{
			while (t0 < (t1-1)) {
				j = (t0 + t1)/2;
#if 0
printf("x[%d]=%g  x[%d]=%g  x[%d]=%g\n", t0, ta->argvec[t0], j, ta->argvec[j], t1, ta->argvec[t1]);
#endif
				if (ta->argvec[j] <= x) {
					t0 = j;
				}else{
					t1 = j;
				}
			}
			j = t0;
#if INTERPOLATE
ta->frac = (x - ta->argvec[j])/(ta->argvec[j+1] - ta->argvec[j]);
#if 0
printf("x[%d]=%g  frac=%g  x=%g\n", j, ta->argvec[j], ta->frac, x);
#endif
#endif
		}
	}else{
#if INTERPOLATE
		ta->frac = 0.;
#endif
		if (x <= ta->min) {
			j = 0;
		}else if (x >= ta->max) {
			j = ta->nsize -1;
		}else{
			double d, x1;
			d = (ta->max - ta->min)/((double)(ta->nsize - 1));
			x1 = (x - ta->min)/d;
			j = (int)x1;
#if INTERPOLATE
			ta->frac = x1 - (double)j;
#endif
		}
	}
	return j;
}

static double inter(frac, tab, j)
	double frac;
	double* tab;
	int j;
{
	if (frac > 0.) {
		return (1. - frac)*tab[j] + frac*tab[j+1];
	}else{
		return tab[j];
	}
}

static double interp(frac, x1, x2)
	double frac, x1, x2;
{
	if (frac > 0.) {
		return (1. - frac)*x1 + frac*x2;
	}else{
		return x1;
	}
}

double hoc_func_table(ft, n, args)
	FuncTable* ft;
	int n;
	double* args;
{
	int i, j;
	double* tab;
	if (!ft) {
		hoc_execerror("table not specified in hoc_func_table", (char*)0);
	}
	tab = ft->table;
	j = 0;
	for (i=0; i < n; ++i) {
		j = (j * ft->targs[i].nsize) + arg_index(ft->targs + i, args[i]);
/*
printf("calculating j: i=%d args=%g j=%d frac=%g\n", i, args[i], j, ft->targs[i].frac);
*/
	}
#if INTERPOLATE
	if (n == 1) {
		return inter(ft->targs[0].frac, tab, j);
	}else if (n == 2) {
		/* adjacent indices ar adjacent values of y dimension */
		double y1, y2;
		double fy = ft->targs[1].frac;
		y1 = inter(fy, tab, j);
/*
printf("calculating y1: fy=%g j=%d y1=%g, tabj=%g\n", fy, j, y1, tab[j]);
*/
		if (ft->targs[0].frac) {	/* x dimension fraction */
			int j1 = j + ft->targs[1].nsize;
			y2 = inter(fy, tab, j1);
/*
printf("calculating y2: fx=%g fy=%g j1=%d y2=%g tabj1=%g\n", ft->targs[0].frac, fy, j1, y2, tab[j1]);
*/
			return interp(ft->targs[0].frac, y1, y2);
		}else{
			return y1;
		}
	}else{
		return tab[j];
	}
#else
	return tab[j];
#endif
}

hoc_spec_table(ppt, n)
	FuncTable** ppt;
	int n;
{
	int i, argcnt;
	FuncTable* ft;
	TableArg* ta;
	if (!*ppt) {
		*ppt = (FuncTable*) ecalloc(1, sizeof(FuncTable));
		(*ppt)->targs = (TableArg*) ecalloc(n, sizeof(TableArg));
	}
	ft = *ppt;
	ta = ft->targs;
	argcnt = 2;
	if (!ifarg(2)) { /* set to constant */
		ft->value = *getarg(1);
		ft->table = &ft->value;
		for (i=0; i < n; ++i) {
			ta[i].nsize = 1;
			ta[i].argvec = (double*)0;
			ta[i].min = 1e20;
			ta[i].max = 1e20;
		}
		return;
	}
  if (hoc_is_object_arg(1)) {
	int ns;
	if (n > 1) {
hoc_execerror("Vector arguments allowed only for functions", "of one variable");
	}
	ns = vector_arg_px(1, &ft->table);
	ta[0].nsize = vector_arg_px(2, &ta[0].argvec);
	if (ns != ta[0].nsize) {
		hoc_execerror("Vector arguments not same size", (char*)0);
	}
  }else{
	for (i = 0; i < n; ++i) {
  		ta[i].nsize = *getarg(argcnt++);
		if (ta[i].nsize < 1) {
			hoc_execerror("size arg < 1 in hoc_spec_table", (char*)0);
		}
		if (hoc_is_double_arg(argcnt)) {
			ta[i].min = *getarg(argcnt++);
			ta[i].max = *getarg(argcnt++);
			if (ta[i].max < ta[i].min) {
				hoc_execerror("min > max in hoc_spec_table", (char*)0);
			}
			ta[i].argvec = (double*)0;
		}else{
			ta[i].argvec = hoc_pgetarg(argcnt++);
		}
	}
	ft->table = hoc_pgetarg(1);
  }
}