CORSIKA add-on package IACT/ATMO:  Version 1.70 (August 2025)
rpolator.h
Go to the documentation of this file.
1 /* ============================================================================
2 
3  Copyright (C) 2017, 2018, 2019 Konrad Bernloehr
4 
5  This file is part of sim_telarray (also known as sim_hessarray)
6  and also part of the IACT/atmo package for CORSIKA.
7 
8  The software in this file is free software: you can redistribute it
9  and/or modify it under the terms of the GNU Lesser General Public
10  License as published by the Free Software Foundation; either
11  version 2.1 of the License, or (at your option) any later version.
12 
13  This program is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU Lesser General Public License for more details.
17 
18  You should have received a copy of the GNU Lesser General Public License
19  along with this program. If not, see <http://www.gnu.org/licenses/>.
20 
21 ============================================================================ */
22 
23 /* =============================================================== */
30 /* =============================================================== */
31 
32 #ifndef RPOLATOR_H__LOADED
33 #define RPOLATOR_H__LOADED 1
34 
35 #define WITH_RPOLATOR 1
36 
37 #ifdef __cplusplus
38 extern "C" {
39 #endif
40 
44 {
45  double a, b, c, d;
46 };
47 typedef struct cubic_params CsplinePar;
48 
51 struct rpol_table
52 {
53  int ndim;
54  size_t nx, ny;
55  double *x;
56  double *y;
57  double *z;
58  double *zxmax;
59  double *zxmin;
60  double aux;
61  double xmin, xmax;
62  double dx, dxi;
63  double xrise;
64  double ymin, ymax;
65  double dy, dyi;
66  double yrise;
67  double zmin, zmax;
68  char *fname;
69  char *options;
71  int scheme;
73  int clipping;
74  int remapped;
77  int zxreq;
78  int logs, xlog, ylog, zlog;
80  int use_count;
81 };
82 typedef struct rpol_table RpolTable;
83 
84 /* rpolator.c */
85 int read_table(const char *fname, int maxrow, double *col1, double *col2);
86 int read_table2(const char *fname, int maxrow, double *col1, double *col2);
87 int read_table3(const char *fname, int maxrow, double *col1, double *col2, double *col3);
88 int read_table4(const char *fname, int maxrow, double *col1, double *col2, double *col3, double *col4);
89 int read_table5(const char *fname, int maxrow, double *col1, double *col2, double *col3, double *col4, double *col5);
90 
91 int read_table_v(const char *fname, FILE *fptr, size_t *nrow, size_t ncol, double ***col, const size_t *selcol);
92 
93 struct rpol_table *read_rpol1d_table(const char *fn);
94 struct rpol_table *read_rpol2d_table(const char *fn, const char *ymarker);
95 struct rpol_table *read_rpol3d_table(const char *fn);
96 struct rpol_table *read_rpol_table(const char *fname, int nd, const char *ymarker, const char *options);
97 
98 struct rpol_table *simple_rpol1d_table(const char *label, double *x, double *y, int n, int clip);
99 
100 int rpol_is_verbose(void);
101 int rpol_set_verbose(int lvl);
102 
103 void rpol_info(struct rpol_table *rpt);
104 void rpol_info_lvl (struct rpol_table *rpt, int lvl);
105 void rpol_free(struct rpol_table *rpt, int removing);
106 void rpol_check_equi_range(struct rpol_table *rpt);
107 
108 double rpol(double *x, double *y, int n, double xp);
109 double rpol_nearest(double *x, double *y, int n, double xp, int eq, int clip);
110 double rpol_linear(double *x, double *y, int n, double xp, int eq, int clip);
111 double rpol_2nd_order(double *x, double *y, int n, double xp, int eq, int clip);
112 CsplinePar* set_1d_cubic_params (double *x, double *y, int n, int clamped);
113 double rpol_cspline ( double *x, double *y, const CsplinePar *csp, int n, double xp, int eq, int clip );
114 double rpol_2d_linear(double *x, double *y, double *z, int nx, int ny, double xp, double yp, int eq, int clip);
115 
116 double rpolate_1d(struct rpol_table *rpt, double x, int scheme);
117 double rpolate_1d_lin(struct rpol_table *rpt, double x);
118 double rpolate_2d(struct rpol_table *rpt, double x, double y, int scheme);
119 double rpolate(struct rpol_table *rpt, double x, double y, int scheme);
120 
121 #ifdef __cplusplus
122 }
123 #endif
124 
125 #endif
CsplinePar * set_1d_cubic_params(double *x, double *y, int n, int clamped)
Set up cubic spline parameters for n-1 intervals resulting from n data points.
Definition: rpolator.c:2740
int rpol_set_verbose(int lvl)
Set the verbosity level for the rpolator code, return old level.
Definition: rpolator.c:636
int rpol_is_verbose(void)
Report what verbosity level is set for the rpolator code.
Definition: rpolator.c:620
void rpol_info_lvl(struct rpol_table *rpt, int lvl)
Report table info at given temporary verbosity level.
Definition: rpolator.c:787
int read_table_v(const char *fname, FILE *fptr, size_t *nrow, size_t ncol, double ***col, const size_t *selcol)
Read tables any length (up to some ridiculous maximum) with the requested columns either in natural o...
Definition: rpolator.c:349
struct rpol_table * read_rpol1d_table(const char *fn)
Simplified version for loading a 1-d table with two columns.
Definition: rpolator.c:574
void rpol_info(struct rpol_table *rpt)
Show information about given interpolation table.
Definition: rpolator.c:649
int read_table3(const char *fname, int maxrow, double *col1, double *col2, double *col3)
read_table3() and so on have more columns than read_table but are still only suitable for 1-D interpo...
Definition: rpolator.c:162
void rpol_free(struct rpol_table *rpt, int removing)
Free a previously allocated interpolation table data structure.
Definition: rpolator.c:806
double rpolate_2d(struct rpol_table *rpt, double x, double y, int scheme)
High-level interpolation function (user code only has to keep a pointer to the allocated object) limi...
Definition: rpolator.c:3236
double rpolate_1d_lin(struct rpol_table *rpt, double x)
High-level interpolation function (user code only has to keep a pointer to the allocated object) limi...
Definition: rpolator.c:3213
double rpol_2d_linear(double *x, double *y, double *z, int nx, int ny, double xp, double yp, int eq, int clip)
Linear interpolation in 2-D.
Definition: rpolator.c:2958
double rpolate(struct rpol_table *rpt, double x, double y, int scheme)
High-level interpolation function (user code only has to keep a pointer to the allocated object) gene...
Definition: rpolator.c:3347
double rpol_linear(double *x, double *y, int n, double xp, int eq, int clip)
Linear interpolation in 1-D with with either direct access for equidistant table or with binary searc...
Definition: rpolator.c:2501
double rpol_2nd_order(double *x, double *y, int n, double xp, int eq, int clip)
Second/third order interpolation in 1-D with clipping option outside range.
Definition: rpolator.c:2581
double rpol_cspline(double *x, double *y, const CsplinePar *csp, int n, double xp, int eq, int clip)
Cubic spline interpolation in 1-D with clipping option outside range.
Definition: rpolator.c:2896
int read_table(const char *fname, int maxrow, double *col1, double *col2)
Low-level reading of 2-column data tables up to given number of data rows.
Definition: rpolator.c:100
struct rpol_table * read_rpol_table(const char *fname, int nd, const char *ymarker, const char *options)
General function for loading interpolation table combines 1-D and 2-D grid case.
Definition: rpolator.c:1141
struct rpol_table * read_rpol2d_table(const char *fn, const char *ymarker)
Simplified version for loading a 2-d table 1+ny columns (ny=available y values in header line).
Definition: rpolator.c:592
struct rpol_table * read_rpol3d_table(const char *fn)
Simplified version for loading a 2-d table with x/y/z in three columns (x/y intervals rectangular).
Definition: rpolator.c:613
struct rpol_table * simple_rpol1d_table(const char *label, double *x, double *y, int n, int clip)
A simplified way of setting up a 1-D rpol table for local use, without reading any files.
Definition: rpolator.c:2217
double rpolate_1d(struct rpol_table *rpt, double x, int scheme)
High-level interpolation function (user code only has to keep a pointer to the allocated object) limi...
Definition: rpolator.c:3083
double rpol(double *x, double *y, int n, double xp)
Linear interpolation with binary search algorithm.
Definition: rpolator.c:2388
double rpol_nearest(double *x, double *y, int n, double xp, int eq, int clip)
Nearest value (not actually interpolation) in 1-D with with either direct access for equidistant tabl...
Definition: rpolator.c:2425
Cubic spline interpolation (natural cubic splines = scheme 3, clampled cubic splines = scheme 4)
Definition: rpolator.h:44
double d
r=xp-x[i], yp = a + b*r + c*r^2 + d*r^3 = ((d*r + c) * r + b) * r + a;
Definition: rpolator.h:45
Structure describing an interpolation table, interpolation scheme and selected options.
Definition: rpolator.h:52
double dxi
Step size in x and inverse of it, for equidistant only.
Definition: rpolator.h:62
double * zxmin
Optional min.
Definition: rpolator.h:59
double xmax
Range covered in x.
Definition: rpolator.h:61
double * y
Supporting points in y (if ndim=2 and ny>1)
Definition: rpolator.h:56
int remapped
If user code remaps x, y, and/or z values w.r.t.
Definition: rpolator.h:74
int use_count
Indicates how often a table is in use, but not safe enough to make it a smart pointer.
Definition: rpolator.h:80
size_t ny
No.
Definition: rpolator.h:54
int ndim
1 or 2 dimension(s), for the independent variable(s) that is.
Definition: rpolator.h:53
char * options
Options used in option parameter or NULL pointer.
Definition: rpolator.h:69
double aux
This auxiliary value may be useful to the application but is not used internally.
Definition: rpolator.h:60
int clipping
0: Extrapolate with edge value, 1: zero outside range.
Definition: rpolator.h:73
double * z
Table values (size nx for 1-dim or nx*ny for 2-dim)
Definition: rpolator.h:57
double xrise
+1 if x values are in ascending order; -1 for descending
Definition: rpolator.h:63
char * fname
Name of the file from which the table was loaded (includes all options).
Definition: rpolator.h:68
int zlog
Log applied to any x/y axis, to x axis, y axis, z axis?
Definition: rpolator.h:78
int zxreq
Flag activated when options indicate that a set of zxmax values should be provided.
Definition: rpolator.h:77
double yrise
+1 if y values are in ascending order; -1 for descending
Definition: rpolator.h:66
double zmax
Range covered in result values.
Definition: rpolator.h:67
double dyi
Step size in y and inverse of it, for equidistant only (2-D).
Definition: rpolator.h:65
double * x
Supporting points in x.
Definition: rpolator.h:55
double * zxmax
Optional max.
Definition: rpolator.h:58
int equidistant
Equidistant support points make life much easier (bit 0: x, bit 1: y).
Definition: rpolator.h:70
double ymax
Range covered in y (if ndim=2)
Definition: rpolator.h:64
int scheme
Requested interpolation scheme (0: nearest, 1: linear, 2: quadratic, 3: natural cubic spline,...
Definition: rpolator.h:71
CsplinePar * csp
Cubic spline parameters (scheme 3 and 4 only), need one-time initialisation.
Definition: rpolator.h:79