PLUTO
Macros | Functions
math_table2D.c File Reference

Miscellaneous functions for handling 2D tables. More...

#include "pluto.h"
Include dependency graph for math_table2D.c:

Go to the source code of this file.

Macros

#define LOG_INTERPOLATION   NO
 

Functions

void PlotCubic (double a, double b, double c, double d)
 
void InitializeTable2D (Table2D *tab, double xmin, double xmax, int nx, double ymin, double ymax, int ny)
 
void FinalizeTable2D (Table2D *tab)
 
int LocateIndex (double *yarr, int beg, int end, double y)
 
int InverseLookupTable2D (Table2D *tab, double y, double f, double *x)
 
int Table2DInterpolate (Table2D *tab, double x, double y, double *f)
 
void WriteBinaryTable2D (char *fname, Table2D *tab)
 

Detailed Description

Miscellaneous functions for handling 2D tables.

Author
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
March 16, 2015

Definition in file math_table2D.c.

Macro Definition Documentation

#define LOG_INTERPOLATION   NO

Definition at line 194 of file math_table2D.c.

Function Documentation

void FinalizeTable2D ( Table2D *  tab)

Definition at line 121 of file math_table2D.c.

125 {
126  int i,j;
127 
128 /* -------------------------------------------------------
129  Compute forward differences
130  ------------------------------------------------------- */
131 
132  for (j = 0; j < tab->ny; j++){
133  for (i = 0; i < tab->nx-1; i++){
134  tab->dfx[j][i] = tab->f[j][i+1] - tab->f[j][i];
135  }}
136  for (j = 0; j < tab->ny-1; j++){
137  for (i = 0; i < tab->nx; i++){
138  tab->dfy[j][i] = tab->f[j+1][i] - tab->f[j][i];
139  }}
140 
141 }
TABLE2D::nx
int nx
Number of columns or points in the x direction.
Definition: structs.h:174
TABLE2D::f
double ** f
Definition: structs.h:186
TABLE2D::dfx
double ** dfx
Definition: structs.h:193
j
int j
Definition: analysis.c:2
TABLE2D::ny
int ny
Number of rows or points in the y direction.
Definition: structs.h:175
TABLE2D::dfy
double ** dfy
Definition: structs.h:194
i
int i
Definition: analysis.c:2

Here is the caller graph for this function:

void InitializeTable2D ( Table2D *  tab,
double  xmin,
double  xmax,
int  nx,
double  ymin,
double  ymax,
int  ny 
)

Allocate memory for the arrays contained in the *tab structure and generate a uniformly spaced grid in log10(x), log10(y) within the range provided by xmin, xmax, ymin, ymax. On output, the function initializes the following structure members:

  • nx, ny
  • lnxmin, lnxmax, lnymin, lnymax
  • lnx[], lny[]
  • dlnx, dlny
  • dlnx_1, dlny_1
  • x[], y[]
  • dx[], dy[]
Parameters
[in,out]tabpointer to a Table2D structure
[in]xminlower column limit.
[in]xmaxupper column limit.
[in]nxnumber of equally-spaced column bins (in log space)
[in]yminlower row limit.
[in]ymaxupper row limit.
[in]nynumber of equally-spaced row bins (in log space)

Definition at line 15 of file math_table2D.c.

