23 #ifndef O2SCL_INTERP2_PLANAR_H 24 #define O2SCL_INTERP2_PLANAR_H 34 #include <gsl/gsl_combination.h> 36 #include <o2scl/err_hnd.h> 37 #include <o2scl/vector.h> 39 #ifndef DOXYGEN_NO_O2NS 136 void set_data(
size_t n_points, vec_t &x, vec_t &y, vec_t &f) {
138 O2SCL_ERR2(
"Must provide at least three points in ",
155 double minx=(*ux)[0], maxx=(*ux)[0];
156 for(
size_t i=1;i<
np;i++) {
157 if ((*
ux)[i]<minx) minx=(*ux)[i];
158 if ((*
ux)[i]>maxx) maxx=(*ux)[i];
165 double miny=(*uy)[0], maxy=(*uy)[0];
166 for(
size_t i=1;i<
np;i++) {
167 if ((*
uy)[i]<miny) miny=(*uy)[i];
168 if ((*
uy)[i]>maxy) maxy=(*uy)[i];
175 if (dx<=0.0 || dy<=0.0) {
184 double eval(
double x,
double y)
const {
185 double x1,
x2,
x3, y1, y2, y3, f;
187 eval_points(x,y,f,i1,x1,y1,i2,x2,y2,i3,x3,y3);
201 return eval(v[0],v[1]);
211 size_t &i1,
double &
x1,
double &y1,
212 size_t &i2,
double &
x2,
double &y2,
213 size_t &i3,
double &
x3,
double &y3)
const {
216 O2SCL_ERR(
"Data not set in interp2_planar::eval_points().",
225 double c1=sqrt(pow((x-(*
ux)[0])/dx,2.0)+pow((y-(*
uy)[0])/dy,2.0));
226 double c2=sqrt(pow((x-(*
ux)[1])/dx,2.0)+pow((y-(*
uy)[1])/dy,2.0));
227 double c3=sqrt(pow((x-(*
ux)[2])/dx,2.0)+pow((y-(*
uy)[2])/dy,2.0));
255 for(
size_t j=3;j<
np;j++) {
257 double c4=sqrt(pow((x-(*
ux)[i4])/dx,2.0)+pow((y-(*
uy)[i4])/dy,2.0));
271 double a, b, c, den, z1, z2, z3;
278 den=(-x2*y1+x3*y1+x1*y2-x3*y2-x1*y3+x2*y3);
289 ubvector_size_t order(np);
290 for(
size_t i=0;i<
np;i++) {
291 dist[i]=sqrt(pow((x-(*
ux)[i])/dx,2.0)+pow((y-(*
uy)[i])/dy,2.0));
296 gsl_combination *cc=gsl_combination_alloc(np,3);
298 gsl_combination_init_first(cc);
301 while (done==
false) {
302 int status=gsl_combination_next(cc);
305 i1=order[gsl_combination_get(cc,0)];
306 i2=order[gsl_combination_get(cc,1)];
307 i3=order[gsl_combination_get(cc,2)];
314 den=(-x2*y1+x3*y1+x1*y2-x3*y2-x1*y3+x2*y3);
330 gsl_combination_free(cc);
341 a=-(-y2*z1+y3*z1+y1*z2-y3*z2-y1*z3+y2*z3)/den;
342 b=-(x2*z1-x3*z1-x1*z2+x3*z2+x1*z3-x2*z3)/den;
343 c=-(x3*y2*z1-x2*y3*z1-x3*y1*z2+x1*y3*z2+x2*y1*z3-x1*y2*z3)/den;
350 #ifndef DOXYGEN_INTERNAL 366 int swap(
size_t &index_1,
double &dist_1,
size_t &index_2,
367 double &dist_2)
const {
372 index_temp=index_1; dist_temp=dist_1;
373 index_1=index_2; dist_1=dist_2;
374 index_2=index_temp; dist_2=dist_temp;
388 #ifndef DOXYGEN_NO_O2NS
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
double operator()(vec2_t &v) const
Perform the planar interpolation using the first two elements of v as input.
void set_data(size_t n_points, vec_t &x, vec_t &y, vec_t &f)
Initialize the data for the planar interpolation and compute the scaling factors. ...
void eval_points(double x, double y, double &f, size_t &i1, double &x1, double &y1, size_t &i2, double &x2, double &y2, size_t &i3, double &x3, double &y3) const
Planar interpolation returning the closest points.
invalid argument supplied by user
double y_scale
The user-specified y scale (default -1)
size_t np
The number of points.
double dy
The scale in the y direction.
double thresh
Threshold for colinearity (default )
double x_scale
The user-specified x scale (default -1)
void vector_sort_index(size_t n, const vec_t &data, vec_size_t &order)
Create a permutation which sorts the first n elements of a vector (in increasing order) ...
static const double x3[11]
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
double dx
The scale in the x direction.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
int swap(size_t &index_1, double &dist_1, size_t &index_2, double &dist_2) const
Swap points 1 and 2.
double eval(double x, double y) const
Perform the planar interpolation.
bool data_set
True if the data has been specified.
void compute_scale()
Find scaling.
Interpolate among two independent variables with planes.
static const double x2[5]
static const double x1[5]
double operator()(double x, double y) const
Perform the planar interpolation.