Improved interpolation

This commit is contained in:
Jim Monte 2020-04-25 19:27:07 +02:00 committed by Holger Vogt
parent 00fa875b9e
commit 554e594f92
4 changed files with 36 additions and 111 deletions

View File

@ -41,9 +41,9 @@ sf_alloc(int n /* number of elements */,
/* Use calloc so that any internal allocations will be set to NULL to /* Use calloc so that any internal allocations will be set to NULL to
* facilitate error recovery */ * facilitate error recovery */
if ((ptr = calloc(n, size)) == NULL) { if ((ptr = calloc((size_t) n, size)) == NULL) {
cm_message_printf("%s: cannot allocate %zd bytes : ", cm_message_printf("%s: cannot allocate %zd bytes : ",
__FILE__, n * size); __FILE__, (size_t) n * size);
return NULL; return NULL;
} }

View File

@ -4,6 +4,8 @@
#include <string.h> #include <string.h>
#include <math.h> #include <math.h>
#include "interp.h"
/*********************/ /*********************/
/* 3d geometry types */ /* 3d geometry types */
/*********************/ /*********************/

View File

@ -73,6 +73,8 @@ NON-STANDARD FEATURES
#include <sys/stat.h> #include <sys/stat.h>
#include "mada/eno2.h" #include "mada/eno2.h"
#include "support/gettokens.h"
#include "support/interp.h"
typedef struct { typedef struct {
int ix; /* size of array in x */ int ix; /* size of array in x */
@ -88,7 +90,6 @@ typedef struct {
double **table; /* f(xi, yj) */ double **table; /* f(xi, yj) */
} Table2_Data_t; } Table2_Data_t;
typedef Table2_Data_t Local_Data_t;
/*=== MACROS ===========================*/ /*=== MACROS ===========================*/
@ -98,33 +99,17 @@ typedef Table2_Data_t Local_Data_t;
#define DIR_PATHSEP "/" #define DIR_PATHSEP "/"
#endif #endif
#if defined(_MSC_VER)
#define strdup _strdup
#endif
/*=== LOCAL VARIABLES & TYPEDEFS =======*/ /*=== LOCAL VARIABLES & TYPEDEFS =======*/
struct filesource_state {
FILE *fp;
long pos;
unsigned char atend;
};
/*=== FUNCTION PROTOTYPE DEFINITIONS ===*/ /*=== FUNCTION PROTOTYPE DEFINITIONS ===*/
double BilinearInterpolation(double x, double y, static void cm_table2D_callback(ARGS, Mif_Callback_Reason_t reason);
int xind, int yind, double **td);
extern char *CNVgettok(char **s);
int cnv_get_spice_value(char *str, double *p_value);
extern int findCrossOver(double arr[], int n, double x);
static void free_local_data(Table2_Data_t *loc); static void free_local_data(Table2_Data_t *loc);
static inline double get_local_diff(int n, double *col, int ind); static inline double get_local_diff(int n, double *col, int ind);
static Table2_Data_t *init_local_data(const char *filename, int order); static Table2_Data_t *init_local_data(const char *filename, int interporder);
@ -174,10 +159,10 @@ static void cm_table2D_callback(ARGS,
{ {
switch (reason) { switch (reason) {
case MIF_CB_DESTROY: { case MIF_CB_DESTROY: {
Table2_Data_t *loc = STATIC_VAR(locdata); Table2_Data_t *loc = (Table2_Data_t *) STATIC_VAR(locdata);
if (loc) { if (loc) {
free_local_data(loc); free_local_data(loc);
STATIC_VAR(locdata) = loc = NULL; loc = (Table2_Data_t *) (STATIC_VAR(locdata) = NULL);
} }
break; break;
} /* end of case MIF_CB_DESTROY */ } /* end of case MIF_CB_DESTROY */
@ -262,7 +247,8 @@ void cm_table2D(ARGS) /* structure holding parms, inputs, outputs, etc. */
} }
/* return immediately if there was an initialization error */ /* return immediately if there was an initialization error */
if ((loc = STATIC_VAR(locdata)) == (Table2_Data_t *) NULL) { if ((loc = (Table2_Data_t *) STATIC_VAR(locdata)) ==
(Table2_Data_t *) NULL) {
return; return;
} }
@ -345,11 +331,10 @@ void cm_table2D(ARGS) /* structure holding parms, inputs, outputs, etc. */
} }
else { else {
Mif_Complex_t ac_gain; Mif_Complex_t ac_gain;
ac_gain.real = PARAM(gain) * derivval[0] / xdiff;
ac_gain.imag= 0.0; ac_gain.imag= 0.0;
ac_gain.real = PARAM(gain) * derivval[0] / xdiff;
AC_GAIN(out, inx) = ac_gain; AC_GAIN(out, inx) = ac_gain;
ac_gain.real = PARAM(gain) * derivval[1] / ydiff; ac_gain.real = PARAM(gain) * derivval[1] / ydiff;
ac_gain.imag= 0.0;
AC_GAIN(out, iny) = ac_gain; AC_GAIN(out, iny) = ac_gain;
} }
} /* end of function cm_table2D */ } /* end of function cm_table2D */
@ -357,7 +342,7 @@ void cm_table2D(ARGS) /* structure holding parms, inputs, outputs, etc. */
/* This function initializes local data */ /* This function initializes local data */
static Table2_Data_t *init_local_data(const char *filename, int order) static Table2_Data_t *init_local_data(const char *filename, int interporder)
{ {
int xrc = 0; int xrc = 0;
int ix = 0, /* elements in a row */ int ix = 0, /* elements in a row */
@ -373,7 +358,6 @@ static Table2_Data_t *init_local_data(const char *filename, int order)
size_t lFileRead; /* Length of file read in */ size_t lFileRead; /* Length of file read in */
int lLineCount; /* Current line number */ int lLineCount; /* Current line number */
size_t lTotalChar; /* Total characters read */ size_t lTotalChar; /* Total characters read */
int interporder; /* order of interpolation for eno */
Table2_Data_t *loc = (Table2_Data_t *) NULL; /* local data */ Table2_Data_t *loc = (Table2_Data_t *) NULL; /* local data */
@ -557,7 +541,6 @@ static Table2_Data_t *init_local_data(const char *filename, int order)
} }
/* generate table core */ /* generate table core */
interporder = order;
/* boundary limits set to param 'order' aren't recognized, /* boundary limits set to param 'order' aren't recognized,
so limit them here */ so limit them here */
if (interporder < 2) { if (interporder < 2) {
@ -570,21 +553,18 @@ static Table2_Data_t *init_local_data(const char *filename, int order)
int n1, int n2 : data dimensions */ int n1, int n2 : data dimensions */
if ((loc->newtable = sf_eno2_init( if ((loc->newtable = sf_eno2_init(
interporder, ix, iy)) == (sf_eno2) NULL) { interporder, ix, iy)) == (sf_eno2) NULL) {
cm_message_printf("eno2 initialization failure."); cm_message_send("eno2 initialization failure.");
xrc = -1; xrc = -1;
goto EXITPOINT; goto EXITPOINT;
} }
/* create table_data in memory */ /* create table_data in memory */
/* data [n2][n1] */ /* data [n2][n1] */
if ((loc->table = table_data = (double **) calloc((size_t) iy, if ((loc->table = table_data = (double **) calloc((size_t) iy,
sizeof(double *))) == (double **) NULL) { sizeof(double *))) == (double **) NULL) {
cm_message_printf("Unable to allocate data table."); cm_message_send("Unable to allocate data table.");
free(cFile); xrc = -1;
free(cThisLine); goto EXITPOINT;
free_local_data(loc);
return (Local_Data_t *) NULL;
} }
{ {
@ -594,16 +574,12 @@ static Table2_Data_t *init_local_data(const char *filename, int order)
sizeof(double))) == (double *) NULL) { sizeof(double))) == (double *) NULL) {
cm_message_printf("Unable to allocate data table " cm_message_printf("Unable to allocate data table "
"row %d", i + 1); "row %d", i + 1);
free(cFile); xrc = -1;
free(cThisLine); goto EXITPOINT;
free_local_data(loc);
return (Local_Data_t *) NULL;
} }
} }
} }
loc->table = table_data; /* give to local data structure */
/* continue reading f(x,y) values from cFile */ /* continue reading f(x,y) values from cFile */
lLineCount = 0; lLineCount = 0;
@ -635,10 +611,8 @@ static Table2_Data_t *init_local_data(const char *filename, int order)
if (lTotalChar >= lFileLen) { if (lTotalChar >= lFileLen) {
cm_message_printf("Not enough data in file %s", cm_message_printf("Not enough data in file %s",
filename); filename);
free(cFile); xrc = -1;
free(cThisLine); goto EXITPOINT;
free_local_data(loc);
return (Local_Data_t *) NULL;
} }
lLineCount--; /* we count only real lines */ lLineCount--; /* we count only real lines */
continue; continue;
@ -698,7 +672,7 @@ EXITPOINT:
/* Free memory allocations in Local_Data_t structure */ /* Free memory allocations in Table2_Data_t structure */
static void free_local_data(Table2_Data_t *loc) static void free_local_data(Table2_Data_t *loc)
{ {
if (loc == (Table2_Data_t *) NULL) { if (loc == (Table2_Data_t *) NULL) {

View File

@ -73,6 +73,7 @@ NON-STANDARD FEATURES
#include <sys/stat.h> #include <sys/stat.h>
#include "support/gettokens.h" #include "support/gettokens.h"
#include "support/interp.h"
#include "mada/eno2.h" #include "mada/eno2.h"
#include "mada/eno3.h" #include "mada/eno3.h"
@ -92,7 +93,6 @@ typedef struct {
double ***table; /* f(xi, yj, zk) */ double ***table; /* f(xi, yj, zk) */
} Table3_Data_t; } Table3_Data_t;
typedef Table3_Data_t Local_Data_t;
/*=== MACROS ===========================*/ /*=== MACROS ===========================*/
@ -102,10 +102,6 @@ typedef Table3_Data_t Local_Data_t;
#define DIR_PATHSEP "/" #define DIR_PATHSEP "/"
#endif #endif
#if defined(_MSC_VER)
#define strdup _strdup
#endif
/*=== LOCAL VARIABLES & TYPEDEFS =======*/ /*=== LOCAL VARIABLES & TYPEDEFS =======*/
@ -113,68 +109,21 @@ typedef Table3_Data_t Local_Data_t;
/*=== FUNCTION PROTOTYPE DEFINITIONS ===*/ /*=== FUNCTION PROTOTYPE DEFINITIONS ===*/
extern char *CNVgettok(char **s); static void cm_table3D_callback(ARGS, Mif_Callback_Reason_t reason);
extern double TrilinearInterpolation(double x, double y, double z, int xind, int yind, int zind, double ***td);
int cnv_get_spice_value(char *str, double *p_value);
extern int findCrossOver(double arr[], int n, double x);
static void free_local_data(Table3_Data_t *loc); static void free_local_data(Table3_Data_t *loc);
static inline double get_local_diff(int n, double *col, int ind); static inline double get_local_diff(int n, double *col, int ind);
static Table3_Data_t *init_local_data(const char *filename, int order); static Table3_Data_t *init_local_data(const char *filename, int order);
/*==============================================================================
FUNCTION cnv_get_spice_value()
AUTHORS
??? Bill Kuhn
MODIFICATIONS
30 Sep 1991 Jeffrey P. Murray
SUMMARY
This function takes as input a string token from a SPICE
deck and returns a floating point equivalent value.
INTERFACES
FILE ROUTINE CALLED
N/A N/A
RETURNED VALUE
Returns the floating point value in pointer *p_value. Also
returns an integer representing successful completion.
GLOBAL VARIABLES
NONE
NON-STANDARD FEATURES
NONE
==============================================================================*/
static void cm_table3D_callback(ARGS, static void cm_table3D_callback(ARGS,
Mif_Callback_Reason_t reason) Mif_Callback_Reason_t reason)
{ {
switch (reason) { switch (reason) {
case MIF_CB_DESTROY: { case MIF_CB_DESTROY: {
Table3_Data_t *loc = STATIC_VAR(locdata); Table3_Data_t *loc = (Table3_Data_t *) STATIC_VAR(locdata);
if (loc) { if (loc) {
free_local_data(loc); free_local_data(loc);
STATIC_VAR(locdata) = loc = NULL; loc = (Table3_Data_t *) (STATIC_VAR(locdata) = NULL);
} }
break; break;
} /* end of case MIF_CB_DESTROY */ } /* end of case MIF_CB_DESTROY */
@ -261,7 +210,8 @@ void cm_table3D(ARGS) /* structure holding parms, inputs, outputs, etc. */
} }
/* return immediately if there was an initialization error */ /* return immediately if there was an initialization error */
if ((loc = STATIC_VAR(locdata)) == (Table3_Data_t *) NULL) { if ((loc = (Table3_Data_t *) STATIC_VAR(locdata)) ==
(Table3_Data_t *) NULL) {
return; return;
} }
@ -362,15 +312,13 @@ void cm_table3D(ARGS) /* structure holding parms, inputs, outputs, etc. */
} }
else { else {
Mif_Complex_t ac_gain; Mif_Complex_t ac_gain;
ac_gain.real = PARAM(gain) * derivval[0] / xdiff;
ac_gain.imag= 0.0; ac_gain.imag= 0.0;
ac_gain.real = PARAM(gain) * derivval[0] / xdiff;
AC_GAIN(out, inx) = ac_gain; AC_GAIN(out, inx) = ac_gain;
ac_gain.real = PARAM(gain) * derivval[1] / ydiff; ac_gain.real = PARAM(gain) * derivval[1] / ydiff;
ac_gain.imag= 0.0;
AC_GAIN(out, iny) = ac_gain; AC_GAIN(out, iny) = ac_gain;
ac_gain.real = PARAM(gain) * derivval[2] / zdiff; ac_gain.real = PARAM(gain) * derivval[2] / zdiff;
ac_gain.imag= 0.0; AC_GAIN(out, inz) = ac_gain;
AC_GAIN(out, iny) = ac_gain;
} }
} /* end of function cm_table3D */ } /* end of function cm_table3D */
@ -619,7 +567,7 @@ static Table3_Data_t *init_local_data(const char *filename, int interporder)
int n1, int n2, int n3 : data dimensions */ int n1, int n2, int n3 : data dimensions */
if ((loc->newtable = sf_eno3_init( if ((loc->newtable = sf_eno3_init(
interporder, ix, iy, iz)) == (sf_eno3) NULL) { interporder, ix, iy, iz)) == (sf_eno3) NULL) {
cm_message_printf("eno3 initialization failure."); cm_message_send("eno3 initialization failure.");
xrc = -1; xrc = -1;
goto EXITPOINT; goto EXITPOINT;
} }
@ -628,7 +576,7 @@ static Table3_Data_t *init_local_data(const char *filename, int interporder)
/* data [n3][n2][n1] */ /* data [n3][n2][n1] */
if ((loc->table = table_data = (double ***) calloc((size_t) iz, if ((loc->table = table_data = (double ***) calloc((size_t) iz,
sizeof(double **))) == (double ***) NULL) { sizeof(double **))) == (double ***) NULL) {
cm_message_printf("Unable to allocate data table."); cm_message_send("Unable to allocate data table.");
xrc = -1; xrc = -1;
goto EXITPOINT; goto EXITPOINT;
} }
@ -752,7 +700,7 @@ EXITPOINT:
/* Free memory allocations in Local_Data_t structure */ /* Free memory allocations in Table3_Data_t structure */
static void free_local_data(Table3_Data_t *loc) static void free_local_data(Table3_Data_t *loc)
{ {
if (loc == (Table3_Data_t *) NULL) { if (loc == (Table3_Data_t *) NULL) {
@ -761,10 +709,11 @@ static void free_local_data(Table3_Data_t *loc)
/* Free data table and related values */ /* Free data table and related values */
if (loc->table) { if (loc->table) {
int i, j; int i;
int n_y = loc->iy; int n_y = loc->iy;
int n_z = loc->iz; int n_z = loc->iz;
for (i = 0; i < n_z; i++) { for (i = 0; i < n_z; i++) {
int j;
for (j = 0; j < n_y; j++) { for (j = 0; j < n_y; j++) {
free(loc->table[i][j]); free(loc->table[i][j]);
} }