40 {
41  int i,j;
42 
43  tab->nx = nx;
44  tab->ny = ny;
45 
46 /* ----------------------------------------------------------------------
47  If you do not initialize a structure variable, the effect depends on
48  whether it is has static storage (see Storage Class Specifiers) or
49  not.
50  If it is, members with integral types are initialized with 0 and
51  pointer members are initialized to NULL; otherwise, the value of
52  the structure's members is indeterminate.
53  ----------------------------------------------------------------------- */
54 
55  tab->f = ARRAY_2D(tab->ny, tab->nx, double); /* function value */
56  tab->a = ARRAY_2D(tab->ny, tab->nx, double); /* cubic coefficient (x^3) */
57  tab->b = ARRAY_2D(tab->ny, tab->nx, double);
58  tab->c = ARRAY_2D(tab->ny, tab->nx, double);
59  tab->d = ARRAY_2D(tab->ny, tab->nx, double);
60 
61  tab->defined = ARRAY_2D(tab->ny, tab->nx, char);
62  for (i = 0; i < tab->nx; i++){
63  for (j = 0; j < tab->ny; j++) tab->f[j][i] = 0.0;
64  for (j = 0; j < tab->ny; j++) tab->defined[j][i] = 1;
65  }
66 
67  tab->interpolation = LINEAR; /* Default */
68 
69  tab->lnx = ARRAY_1D(tab->nx, double);
70  tab->lny = ARRAY_1D(tab->ny, double);
71 
72  tab->x = ARRAY_1D(tab->nx, double);
73  tab->y = ARRAY_1D(tab->ny, double);
74 
75  tab->dx = ARRAY_1D(tab->nx, double);
76  tab->dy = ARRAY_1D(tab->ny, double);
77 
78  tab->dfx = ARRAY_2D(tab->ny, tab->nx, double);
79  tab->dfy = ARRAY_2D(tab->ny, tab->nx, double);
80 /*
81  -- (beta) function index bookeeping --
82  tab->fmin = ARRAY_1D(tab->ny, double);
83  tab->fmax = ARRAY_1D(tab->ny, double);
84  tab->df = ARRAY_1D(tab->ny, double);
85  tab->nf = 256;
86  tab->i = ARRAY_2D(tab->ny, tab->nf, int);
87 */
88 
89 /* ---------------------------------------------------------------
90  Compute table bounds
91  --------------------------------------------------------------- */
92 
93  tab->lnxmin = log10(xmin);
94  tab->lnxmax = log10(xmax);
95  tab->lnymin = log10(ymin);
96  tab->lnymax = log10(ymax);
97 
98  tab->dlnx = (tab->lnxmax - tab->lnxmin)/(double)(tab->nx - 1.0);
99  tab->dlny = (tab->lnymax - tab->lnymin)/(double)(tab->ny - 1.0);
100 
101  tab->dlnx_1 = 1.0/tab->dlnx;
102  tab->dlny_1 = 1.0/tab->dlny;
103 
104  for (i = 0; i < tab->nx; i++) {
105  tab->lnx[i] = tab->lnxmin + i*tab->dlnx;
106  tab->x[i] = pow(10.0, tab->lnx[i]);
107  }
108 
109  for (j = 0; j < tab->ny; j++) {
110  tab->lny[j] = tab->lnymin + j*tab->dlny;
111  tab->y[j] = pow(10.0, tab->lny[j]);
112  }
113 
114 /* -- table will not be uniform in x and y: compute spacings -- */
115 
116  for (i = 0; i < tab->nx-1; i++) tab->dx[i] = tab->x[i+1] - tab->x[i];
117  for (j = 0; j < tab->ny-1; j++) tab->dy[j] = tab->y[j+1] - tab->y[j];
118 
119 }
TABLE2D::y
double * y
array of y-values (not uniform)
Definition: structs.h:181
TABLE2D::nx
int nx
Number of columns or points in the x direction.
Definition: structs.h:174
TABLE2D::dlny_1
double dlny_1
Definition: structs.h:205
LINEAR
#define LINEAR
Definition: pluto.h:212
TABLE2D::f
double ** f
Definition: structs.h:186
TABLE2D::a
double ** a
Spline coefficient (x^3)
Definition: structs.h:188
TABLE2D::c
double ** c
Spline coefficient (x)
Definition: structs.h:190
TABLE2D::b
double ** b
Spline coefficient (x^2)
Definition: structs.h:189
TABLE2D::dy
double * dy
grid spacing array in the y direction (not uniform)
Definition: structs.h:183
TABLE2D::dfx
double ** dfx
Definition: structs.h:193
TABLE2D::dlnx_1
double dlnx_1
Definition: structs.h:204
j
int j
Definition: analysis.c:2
TABLE2D::lnxmax
double lnxmax
upper limit (in log10) in the x-direction
Definition: structs.h:199
TABLE2D::ny
int ny
Number of rows or points in the y direction.
Definition: structs.h:175
TABLE2D::dfy
double ** dfy
Definition: structs.h:194
TABLE2D::lny
double * lny
array of log10(y) values (uniform)
Definition: structs.h:185
TABLE2D::dlnx
double dlnx
uniform spacing in log10(x)
Definition: structs.h:202
TABLE2D::interpolation
int interpolation
LINEAR/SPLINE1.
Definition: structs.h:177
ARRAY_1D
#define ARRAY_1D(nx, type)
Definition: prototypes.h:170
TABLE2D::lnymax
double lnymax
upper limit (in log10) in the y-direction
Definition: structs.h:201
i
int i
Definition: analysis.c:2
TABLE2D::lnx
double * lnx
array of log10(x) values (uniform)
Definition: structs.h:184
TABLE2D::dlny
double dlny
uniform spacing in log10(y)
Definition: structs.h:203
TABLE2D::lnxmin
double lnxmin
lower limit (in log10) in the x-direction
Definition: structs.h:198
TABLE2D::lnymin
double lnymin
lower limit (in log10) in the y-direction
Definition: structs.h:200
ARRAY_2D
#define ARRAY_2D(nx, ny, type)
Definition: prototypes.h:171
TABLE2D::defined
char ** defined
Definition: structs.h:173
TABLE2D::d
double ** d
Spline coefficiten (1)
Definition: structs.h:191
TABLE2D::dx
double * dx
grid spacing array in the x direction (not uniform)
Definition: structs.h:182
TABLE2D::x
double * x
array of x-values (not uniform)
Definition: structs.h:180

Here is the caller graph for this function:

int InverseLookupTable2D ( Table2D *  tab,
double  y,
double  f,
double *  x 
)

Perform inverse lookup table interpolation: given a 2D table f(i,j)=f(x(i), y(j)) and the value of y, find the value of x. The algorithm proceeds by a combination of 1D lookup table algorithms:

  • Find the value of j such that y[j]<y<y[j+1];
  • The solution must then lie on the intermediate line between j and j+1 where the value of y is given.
  • To find the horizontal coordinate, we first locate initial and final indices ib and ie at j and j+1 using 1D lookup table.
  • Then construct a 1D array between these two indices values, perform again 1D lookup table to find the actual index i and invert the bilinear interpolant by solving

    \[ f = f_{i,j}(1 - x_n)(1 - y_n) + f_{i+1,j}x_n(1 - y_n) + f_{i,j+1}(1 - x_n)y_n + f_{i+1,j+1}x_ny_n; \]

    for $ x_n$ (normlized coordinate between 0 and 1).
Note
The table must be monotonic in both and y.
Parameters
[in]taba pointer to a Table2D structure
[in]ythe ordinata
[in]fthe value of the function
[out]*xthe value of the abscissa.
Returns
Returns 0 on success, otherwise:
  • -2 if y is below range;
  • 2 if y is above range.

Definition at line 196 of file math_table2D.c.

