interp.c File Reference

#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#include <math.h>
#include "interp.h"
#include "mymalloc.h"

Include dependency graph for interp.c:

Go to the source code of this file.

Functions
void	interp_init (Interp *obj, int Ndim, int64_t *dims)
void	interp_init_dim (Interp *obj, int d, double Min, double Max)
static ptrdiff_t	linearindex (ptrdiff_t *strides, int *xi, int Ndim)
double	interp_eval (Interp *obj, double *x, double *ydata, int *status)
double	interp_eval_periodic (Interp *obj, double *x, double *ydata)
void	interp_destroy (Interp *obj)

Function Documentation

interp_destroy()

void interp_destroy

(

Interp *

obj

)

Definition at line 172 of file interp.c.

                                  {
    myfree(obj->data);
}

References Interp::data, and myfree.

Referenced by test_interp().

Here is the caller graph for this function:

interp_eval()

double interp_eval	(	Interp *	obj,
		double *	x,
		double *	ydata,
		int *	status
	)

Definition at line 72 of file interp.c.

                                                                           {
    int d;
    if(status == NULL) {
        status = (int *) alloca(sizeof(int) * obj->Ndim);
    }
    int * xi = (int *) alloca(sizeof(int) * obj->Ndim);
    double * f = (double *) alloca(sizeof(double) * obj->Ndim);
 
    for(d = 0; d < obj->Ndim; d++) {
        double xd = (x[d] - obj->Min[d]) / obj->Step[d];
        if (x[d] < obj->Min[d]) {
            status[d] = -1;
            xi[d] = 0;
            f[d] = 0;
        } else
        if (x[d] > obj->Max[d]) {
            status[d] = 1;
            xi[d] = obj->dims[d] - 1;
            f[d] = 0;
        } else {
            xi[d] = floor(xd);
            f[d] = xd - xi[d];
            status[d] = 0;
        }
    }
 
    double ret = 0;
    /* the origin, "this point" */
    ptrdiff_t l0 = linearindex(obj->strides, xi, obj->Ndim);
 
    int i;
    /* for each point covered by the filter */
    for(i = 0; i < obj->fsize; i ++) {
        double filter = 1.0;
        ptrdiff_t l = l0;
        int skip = 0;
        for(d = 0; d < obj->Ndim; d++ ) {
            int foffset = (i & (1 << d))?1:0;
            if(f[d] == 0 && foffset == 1) {
                /* on this dimension the second data point 
                 * is not needed */
                skip = 1;
                break;
            }
 
            /* 
             * are we on this point or next point?
             *
             * weight on next point is f[d]
             * weight on this point is 1 - f[d] 
             * */
            filter *= foffset?f[d] : (1 - f[d]);
            l += foffset * obj->strides[d];
        }
        if(!skip) {
            ret += ydata[l] * filter;
        }
    }
    return ret;
}

References Interp::dims, Interp::fsize, linearindex(), Interp::Max, Interp::Min, Interp::Ndim, Interp::Step, and Interp::strides.

Referenced by TableMetalCoolingRate(), and test_interp().

Here is the call graph for this function:

Here is the caller graph for this function:

interp_eval_periodic()

double interp_eval_periodic	(	Interp *	obj,
		double *	x,
		double *	ydata
	)

Definition at line 134 of file interp.c.

                                                                      {
    int * xi = (int *) alloca(sizeof(int) * obj->Ndim);
    int * xi1 = (int *) alloca(sizeof(int) * obj->Ndim);
    double * f = (double *) alloca(sizeof(double) * obj->Ndim);
 
    int d;
    for(d = 0; d < obj->Ndim; d++) {
        double xd = (x[d] - obj->Min[d]) / obj->Step[d];
        xi[d] = floor(xd);
        f[d] = xd - xi[d];
    }
 
    double ret = 0;
    /* the origin, "this point" */
 
    int i;
    /* for each point covered by the filter */
    for(i = 0; i < obj->fsize; i ++) {
        double filter = 1.0;
        for(d = 0; d < obj->Ndim; d++ ) {
            int foffset = (i & (1 << d))?1:0;
            xi1[d] = xi[d] + foffset;
            while(xi1[d] >= obj->dims[d]) xi1[d] -= obj->dims[d];
            while(xi1[d] < 0 ) xi1[d] += obj->dims[d];
            /* 
             * are we on this point or next point?
             *
             * weight on next point is f[d]
             * weight on this point is 1 - f[d] 
             * */
            filter *= foffset?f[d] : (1 - f[d]);
        }
        ptrdiff_t l = linearindex(obj->strides, xi1, obj->Ndim);
        ret += ydata[l] * filter;
    }
    return ret;
}

References Interp::dims, Interp::fsize, linearindex(), Interp::Min, Interp::Ndim, Interp::Step, and Interp::strides.

Referenced by test_interp().

Here is the call graph for this function:

Here is the caller graph for this function:

interp_init()

void interp_init	(	Interp *	obj,
		int	Ndim,
		int64_t *	dims
	)

Definition at line 9 of file interp.c.

                                                         {
    ptrdiff_t N = 1;
    int d;
    for(d = 0 ; d < Ndim; d ++) {
        N *= dims[d];     
    }
 
    /* alloc memory */
    obj->Ndim = Ndim;
    obj->data = mymalloc("interp_data", 0
        +   sizeof(double) * Ndim * 3
        +   sizeof(ptrdiff_t) * Ndim
        +   sizeof(int) * Ndim);
 
    obj->dims = (int *) obj->data;
    obj->strides = (ptrdiff_t*) (obj->dims + Ndim);
    obj->Min = (double*) (obj->strides + Ndim);
    obj->Step = obj->Min + Ndim;
    obj->Max = obj->Step + Ndim;
 
    /* fillin strides */
    N = 1;
    for(d = Ndim - 1 ; d >= 0; d --) {
        obj->dims[d] = dims[d];
        /* column first, C ordering */
        obj->strides[d] = N;
        N *= dims[d];     
    }
 
    int fsize = 1;
    for(d = 0; d < obj->Ndim; d++) {
        fsize *= 2;
    }
    obj->fsize = fsize;
}

References Interp::data, Interp::dims, Interp::fsize, Interp::Max, Interp::Min, mymalloc, Interp::Ndim, Interp::Step, and Interp::strides.

Referenced by get_local_UVBG(), InitMetalCooling(), and test_interp().

Here is the caller graph for this function:

interp_init_dim()

void interp_init_dim	(	Interp *	obj,
		int	d,
		double	Min,
		double	Max
	)

Definition at line 50 of file interp.c.

                                                                  {
    obj->Min[d] = Min;
    obj->Max[d] = Max;
    obj->Step[d] = (Max - Min) / (obj->dims[d] - 1);
}

References Interp::dims, Interp::Max, Interp::Min, and Interp::Step.

Referenced by get_local_UVBG(), InitMetalCooling(), and test_interp().

Here is the caller graph for this function:

linearindex()

static ptrdiff_t linearindex ( ptrdiff_t *  strides, int *  xi, int  Ndim  )

static

Definition at line 57 of file interp.c.

                                                                      {
    int d;
    ptrdiff_t rt = 0;
    for(d = 0; d < Ndim; d++) {
        rt += strides[d] * xi[d] ;
    }
    return rt;
}

Referenced by interp_eval(), and interp_eval_periodic().

Here is the caller graph for this function: