functabl.cpp

Go to the documentation of this file.

#include <../../nrnconf.h>
#include "hoc.h"

#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(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(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(double frac, double x1, double x2)
{
    if (frac > 0.) {
        return (1. - frac)*x1 + frac*x2;
    }else{
        return x1;
    }
}

double hoc_func_table(void* vpft, int n, double* args)
{
    int i, j;
    double* tab;
    FuncTable* ft = (FuncTable*)vpft;
    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
}

void hoc_spec_table(void** vppt, int n)
{
    int i, argcnt;
    FuncTable** ppt = (FuncTable**)vppt;
    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);
  }
}