229 {
230  int i, j, k, i0, i1, ib, ie;
231  double lny, xn, yn;
232  double **ftab, **dfx, **dfy;
233  static double *f1;
234 
235  if (f1 == NULL) f1 = ARRAY_1D(8192, double);
236 
237  ftab = tab->f;
238  dfx = tab->dfx;
239  dfy = tab->dfy;
240 
241  lny = log10(y);
242  if (lny < tab->lnymin){
243  print ("! InverseLookupTable2D: lny outside range: %12.6e < %12.6e\n",
244  lny, tab->lnymin);
245  return -2;
246  }
247  if (lny > tab->lnymax){
248  print ("! InverseLookupTable2D: lny outside range: %12.6e > %12.6e\n",
249  lny, tab->lnymax);
250  return 2;
251  }
252 
253  j = INT_FLOOR((lny - tab->lnymin)*tab->dlny_1); /* Find row */
254  #if LOG_INTERPOLATION == YES
255  yn = (lny - tab->lny[j])/tab->dlny;
256  #else
257  yn = (y - tab->y[j])/tab->dy[j];
258  #endif
259 
260  i0 = LocateIndex(tab->f[j], 0, tab->nx-1, f);
261  i1 = LocateIndex(tab->f[j+1], 0, tab->nx-1, f);
262 
263  if (i0 < 0 || i1 < 0) return 2;
264 
265  ib = MIN(i0,i1);
266  ie = MAX(i0,i1) + 1;
267 
268  if (ie > (ib+1)){
269  for (i = ib; i <= ie; i++) {
270  f1[i] = tab->f[j][i]*(1.0 - yn) + tab->f[j+1][i]*yn;
271  }
272  i = LocateIndex(f1, ib, ie, f);
273  }else{
274  i = ib;
275  }
276 
277 /* -----------------------------------------------
278  Find xn in [0,1]
279  ----------------------------------------------- */
280 
281  if (tab->interpolation == LINEAR || tab->a[j][i] == 0.0){
282 
283  xn = f - ftab[j][i] - yn*dfy[j][i];
284  xn /= yn*(dfx[j+1][i] - dfx[j][i]) + dfx[j][i];
285 
286  }else if (tab->interpolation == SPLINE1 || tab->interpolation == SPLINE2){
287  double a,b,c,d;
288 
289 #if 1 /* Use Newton-Raphson */
290  int kmax = 6;
291  double fn, dfn, d2fn, dxn;
292 
293  a = tab->a[j][i]*(1.0 - yn) + tab->a[j+1][i]*yn;
294  b = tab->b[j][i]*(1.0 - yn) + tab->b[j+1][i]*yn;
295  c = tab->c[j][i]*(1.0 - yn) + tab->c[j+1][i]*yn;
296  d = tab->d[j][i]*(1.0 - yn) + tab->d[j+1][i]*yn - f;
297 
298  /* -- Initial guess is provided using linear interpolation -- */
299 
300  xn = f - ftab[j][i] - yn*dfy[j][i];
301  xn /= yn*(dfx[j+1][i] - dfx[j][i]) + dfx[j][i];
302 
303  /* -- solve cubic using Newton's method -- */
304 
305  for (k = 0; k <= kmax; k++){
306  fn = d + xn*(c + xn*(b + xn*a));
307  dfn = c + xn*(2.0*b + 3.0*a*xn);
308  dxn = fn/dfn;
309  xn -= dxn;
310  if (fabs(dxn) < 1.e-11) break;
311  if (k == kmax){
312  printf ("! InverseLookupTable2D(): too many iterations in Newton-cubic\n");
313  QUIT_PLUTO(1);
314  }
315  }
316 
317  /* -- solve cubic using Halley's method -- */
318 /*
319  for (k = 0; k <= kmax; k++){
320  fn = d + xn*(c + xn*(b + xn*a));
321  dfn = c + xn*(2.0*b + 3.0*a*xn);
322  d2fn = 2.0*b + 6.0*a*xn;
323  dxn = fn*dfn/(dfn*dfn - 0.5*d2fn*fn);
324  xn -= dxn;
325  if (fabs(dxn) < 1.e-11) break;
326  if (k == kmax){
327  printf ("! InverseLookupTable2D(): too many iterations in Newton-cubic\n");
328  QUIT_PLUTO(1);
329  }
330  }
331 */
332  if (xn < 0.0 || xn > 1.0) {
333  PlotCubic(a,b,c,d);
334  print ("! InverseLookupTable2D(): no root in [0,1]; i,j = %d,%d\n",i,j);
335  print ("! xn = %8.3e\n",xn);
336  print ("! Cubic tabulated in 'cubic.dat'\n");
337  QUIT_PLUTO(1);
338  }
339 
340 #else /* Use standard cubic solver (Numerical recipe, Sect. 5.6) */
341  double Q, R, A, B, R_Q;
342  double sn, cn, tn, sQ, th, xf;
343 
344  /* Re-write the cubic as xn^3 + aa*xn^2 + bb*xn + cc */
345 
346  a = tab->a[j][i]*(1.0 - yn) + tab->a[j+1][i]*yn;
347  b = tab->b[j][i]*(1.0 - yn) + tab->b[j+1][i]*yn;
348  c = tab->c[j][i]*(1.0 - yn) + tab->c[j+1][i]*yn;
349  d = tab->d[j][i]*(1.0 - yn) + tab->d[j+1][i]*yn - f;
350 
351  b /= a;
352  c /= a;
353  d /= a;
354 
355  Q = (b*b - 3.0*c)/9.0; /* Eq. [5.6.10] */
356  R = (2.0*b*b*b - 9.0*b*c + 27.0*d)/54.0; /* Eq. [5.6.10] */
357  R_Q = R/Q;
358  if (R_Q*R_Q <= Q){ /*** First case: 3 real roots ***/
359  sQ = sqrt(Q);
360  th = acos(R_Q/sQ); /* Eq. [5.6.11] */
361  tn = tan(th/3.0);
362  cn = 1.0/sqrt(1.0 + tn*tn);
363  sn = cn*tn;
364 
365  /* The three roots are given by Eq. (5.6.12) of Numerical Recipe
366  (we switch x3 <-> x2 so that one can verify that x1 < x2 < x3 always):
367 
368  x1 = -2.0*sQ*cos(th/3.0) - aa/3.0;
369  x2 = -2.0*sQ*cos((th - 2.0*CONST_PI)/3.0) - aa/3.0;
370  x3 = -2.0*sQ*cos((th + 2.0*CONST_PI)/3.0) - aa/3.0;
371 
372  To avoid loss of precsion we use
373 
374  cos((th + 2pi/3) = cos(th/3)*cos(2pi/3) - sin(th/3)*sin(2pi/3)
375  cos((th - 2pi/3) = cos(th/3)*cos(2pi/3) + sin(th/3)*sin(2pi/3)
376 
377  Using the fact that the roots must lie in [0,1] and that the spline
378  is monotonically increasing, this is how we pick up the solution: */
379 
380 
381  if ( a > 0.0 ){ /* Cubic goes from -inf to +inf: only x1 or x3 */
382  /* can be the correct root. */
383  xf = -b/3.0; /* Check inflection point to decide which one. */
384  if (xf < 0.0) xn = sQ*(cn + sqrt(3.0)*sn) - b/3.0; /* = x3 (Eq. [5.6.12])*/
385  else xn = -2.0*sQ*cn - b/3.0; /* = x1 (Eq. [5.6.12]) */
386  } else { /* Cubic goes from +inf to -inf: only x2 */
387  /* can be the correct root. */
388  xn = sQ*(cn - sqrt(3.0)*sn) - b/3.0; /* = x2 (Eq. [5.6.12]) */
389  }
390 
391  }else{ /*** Second case: one root only ***/
392  A = Q/R;
393  A = fabs(R)*(1.0 + sqrt(1.0 - A*A*Q));
394  A = -DSIGN(R)*pow(A, 1.0/3.0); /* Eq. [5.6.15] */
395  if (A == 0) B = 0.0;
396  else B = Q/A;
397  xn = (A + B) - b/3.0; /* Eq. [5.6.17] */
398  }
399 
400  if (xn < 0.0 || xn > 1.0) {
401  PlotCubic(a,b,c,d);
402  print ("! InverseLookupTable2D(): no root in [0,1]; i,j = %d,%d\n",i,j);
403  print ("! xn = %8.3e\n",xn);
404  print ("! R = %8.3e; Q = %8.3e\n",R,Q);
405  print ("! R^2 - Q^3 = %8.3e\n",R*R-Q*Q*Q);
406  print ("! |R/Q|/sqrt(Q) = %8.3e\n",fabs(R/Q)/sqrt(fabs(Q)));
407  print ("! a = %8.3e; b = %8.3e; c = %8.3e; d = %8.3e\n",a,b,c,d);
408  print ("! Cubic tabulated in 'cubic.dat'\n");
409  QUIT_PLUTO(1);
410  }
411 #endif
412 
413  } /* end if tab->interpolation == SPLINE1 */
414 
415 /* ----------------------------------------
416  Compute x = xn*dx + x[i]
417  ---------------------------------------- */
418 
419  #if LOG_INTERPOLATION == YES
420  (*x) = xn*tab->dlnx + tab->lnx[i];
421  (*x) = pow(10.0, *x);
422  #else
423  (*x) = xn*tab->dx[i] + tab->x[i];
424  #endif
425 
426  return 0;
427 }
a
static double a
Definition: init.c:135
MAX
#define MAX(a, b)
Definition: macros.h:101
TABLE2D::y
double * y
array of y-values (not uniform)
Definition: structs.h:181
SPLINE2
#define SPLINE2
Definition: pluto.h:220
LocateIndex
int LocateIndex(double *yarr, int beg, int end, double y)
Definition: math_table2D.c:144
SPLINE1
#define SPLINE1
Definition: pluto.h:219
TABLE2D::nx
int nx
Number of columns or points in the x direction.
Definition: structs.h:174
TABLE2D::dlny_1
double dlny_1
Definition: structs.h:205
LINEAR
#define LINEAR
Definition: pluto.h:212
DSIGN
#define DSIGN(x)
Definition: macros.h:110
TABLE2D::f
double ** f
Definition: structs.h:186
TABLE2D::a
double ** a
Spline coefficient (x^3)
Definition: structs.h:188
TABLE2D::c
double ** c
Spline coefficient (x)
Definition: structs.h:190
TABLE2D::b
double ** b
Spline coefficient (x^2)
Definition: structs.h:189
MIN
#define MIN(a, b)
Definition: macros.h:104
TABLE2D::dy
double * dy
grid spacing array in the y direction (not uniform)
Definition: structs.h:183
PlotCubic
void PlotCubic(double a, double b, double c, double d)
Definition: math_table2D.c:571
TABLE2D::dfx
double ** dfx
Definition: structs.h:193
INT_FLOOR
#define INT_FLOOR(z)
Definition: macros.h:98
j
int j
Definition: analysis.c:2
k
int k
Definition: analysis.c:2
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
TABLE2D::dfy
double ** dfy
Definition: structs.h:194
menu.c
tuple c
Definition: menu.py:375
TABLE2D::lny
double * lny
array of log10(y) values (uniform)
Definition: structs.h:185
TABLE2D::dlnx
double dlnx
uniform spacing in log10(x)
Definition: structs.h:202
TABLE2D::interpolation
int interpolation
LINEAR/SPLINE1.
Definition: structs.h:177
ARRAY_1D
#define ARRAY_1D(nx, type)
Definition: prototypes.h:170
sod.f1
tuple f1
Definition: sod.py:13
TABLE2D::lnymax
double lnymax
upper limit (in log10) in the y-direction
Definition: structs.h:201
i
int i
Definition: analysis.c:2
TABLE2D::lnx
double * lnx
array of log10(x) values (uniform)
Definition: structs.h:184
TABLE2D::dlny
double dlny
uniform spacing in log10(y)
Definition: structs.h:203
TABLE2D::lnymin
double lnymin
lower limit (in log10) in the y-direction
Definition: structs.h:200
TABLE2D::d
double ** d
Spline coefficiten (1)
Definition: structs.h:191
TABLE2D::dx
double * dx
grid spacing array in the x direction (not uniform)
Definition: structs.h:182
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
TABLE2D::x
double * x
array of x-values (not uniform)
Definition: structs.h:180

Here is the call graph for this function:

Here is the caller graph for this function:

int LocateIndex ( double *  yarr,
int  beg,
int  end,
double  y 
)

Given an array yarr[beg..end] and a given value y, returns the array index i such that yarr[i] < y < yarr[i+1] (for increasing arrays) or yarr[i+1] < y < yarr[i] (for decreasing arrays). yarr must be monotonically increasing or monotonically decreasing.

Parameters
[in]yarrthe 1D array
[in]beginitial index of the array
[in]endfinal index of the array
[in]ythe value to be searched for
Returns
On success return the index of the array. Otherwise returns -1.

Reference

  • "Numerical Recipes in C", Sect. 3.4 "How to Search an Ordered Table"

Definition at line 144 of file math_table2D.c.

165 {
166  int iu,im,il, ascnd;
167 
168  ascnd = (yarr[end] >= yarr[beg]); /* ascnd = 1 if table is increasing */
169 
170 /* --------------------------------------------
171  Check if y is bounded between max and min
172  -------------------------------------------- */
173 
174  if ( (ascnd && (y < yarr[beg] || y > yarr[end])) ||
175  (!ascnd && (y > yarr[beg] || y < yarr[end]))){
176 /* print ("! LocatIndex: element outside range ");
177  print ("(%12.6e not in [%12.6e, %12.6e]\n",y,yarr[beg],yarr[end]);*/
178  return -1;
179  }
180 
181  il = beg; /* Initialize lower */
182  iu = end; /* and upper limits. */
183 
184  while ( (iu - il) > 1) {
185  im = (iu + il) >> 1; /* Compute midpoint */
186  if ((y >= yarr[im]) == ascnd) il = im; /* Replace lower limit */
187  else iu = im; /* or upper limit */
188  }
189  if (y == yarr[beg]) return beg;
190  else if (y == yarr[end]) return end;
191  else return il;
192 }

Here is the caller graph for this function:

void PlotCubic ( double  a,
double  b,
double  c,
double  d 
)

Definition at line 571 of file math_table2D.c.

572 {
573  double t, f, dt = 1.e-1;
574  FILE *fp;
575 
576  fp = fopen("cubic.dat","w");
577  for (t = -50.0; t <= 50.0; t += dt){
578  f = d + t*(c + t*(b + t*a));
579  fprintf (fp, "%12.6e %12.6e\n",t,f);
580  }
581  fclose(fp);
582 }
a
static double a
Definition: init.c:135
menu.c
tuple c
Definition: menu.py:375
fp
FILE * fp
Definition: analysis.c:7

Here is the caller graph for this function:

int Table2DInterpolate ( Table2D *  tab,
double  x,
double  y,
double *  f 
)

Use bilinear interpolation to find the function f(x,y) from the 2D table *tab. Since the grid is equally spaced in log(x) and log(y), the (i,j) indices such that x[i] <= x < x[i+1] and y[j] <= y < y[j+1] are found by a simple division. Then bilinear interpolation can be done in either log or linear coordinates. The latter is expected to be slightly faster since it avoids one pow() operation.

Parameters
[in]*taba pointer to a Table2D structure
[in]xthe abscissa where interpolation is needed.
[in]ythe ordinata where interpolation is needed.
[out]*fthe interpolated value
Returns
- Return 0 on success, otherwise:
  • -1 if x is below column range;
  • 1 if x is above column range;
  • -2 if y is below row range;
  • 2 if y is above row range.

Definition at line 430 of file math_table2D.c.

454 {
455  int i,j;
456  double lnx, lny, xn, yn;
457  double s1, s2;
458 
459  lnx = log10(x); /* Take the log to find the indices */
460  lny = log10(y);
461 
462 /* -------------------------------------------------------
463  Check bounds
464  ------------------------------------------------------- */
465 
466  if (lnx < tab->lnxmin){
467  WARNING(
468  print ("! Table2DInterpolate: lnx outside range %12.6e < %12.6e\n",
469  lnx, tab->lnxmin);
470  )
471  return -1;
472  }
473 
474  if (lnx > tab->lnxmax){
475  WARNING(
476  print ("! Table2DInterpolate: lnx outside range: %12.6e > %12.6e\n",
477  lnx, tab->lnxmax);
478  )
479  return 1;
480  }
481 
482  if (lny < tab->lnymin){
483  WARNING(
484  print ("! Table2DInterpolate: lny outside range: %12.6e < %12.6e\n",
485  lny, tab->lnymin);
486  )
487  return -2;
488  }
489 
490  if (lny > tab->lnymax){
491  WARNING(
492  print ("! Table2DInterpolate: lny outside range: %12.6e > %12.6e\n",
493  lny, tab->lnymax);
494  )
495  return 2;
496  }
497 
498 /* ------------------------------------------------
499  Find column and row indices i and j
500  ------------------------------------------------ */
501 
502  i = INT_FLOOR((lnx - tab->lnxmin)*tab->dlnx_1);
503  j = INT_FLOOR((lny - tab->lnymin)*tab->dlny_1);
504 
505  #if LOG_INTERPOLATION == YES
506  xn = (lnx - tab->lnx[i])/tab->dlnx; /* Compute normalized log coordinates */
507  yn = (lny - tab->lny[j])/tab->dlny; /* on the unit square [0,1]x[0,1] */
508  #else
509  xn = (x - tab->x[i])/tab->dx[i]; /* Compute normalized linear coordinates */
510  yn = (y - tab->y[j])/tab->dy[j]; /* on the unit square [0,1]x[0,1] */
511  #endif
512 
513  if (tab->interpolation == LINEAR || tab->a[j][i] == 0.0){
514  (*f) = (tab->f[j][i] + (tab->f[j][i+1] - tab->f[j][i])*xn)*(1.0 - yn)
515  + (tab->f[j+1][i] + (tab->f[j+1][i+1] - tab->f[j+1][i])*xn)*yn;
516  }else if (tab->interpolation == SPLINE1 || tab->interpolation == SPLINE2){
517  s1 = tab->d[j][i] + xn*(tab->c[j][i] + xn*(tab->b[j][i] + xn*tab->a[j][i]));
518  s2 = tab->d[j+1][i] + xn*(tab->c[j+1][i] + xn*(tab->b[j+1][i] + xn*tab->a[j+1][i]));
519  *f = s1*(1.0 - yn) + s2*yn;
520  }else{
521  print1 ("! Table2DInterpolate(): table interpolation not corretly defined\n");
522  QUIT_PLUTO(1);
523  }
524  return 0; /* success */
525 }
TABLE2D::y
double * y
array of y-values (not uniform)
Definition: structs.h:181
SPLINE2
#define SPLINE2
Definition: pluto.h:220
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
SPLINE1
#define SPLINE1
Definition: pluto.h:219
WARNING
#define WARNING(a)
Definition: macros.h:217
LINEAR
#define LINEAR
Definition: pluto.h:212
TABLE2D::f
double ** f
Definition: structs.h:186
TABLE2D::a
double ** a
Spline coefficient (x^3)
Definition: structs.h:188
TABLE2D::c
double ** c
Spline coefficient (x)
Definition: structs.h:190
TABLE2D::b
double ** b
Spline coefficient (x^2)
Definition: structs.h:189
TABLE2D::dy
double * dy
grid spacing array in the y direction (not uniform)
Definition: structs.h:183
INT_FLOOR
#define INT_FLOOR(z)
Definition: macros.h:98
j
int j
Definition: analysis.c:2
if
if(divB==NULL)
Definition: analysis.c:10
TABLE2D::lnxmax
double lnxmax
upper limit (in log10) in the x-direction
Definition: structs.h:199
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
TABLE2D::lny
double * lny
array of log10(y) values (uniform)
Definition: structs.h:185
TABLE2D::dlnx
double dlnx
uniform spacing in log10(x)
Definition: structs.h:202
TABLE2D::interpolation
int interpolation
LINEAR/SPLINE1.
Definition: structs.h:177
TABLE2D::lnymax
double lnymax
upper limit (in log10) in the y-direction
Definition: structs.h:201
i
int i
Definition: analysis.c:2
TABLE2D::lnx
double * lnx
array of log10(x) values (uniform)
Definition: structs.h:184
TABLE2D::dlny
double dlny
uniform spacing in log10(y)
Definition: structs.h:203
TABLE2D::lnxmin
double lnxmin
lower limit (in log10) in the x-direction
Definition: structs.h:198
TABLE2D::lnymin
double lnymin
lower limit (in log10) in the y-direction
Definition: structs.h:200
TABLE2D::d
double ** d
Spline coefficiten (1)
Definition: structs.h:191
TABLE2D::dx
double * dx
grid spacing array in the x direction (not uniform)
Definition: structs.h:182
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
TABLE2D::x
double * x
array of x-values (not uniform)
Definition: structs.h:180

Here is the call graph for this function:

Here is the caller graph for this function:

void WriteBinaryTable2D ( char *  fname,
Table2D *  tab 
)

The binary table is a compact format used to write a 2D array together with simple structured coordinates. The file consists of the following information: 

 nx
 ny
 <x[0]..x[nx-1]>
 <y[0]..y[ny-1]>
 <q[0][0]..q[0][nx-1]
  q[1][0]..q[1][nx-1]
    .......
  q[ny-1][0]..q[ny-1][nx-1]>

All fields are written in binary format using double precision arithmetic.

Definition at line 528 of file math_table2D.c.

547 {
548  int i,j;
549  double scrh;
550  FILE *fp;
551 
552  fp = fopen(fname,"wb");
553 
554  scrh = (double)tab->nx;
555  fwrite(&scrh, sizeof(double), 1, fp);
556 
557  scrh = (double)tab->ny;
558  fwrite(&scrh, sizeof(double), 1, fp);
559 
560  fwrite(tab->lnx, sizeof(double), tab->nx, fp);
561  fwrite(tab->lny, sizeof(double), tab->ny, fp);
562 
563  for (j = 0; j < tab->ny; j++){
564  fwrite (tab->f[j], sizeof(double), tab->nx, fp);
565  }
566 
567  fprintf (fp,"\n");
568  fclose(fp);
569 }
TABLE2D::nx
int nx
Number of columns or points in the x direction.
Definition: structs.h:174
configure.scrh
tuple scrh
Definition: configure.py:200
TABLE2D::f
double ** f
Definition: structs.h:186
j
int j
Definition: analysis.c:2
TABLE2D::ny
int ny
Number of rows or points in the y direction.
Definition: structs.h:175
TABLE2D::lny
double * lny
array of log10(y) values (uniform)
Definition: structs.h:185
i
int i
Definition: analysis.c:2
TABLE2D::lnx
double * lnx
array of log10(x) values (uniform)
Definition: structs.h:184
fp
FILE * fp
Definition: analysis.c:7

Here is the caller graph for this function: