interp2_direct.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2013-2018, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl 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 General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /*
24  * Some of the code in this file was originally part of interp2d, a
25  * GSL-compatible two-dimensional interpolation library.
26  * <http://www.ellipsix.net/interp2d.html>
27  *
28  * Copyright 2012 Thomas Beutlich, David Zaslavsky
29  * Portions based on alglib 3.6 interpolation code,
30  * copyright Sergey Bochkanov
31  * Portions based on GNU GSL interpolation code,
32  * copyright 1996, 1997, 1998, 1999, 2000, 2004 Gerard Jungman
33  *
34  * This program is free software: you can redistribute it and/or modify
35  * it under the terms of the GNU General Public License as published by
36  * the Free Software Foundation, either version 3 of the License, or
37  * (at your option) any later version.
38  *
39  * This program is distributed in the hope that it will be useful,
40  * but WITHOUT ANY WARRANTY; without even the implied warranty of
41  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
42  * GNU General Public License for more details.
43  *
44  * You should have received a copy of the GNU General Public License
45  * along with this program. If not, see <http://www.gnu.org/licenses/>.
46  */
47 /** \file interp2_direct.h
48  \brief File defining \ref o2scl::interp2_direct
49 */
50 #ifndef O2SCL_INTERP2_DIRECT_H
51 #define O2SCL_INTERP2_DIRECT_H
52 
53 #include <iostream>
54 #include <string>
55 #include <vector>
56 
57 #include <boost/numeric/ublas/vector.hpp>
58 #include <boost/numeric/ublas/matrix.hpp>
59 #include <boost/numeric/ublas/matrix_proxy.hpp>
60 
61 #include <o2scl/interp.h>
62 #include <o2scl/interp2.h>
63 #include <o2scl/search_vec.h>
64 
65 #ifndef DOXYGEN_NO_O2NS
66 namespace o2scl {
67 #endif
68 
69  /** \brief Bilinear or bicubic two-dimensional interpolation
70 
71  This class implements two-dimensional interpolation. First and
72  second derivatives along both x- and y-directions can be
73  computed. This class is likely a bit faster than \ref
74  o2scl::interp2_seq but less flexible.
75 
76  The convention used by this class is that the first (row) index
77  of the matrix enumerates the x coordinate and that the second
78  (column) index enumerates the y coordinate. See the discussion
79  in the User's guide in the section called \ref rowcol_subsect.
80 
81  The function set_data() does not copy the data, it stores
82  pointers to the data. If the data is modified, then the function
83  reset_interp() must be called to reset the interpolation
84  information with the original pointer information. The storage
85  for the data, including the arrays \c x_grid and \c y_grid are
86  all managed by the user.
87 
88  By default, cubic spline interpolation with natural boundary
89  conditions is used. This can be changed by calling set_interp()
90  again with the same data and the new interpolation type.
91  Only cubic spline and linear interpolation are supported.
92 
93  Based on D. Zaslavsky's routines at
94  https://github.com/diazona/interp2d (licensed under GPLv3).
95  */
96  template<class vec_t=boost::numeric::ublas::vector<double>,
97  class mat_t=boost::numeric::ublas::matrix<double>,
98  class mat_row_t=boost::numeric::ublas::matrix_row<mat_t>,
99  class mat_column_t=boost::numeric::ublas::matrix_column<mat_t> >
100  class interp2_direct : public interp2_base<vec_t,mat_t> {
101 
102 #ifdef O2SCL_NEVER_DEFINED
103  }{
104 #endif
105 
106  public:
107 
109  typedef boost::numeric::ublas::matrix_column<ubmatrix> ubmatrix_col;
110 
111  interp2_direct() {
112  data_set=false;
114  }
115 
116  /** \brief Initialize the data for the 2-dimensional interpolation
117  */
118  void set_data(size_t n_x, size_t n_y, vec_t &x_grid,
119  vec_t &y_grid, mat_t &data,
120  size_t interp_type=itp_cspline) {
121 
122  if (interp_type!=itp_linear && interp_type!=itp_cspline &&
123  interp_type!=itp_cspline_peri) {
124  O2SCL_ERR2("Unsupported interpolation type in ",
125  "interp2_direct::set_data().",exc_eunimpl);
126  }
127 
128  this->nx=n_x;
129  this->ny=n_y;
130  this->xfun=&x_grid;
131  this->yfun=&y_grid;
132  this->datap=&data;
133  itype=interp_type;
134 
135  svx.set_vec(n_x,x_grid);
136  svy.set_vec(n_y,y_grid);
137 
138  if (interp_type==itp_cspline || interp_type==itp_cspline_peri) {
139 
140  zx.resize(n_x,n_y);
141  zy.resize(n_x,n_y);
142  zxy.resize(n_x,n_y);
143 
144  // Partial derivative with respect to x
145  for(size_t j=0;j<n_y;j++) {
146  mat_column_t col=
147  o2scl::matrix_column<mat_t,mat_column_t>(data,j);
148  interp_vec<vec_t,mat_column_t> itp(n_x,x_grid,col,interp_type);
149  for(size_t i=0;i<n_x;i++) {
150  zx(i,j)=itp.deriv(x_grid[i]);
151  }
152  }
153 
154  // Partial derivative with respect to y
155  for(size_t i=0;i<n_x;i++) {
156  mat_row_t row=
157  o2scl::matrix_row<mat_t,mat_row_t>(data,i);
158  interp_vec<vec_t,mat_row_t> itp(n_y,y_grid,row,interp_type);
159  for(size_t j=0;j<n_y;j++) {
160  zy(i,j)=itp.deriv(y_grid[j]);
161  }
162  }
163 
164  // Mixed partial derivative
165  for(size_t j=0;j<n_y;j++) {
166  ubmatrix_col col=
167  o2scl::matrix_column<ubmatrix,ubmatrix_col>(zy,j);
168  interp_vec<vec_t,ubmatrix_col> itp(n_x,x_grid,col,interp_type);
169  for(size_t i=0;i<n_x;i++) {
170  zxy(i,j)=itp.deriv(x_grid[i]);
171  }
172  }
173 
174  }
175 
176  data_set=true;
177 
178  return;
179  }
180 
181  /** \brief Perform the 2-d interpolation
182  */
183  virtual double eval(double x, double y) const {
184 
185  if (!data_set) {
186  O2SCL_ERR("Data not set in interp2_direct::eval().",exc_einval);
187  }
188 
189  size_t cache=0;
190  size_t xi=svx.find_const(x,cache);
191  cache=0;
192  size_t yi=svy.find_const(y,cache);
193 
194  double xmin=(*this->xfun)[xi];
195  double xmax=(*this->xfun)[xi+1];
196  double ymin=(*this->yfun)[yi];
197  double ymax=(*this->yfun)[yi+1];
198 
199  double zminmin=(*this->datap)(xi,yi);
200  double zminmax=(*this->datap)(xi,yi+1);
201  double zmaxmin=(*this->datap)(xi+1,yi);
202  double zmaxmax=(*this->datap)(xi+1,yi+1);
203 
204  double dx=xmax-xmin;
205  double dy=ymax-ymin;
206 
207  double t=(x-xmin)/dx;
208  double u=(y-ymin)/dy;
209  double dt=1.0/dx;
210  double du=1.0/dy;
211 
212  if (itype==itp_linear) {
213  return (1.0-t)*(1.0-u)*zminmin+t*(1.0-u)*zmaxmin+
214  (1.0-t)*u*zminmax+t*u*zmaxmax;
215  }
216 
217  double zxminmin=zx(xi,yi)/dt;
218  double zxminmax=zx(xi,yi+1)/dt;
219  double zxmaxmin=zx(xi+1,yi)/dt;
220  double zxmaxmax=zx(xi+1,yi+1)/dt;
221 
222  double zyminmin=zy(xi,yi)/du;
223  double zyminmax=zy(xi,yi+1)/du;
224  double zymaxmin=zy(xi+1,yi)/du;
225  double zymaxmax=zy(xi+1,yi+1)/du;
226 
227  double zxyminmin=zxy(xi,yi)/du/dt;
228  double zxyminmax=zxy(xi,yi+1)/du/dt;
229  double zxymaxmin=zxy(xi+1,yi)/du/dt;
230  double zxymaxmax=zxy(xi+1,yi+1)/du/dt;
231 
232  double t0=1.0;
233  double t1=t;
234  double t2=t*t;
235  double t3=t*t2;
236  double u0=1.0;
237  double u1=u;
238  double u2=u*u;
239  double u3=u*u2;
240 
241  double z=0.0;
242  double v=zminmin;
243  z+=v*t0*u0;
244  v=zyminmin;
245  z+=v*t0*u1;
246  v=-3*zminmin+3*zminmax-2*zyminmin-zyminmax;
247  z+=v*t0*u2;
248  v=2*zminmin-2*zminmax+zyminmin+zyminmax;
249  z+=v*t0*u3;
250  v=zxminmin;
251  z+=v*t1*u0;
252  v=zxyminmin;
253  z+=v*t1*u1;
254  v=-3*zxminmin+3*zxminmax-2*zxyminmin-zxyminmax;
255  z+=v*t1*u2;
256  v=2*zxminmin-2*zxminmax+zxyminmin+zxyminmax;
257  z+=v*t1*u3;
258  v=-3*zminmin+3*zmaxmin-2*zxminmin-zxmaxmin;
259  z+=v*t2*u0;
260  v=-3*zyminmin+3*zymaxmin-2*zxyminmin-zxymaxmin;
261  z+=v*t2*u1;
262  v=9*zminmin-9*zmaxmin+9*zmaxmax-9*zminmax+6*zxminmin+3*zxmaxmin-
263  3*zxmaxmax-6*zxminmax+6*zyminmin-6*zymaxmin-3*zymaxmax+
264  3*zyminmax+4*zxyminmin+2*zxymaxmin+zxymaxmax+2*zxyminmax;
265  z+=v*t2*u2;
266  v=-6*zminmin+6*zmaxmin-6*zmaxmax+6*zminmax-4*zxminmin-2*zxmaxmin+
267  2*zxmaxmax+4*zxminmax-3*zyminmin+3*zymaxmin+3*zymaxmax-
268  3*zyminmax-2*zxyminmin-zxymaxmin-zxymaxmax-2*zxyminmax;
269  z+=v*t2*u3;
270  v=2*zminmin-2*zmaxmin+zxminmin+zxmaxmin;
271  z+=v*t3*u0;
272  v=2*zyminmin-2*zymaxmin+zxyminmin+zxymaxmin;
273  z+=v*t3*u1;
274  v=-6*zminmin+6*zmaxmin-6*zmaxmax+6*zminmax-3*zxminmin-3*zxmaxmin+
275  3*zxmaxmax+3*zxminmax-4*zyminmin+4*zymaxmin+2*zymaxmax-
276  2*zyminmax-2*zxyminmin-2*zxymaxmin-zxymaxmax-zxyminmax;
277  z+=v*t3*u2;
278  v=4*zminmin-4*zmaxmin+4*zmaxmax-4*zminmax+2*zxminmin+2*zxmaxmin-
279  2*zxmaxmax-2*zxminmax+2*zyminmin-2*zymaxmin-2*zymaxmax+
280  2*zyminmax+zxyminmin+zxymaxmin+zxymaxmax+zxyminmax;
281  z+=v*t3*u3;
282 
283  return z;
284  }
285 
286  /** \brief Compute the partial derivative in the x-direction
287  */
288  virtual double deriv_x(double x, double y) const {
289 
290  if (!data_set) {
291  O2SCL_ERR("Data not set in interp2_direct::deriv_x().",exc_einval);
292  }
293 
294  size_t cache=0;
295  size_t xi=svx.find_const(x,cache);
296  cache=0;
297  size_t yi=svy.find_const(y,cache);
298 
299  double xmin=(*this->xfun)[xi];
300  double xmax=(*this->xfun)[xi+1];
301  double ymin=(*this->yfun)[yi];
302  double ymax=(*this->yfun)[yi+1];
303 
304  double zminmin=(*this->datap)(xi,yi);
305  double zminmax=(*this->datap)(xi,yi+1);
306  double zmaxmin=(*this->datap)(xi+1,yi);
307  double zmaxmax=(*this->datap)(xi+1,yi+1);
308 
309  double dx=xmax-xmin;
310  double dy=ymax-ymin;
311 
312  double t=(x-xmin)/dx;
313  double u=(y-ymin)/dy;
314  double dt=1.0/dx;
315  double du=1.0/dy;
316 
317  if (itype==itp_linear) {
318  return dt*(-(1.0-u)*zminmin+(1.0-u)*zmaxmin-u*zminmax+u*zmaxmax);
319  }
320 
321  double zxminmin=zx(xi,yi)/dt;
322  double zxminmax=zx(xi,yi+1)/dt;
323  double zxmaxmin=zx(xi+1,yi)/dt;
324  double zxmaxmax=zx(xi+1,yi+1)/dt;
325 
326  double zyminmin=zy(xi,yi)/du;
327  double zyminmax=zy(xi,yi+1)/du;
328  double zymaxmin=zy(xi+1,yi)/du;
329  double zymaxmax=zy(xi+1,yi+1)/du;
330 
331  double zxyminmin=zxy(xi,yi)/du/dt;
332  double zxyminmax=zxy(xi,yi+1)/du/dt;
333  double zxymaxmin=zxy(xi+1,yi)/du/dt;
334  double zxymaxmax=zxy(xi+1,yi+1)/du/dt;
335 
336  double t0=1.0;
337  double t1=t;
338  double t2=t*t;
339  double t3=t*t2;
340  double u0=1.0;
341  double u1=u;
342  double u2=u*u;
343  double u3=u*u2;
344 
345  double z_p=0.0;
346  double v=zxminmin;
347  z_p+=v*t0*u0;
348  v=zxyminmin;
349  z_p+=v*t0*u1;
350  v=-3*zxminmin+3*zxminmax-2*zxyminmin-zxyminmax;
351  z_p+=v*t0*u2;
352  v=2*zxminmin-2*zxminmax+zxyminmin+zxyminmax;
353  z_p+=v*t0*u3;
354  v=-3*zminmin+3*zmaxmin-2*zxminmin-zxmaxmin;
355  z_p+=2*v*t1*u0;
356  v=-3*zyminmin+3*zymaxmin-2*zxyminmin-zxymaxmin;
357  z_p+=2*v*t1*u1;
358  v=9*zminmin-9*zmaxmin+9*zmaxmax-9*zminmax+6*zxminmin+3*zxmaxmin-
359  3*zxmaxmax-6*zxminmax+6*zyminmin-6*zymaxmin-3*zymaxmax+
360  3*zyminmax+4*zxyminmin+2*zxymaxmin+zxymaxmax+2*zxyminmax;
361  z_p+=2*v*t1*u2;
362  v=-6*zminmin+6*zmaxmin-6*zmaxmax+6*zminmax-4*zxminmin-2*zxmaxmin+
363  2*zxmaxmax+4*zxminmax-3*zyminmin+3*zymaxmin+3*zymaxmax-
364  3*zyminmax-2*zxyminmin-zxymaxmin-zxymaxmax-2*zxyminmax;
365  z_p+=2*v*t1*u3;
366  v=2*zminmin-2*zmaxmin+zxminmin+zxmaxmin;
367  z_p+=3*v*t2*u0;
368  v=2*zyminmin-2*zymaxmin+zxyminmin+zxymaxmin;
369  z_p+=3*v*t2*u1;
370  v=-6*zminmin+6*zmaxmin-6*zmaxmax+6*zminmax-3*zxminmin-3*zxmaxmin+
371  3*zxmaxmax+3*zxminmax-4*zyminmin+4*zymaxmin+2*zymaxmax-
372  2*zyminmax-2*zxyminmin-2*zxymaxmin-zxymaxmax-zxyminmax;
373  z_p+=3*v*t2*u2;
374  v=4*zminmin-4*zmaxmin+4*zmaxmax-4*zminmax+2*zxminmin+2*zxmaxmin-
375  2*zxmaxmax-2*zxminmax+2*zyminmin-2*zymaxmin-2*zymaxmax+
376  2*zyminmax+zxyminmin+zxymaxmin+zxymaxmax+zxyminmax;
377  z_p+=3*v*t2*u3;
378  z_p*=dt;
379 
380  return z_p;
381  }
382 
383  /** \brief Compute the partial second derivative in the x-direction
384  */
385  virtual double deriv_xx(double x, double y) const {
386 
387  if (!data_set) {
388  O2SCL_ERR("Data not set in interp2_direct::deriv_xx().",exc_einval);
389  }
390 
391  if (itype==itp_linear) {
392  return 0.0;
393  }
394 
395  size_t cache=0;
396  size_t xi=svx.find_const(x,cache);
397  cache=0;
398  size_t yi=svy.find_const(y,cache);
399 
400  double xmin=(*this->xfun)[xi];
401  double xmax=(*this->xfun)[xi+1];
402  double ymin=(*this->yfun)[yi];
403  double ymax=(*this->yfun)[yi+1];
404 
405  double zminmin=(*this->datap)(xi,yi);
406  double zminmax=(*this->datap)(xi,yi+1);
407  double zmaxmin=(*this->datap)(xi+1,yi);
408  double zmaxmax=(*this->datap)(xi+1,yi+1);
409 
410  double dx=xmax-xmin;
411  double dy=ymax-ymin;
412 
413  double t=(x-xmin)/dx;
414  double u=(y-ymin)/dy;
415  double dt=1.0/dx;
416  double du=1.0/dy;
417 
418  double zxminmin=zx(xi,yi)/dt;
419  double zxminmax=zx(xi,yi+1)/dt;
420  double zxmaxmin=zx(xi+1,yi)/dt;
421  double zxmaxmax=zx(xi+1,yi+1)/dt;
422 
423  double zyminmin=zy(xi,yi)/du;
424  double zyminmax=zy(xi,yi+1)/du;
425  double zymaxmin=zy(xi+1,yi)/du;
426  double zymaxmax=zy(xi+1,yi+1)/du;
427 
428  double zxyminmin=zxy(xi,yi)/du/dt;
429  double zxyminmax=zxy(xi,yi+1)/du/dt;
430  double zxymaxmin=zxy(xi+1,yi)/du/dt;
431  double zxymaxmax=zxy(xi+1,yi+1)/du/dt;
432 
433  double t0=1.0;
434  double t1=t;
435  double t2=t*t;
436  double t3=t*t2;
437  double u0=1.0;
438  double u1=u;
439  double u2=u*u;
440  double u3=u*u2;
441 
442  double z_pp=0.0;
443  double v=-3*zminmin+3*zmaxmin-2*zxminmin-zxmaxmin;
444  z_pp+=2*v*t0*u0;
445  v=-3*zyminmin+3*zymaxmin-2*zxyminmin-zxymaxmin;
446  z_pp+=2*v*t0*u1;
447  v=9*zminmin-9*zmaxmin+9*zmaxmax-9*zminmax+6*zxminmin+3*zxmaxmin-
448  3*zxmaxmax-6*zxminmax+6*zyminmin-6*zymaxmin-3*zymaxmax+
449  3*zyminmax+4*zxyminmin+2*zxymaxmin+zxymaxmax+2*zxyminmax;
450  z_pp+=2*v*t0*u2;
451  v=-6*zminmin+6*zmaxmin-6*zmaxmax+6*zminmax-4*zxminmin-2*zxmaxmin+
452  2*zxmaxmax+4*zxminmax-3*zyminmin+3*zymaxmin+3*zymaxmax-
453  3*zyminmax-2*zxyminmin-zxymaxmin-zxymaxmax-2*zxyminmax;
454  z_pp+=2*v*t0*u3;
455  v=2*zminmin-2*zmaxmin+zxminmin+zxmaxmin;
456  z_pp+=6*v*t1*u0;
457  v=2*zyminmin-2*zymaxmin+zxyminmin+zxymaxmin;
458  z_pp+=6*v*t1*u1;
459  v=-6*zminmin+6*zmaxmin-6*zmaxmax+6*zminmax-3*zxminmin-3*zxmaxmin+
460  3*zxmaxmax+3*zxminmax-4*zyminmin+4*zymaxmin+2*zymaxmax-
461  2*zyminmax-2*zxyminmin-2*zxymaxmin-zxymaxmax-zxyminmax;
462  z_pp+=6*v*t1*u2;
463  v=4*zminmin-4*zmaxmin+4*zmaxmax-4*zminmax+2*zxminmin+2*zxmaxmin-
464  2*zxmaxmax-2*zxminmax+2*zyminmin-2*zymaxmin-2*zymaxmax+
465  2*zyminmax+zxyminmin+zxymaxmin+zxymaxmax+zxyminmax;
466  z_pp+=6*v*t1*u3;
467  z_pp*=dt*dt;
468 
469  return z_pp;
470  }
471 
472  /** \brief Compute the partial derivative in the y-direction
473  */
474  virtual double deriv_y(double x, double y) const {
475 
476  if (!data_set) {
477  O2SCL_ERR("Data not set in interp2_direct::deriv_y().",exc_einval);
478  }
479 
480  size_t cache=0;
481  size_t xi=svx.find_const(x,cache);
482  cache=0;
483  size_t yi=svy.find_const(y,cache);
484 
485  double xmin=(*this->xfun)[xi];
486  double xmax=(*this->xfun)[xi+1];
487  double ymin=(*this->yfun)[yi];
488  double ymax=(*this->yfun)[yi+1];
489 
490  double zminmin=(*this->datap)(xi,yi);
491  double zminmax=(*this->datap)(xi,yi+1);
492  double zmaxmin=(*this->datap)(xi+1,yi);
493  double zmaxmax=(*this->datap)(xi+1,yi+1);
494 
495  double dx=xmax-xmin;
496  double dy=ymax-ymin;
497 
498  double t=(x-xmin)/dx;
499  double u=(y-ymin)/dy;
500  double dt=1.0/dx;
501  double du=1.0/dy;
502 
503  if (itype==itp_linear) {
504  return du*(-(1.0-t)*zminmin-t*zmaxmin+(1.0-t)*zminmax+t*zmaxmax);
505  }
506 
507  double zxminmin=zx(xi,yi)/dt;
508  double zxminmax=zx(xi,yi+1)/dt;
509  double zxmaxmin=zx(xi+1,yi)/dt;
510  double zxmaxmax=zx(xi+1,yi+1)/dt;
511 
512  double zyminmin=zy(xi,yi)/du;
513  double zyminmax=zy(xi,yi+1)/du;
514  double zymaxmin=zy(xi+1,yi)/du;
515  double zymaxmax=zy(xi+1,yi+1)/du;
516 
517  double zxyminmin=zxy(xi,yi)/du/dt;
518  double zxyminmax=zxy(xi,yi+1)/du/dt;
519  double zxymaxmin=zxy(xi+1,yi)/du/dt;
520  double zxymaxmax=zxy(xi+1,yi+1)/du/dt;
521 
522  double t0=1.0;
523  double t1=t;
524  double t2=t*t;
525  double t3=t*t2;
526  double u0=1.0;
527  double u1=u;
528  double u2=u*u;
529  double u3=u*u2;
530 
531  double z_p=0.0;
532  double v=zyminmin;
533  z_p+=v*t0*u0;
534  v=-3*zminmin+3*zminmax-2*zyminmin-zyminmax;
535  z_p+=2*v*t0*u1;
536  v=2*zminmin-2*zminmax+zyminmin+zyminmax;
537  z_p+=3*v*t0*u2;
538  v=zxyminmin;
539  z_p+=v*t1*u0;
540  v=-3*zxminmin+3*zxminmax-2*zxyminmin-zxyminmax;
541  z_p+=2*v*t1*u1;
542  v=2*zxminmin-2*zxminmax+zxyminmin+zxyminmax;
543  z_p+=3*v*t1*u2;
544  v=-3*zyminmin+3*zymaxmin-2*zxyminmin-zxymaxmin;
545  z_p+=v*t2*u0;
546  v=9*zminmin-9*zmaxmin+9*zmaxmax-9*zminmax+6*zxminmin+3*zxmaxmin-
547  3*zxmaxmax-6*zxminmax+6*zyminmin-6*zymaxmin-3*zymaxmax+
548  3*zyminmax+4*zxyminmin+2*zxymaxmin+zxymaxmax+2*zxyminmax;
549  z_p+=2*v*t2*u1;
550  v=-6*zminmin+6*zmaxmin-6*zmaxmax+6*zminmax-4*zxminmin-2*zxmaxmin+
551  2*zxmaxmax+4*zxminmax-3*zyminmin+3*zymaxmin+3*zymaxmax-
552  3*zyminmax-2*zxyminmin-zxymaxmin-zxymaxmax-2*zxyminmax;
553  z_p+=3*v*t2*u2;
554  v=2*zyminmin-2*zymaxmin+zxyminmin+zxymaxmin;
555  z_p+=v*t3*u0;
556  v=-6*zminmin+6*zmaxmin-6*zmaxmax+6*zminmax-3*zxminmin-3*zxmaxmin+
557  3*zxmaxmax+3*zxminmax-4*zyminmin+4*zymaxmin+2*zymaxmax-
558  2*zyminmax-2*zxyminmin-2*zxymaxmin-zxymaxmax-zxyminmax;
559  z_p+=2*v*t3*u1;
560  v=4*zminmin-4*zmaxmin+4*zmaxmax-4*zminmax+2*zxminmin+2*zxmaxmin-
561  2*zxmaxmax-2*zxminmax+2*zyminmin-2*zymaxmin-2*zymaxmax+
562  2*zyminmax+zxyminmin+zxymaxmin+zxymaxmax+zxyminmax;
563  z_p+=3*v*t3*u2;
564  z_p*=du;
565 
566  return z_p;
567  }
568 
569  /** \brief Compute the partial second derivative in the y-direction
570  */
571  virtual double deriv_yy(double x, double y) const {
572 
573  if (!data_set) {
574  O2SCL_ERR("Data not set in interp2_direct::deriv_yy().",exc_einval);
575  }
576 
577  if (itype==itp_linear) {
578  return 0.0;
579  }
580 
581  size_t cache=0;
582  size_t xi=svx.find_const(x,cache);
583  cache=0;
584  size_t yi=svy.find_const(y,cache);
585 
586  double xmin=(*this->xfun)[xi];
587  double xmax=(*this->xfun)[xi+1];
588  double ymin=(*this->yfun)[yi];
589  double ymax=(*this->yfun)[yi+1];
590 
591  double zminmin=(*this->datap)(xi,yi);
592  double zminmax=(*this->datap)(xi,yi+1);
593  double zmaxmin=(*this->datap)(xi+1,yi);
594  double zmaxmax=(*this->datap)(xi+1,yi+1);
595 
596  double dx=xmax-xmin;
597  double dy=ymax-ymin;
598 
599  double t=(x-xmin)/dx;
600  double u=(y-ymin)/dy;
601  double dt=1.0/dx;
602  double du=1.0/dy;
603 
604  double zxminmin=zx(xi,yi)/dt;
605  double zxminmax=zx(xi,yi+1)/dt;
606  double zxmaxmin=zx(xi+1,yi)/dt;
607  double zxmaxmax=zx(xi+1,yi+1)/dt;
608 
609  double zyminmin=zy(xi,yi)/du;
610  double zyminmax=zy(xi,yi+1)/du;
611  double zymaxmin=zy(xi+1,yi)/du;
612  double zymaxmax=zy(xi+1,yi+1)/du;
613 
614  double zxyminmin=zxy(xi,yi)/du/dt;
615  double zxyminmax=zxy(xi,yi+1)/du/dt;
616  double zxymaxmin=zxy(xi+1,yi)/du/dt;
617  double zxymaxmax=zxy(xi+1,yi+1)/du/dt;
618 
619  double t0=1.0;
620  double t1=t;
621  double t2=t*t;
622  double t3=t*t2;
623  double u0=1.0;
624  double u1=u;
625  double u2=u*u;
626  double u3=u*u2;
627 
628  double z_pp=0.0;
629  double v=-3*zminmin+3*zminmax-2*zyminmin-zyminmax;
630  z_pp+=2*v*t0*u0;
631  v=2*zminmin-2*zminmax+zyminmin+zyminmax;
632  z_pp+=6*v*t0*u1;
633  v=-3*zxminmin+3*zxminmax-2*zxyminmin-zxyminmax;
634  z_pp+=2*v*t1*u0;
635  v=2*zxminmin-2*zxminmax+zxyminmin+zxyminmax;
636  z_pp+=6*v*t1*u1;
637  v=9*zminmin-9*zmaxmin+9*zmaxmax-9*zminmax+6*zxminmin+3*zxmaxmin-
638  3*zxmaxmax-6*zxminmax+6*zyminmin-6*zymaxmin-3*zymaxmax+
639  3*zyminmax+4*zxyminmin+2*zxymaxmin+zxymaxmax+2*zxyminmax;
640  z_pp+=2*v*t2*u0;
641  v=-6*zminmin+6*zmaxmin-6*zmaxmax+6*zminmax-4*zxminmin-2*zxmaxmin+
642  2*zxmaxmax+4*zxminmax-3*zyminmin+3*zymaxmin+3*zymaxmax-
643  3*zyminmax-2*zxyminmin-zxymaxmin-zxymaxmax-2*zxyminmax;
644  z_pp+=6*v*t2*u1;
645  v=-6*zminmin+6*zmaxmin-6*zmaxmax+6*zminmax-3*zxminmin-3*zxmaxmin+
646  3*zxmaxmax+3*zxminmax-4*zyminmin+4*zymaxmin+2*zymaxmax-
647  2*zyminmax-2*zxyminmin-2*zxymaxmin-zxymaxmax-zxyminmax;
648  z_pp+=2*v*t3*u0;
649  v=4*zminmin-4*zmaxmin+4*zmaxmax-4*zminmax+2*zxminmin+2*zxmaxmin-
650  2*zxmaxmax-2*zxminmax+2*zyminmin-2*zymaxmin-2*zymaxmax+
651  2*zyminmax+zxyminmin+zxymaxmin+zxymaxmax+zxyminmax;
652  z_pp+=6*v*t3*u1;
653  z_pp*=du*du;
654 
655  return z_pp;
656  }
657 
658  /** \brief Compute the mixed partial derivative
659  \f$ \frac{\partial^2 f}{\partial x \partial y} \f$
660  */
661  virtual double deriv_xy(double x, double y) const {
662 
663  if (!data_set) {
664  O2SCL_ERR("Data not set in interp2_direct::deriv_xy().",exc_einval);
665  }
666 
667  size_t cache=0;
668  size_t xi=svx.find_const(x,cache);
669  cache=0;
670  size_t yi=svy.find_const(y,cache);
671 
672  double xmin=(*this->xfun)[xi];
673  double xmax=(*this->xfun)[xi+1];
674  double ymin=(*this->yfun)[yi];
675  double ymax=(*this->yfun)[yi+1];
676 
677  double zminmin=(*this->datap)(xi,yi);
678  double zminmax=(*this->datap)(xi,yi+1);
679  double zmaxmin=(*this->datap)(xi+1,yi);
680  double zmaxmax=(*this->datap)(xi+1,yi+1);
681 
682  double dx=xmax-xmin;
683  double dy=ymax-ymin;
684 
685  double t=(x-xmin)/dx;
686  double u=(y-ymin)/dy;
687  double dt=1.0/dx;
688  double du=1.0/dy;
689 
690  if (itype==itp_linear) {
691  return dt*du*(zminmin-zmaxmin-zminmax+zmaxmax);
692  }
693 
694  double zxminmin=zx(xi,yi)/dt;
695  double zxminmax=zx(xi,yi+1)/dt;
696  double zxmaxmin=zx(xi+1,yi)/dt;
697  double zxmaxmax=zx(xi+1,yi+1)/dt;
698 
699  double zyminmin=zy(xi,yi)/du;
700  double zyminmax=zy(xi,yi+1)/du;
701  double zymaxmin=zy(xi+1,yi)/du;
702  double zymaxmax=zy(xi+1,yi+1)/du;
703 
704  double zxyminmin=zxy(xi,yi)/du/dt;
705  double zxyminmax=zxy(xi,yi+1)/du/dt;
706  double zxymaxmin=zxy(xi+1,yi)/du/dt;
707  double zxymaxmax=zxy(xi+1,yi+1)/du/dt;
708 
709  double t0=1.0;
710  double t1=t;
711  double t2=t*t;
712  double t3=t*t2;
713  double u0=1.0;
714  double u1=u;
715  double u2=u*u;
716  double u3=u*u2;
717 
718  double z_pp=0.0;
719  double v=zxyminmin;
720  z_pp+=v*t0*u0;
721  v=-3*zxminmin+3*zxminmax-2*zxyminmin-zxyminmax;
722  z_pp+=2*v*t0*u1;
723  v=2*zxminmin-2*zxminmax+zxyminmin+zxyminmax;
724  z_pp+=3*v*t0*u2;
725  v=-3*zyminmin+3*zymaxmin-2*zxyminmin-zxymaxmin;
726  z_pp+=2*v*t1*u0;
727  v=9*zminmin-9*zmaxmin+9*zmaxmax-9*zminmax+6*zxminmin+3*zxmaxmin-
728  3*zxmaxmax-6*zxminmax+6*zyminmin-6*zymaxmin-3*zymaxmax+
729  3*zyminmax+4*zxyminmin+2*zxymaxmin+zxymaxmax+2*zxyminmax;
730  z_pp+=4*v*t1*u1;
731  v=-6*zminmin+6*zmaxmin-6*zmaxmax+6*zminmax-4*zxminmin-2*zxmaxmin+
732  2*zxmaxmax+4*zxminmax-3*zyminmin+3*zymaxmin+3*zymaxmax-
733  3*zyminmax-2*zxyminmin-zxymaxmin-zxymaxmax-2*zxyminmax;
734  z_pp+=6*v*t1*u2;
735  v=2*zyminmin-2*zymaxmin+zxyminmin+zxymaxmin;
736  z_pp+=3*v*t2*u0;
737  v=-6*zminmin+6*zmaxmin-6*zmaxmax+6*zminmax-3*zxminmin-3*zxmaxmin+
738  3*zxmaxmax+3*zxminmax-4*zyminmin+4*zymaxmin+2*zymaxmax-
739  2*zyminmax-2*zxyminmin-2*zxymaxmin-zxymaxmax-zxyminmax;
740  z_pp+=6*v*t2*u1;
741  v=4*zminmin-4*zmaxmin+4*zmaxmax-4*zminmax+2*zxminmin+2*zxmaxmin-
742  2*zxmaxmax-2*zxminmax+2*zyminmin-2*zymaxmin-2*zymaxmax+
743  2*zyminmax+zxyminmin+zxymaxmin+zxymaxmax+zxyminmax;
744  z_pp+=9*v*t2*u2;
745  z_pp*=dt*du;
746 
747  return z_pp;
748  }
749 
750  virtual double integ_x(double x0, double x1, double y) const {
751  O2SCL_ERR("Integration unimplemented in interp2_direct.",
752  exc_eunimpl);
753  return 0.0;
754  }
755 
756  virtual double integ_y(double x, double y0, double y1) const {
757  O2SCL_ERR("Integration unimplemented in interp2_direct.",
758  exc_eunimpl);
759  return 0.0;
760  }
761 
762  virtual double eval_gen(int ix, int iy, double x0, double x1,
763  double y0, double y1) const {
764  O2SCL_ERR("Function eval_gen() unimplemented in interp2_direct.",
765  exc_eunimpl);
766  return 0.0;
767  }
768 
769 
770 #ifndef DOXYGEN_NO_O2NS
771 
772  protected:
773 
774  /// True if the data has been specified by the user
775  bool data_set;
776 
777  /// Interpolation type
778  size_t itype;
779 
780  /// Partial derivative with respect to x
781  ubmatrix zx;
782 
783  /// Partial derivative with respect to y
784  ubmatrix zy;
785 
786  /// Mixed partial derivative
787  ubmatrix zxy;
788 
789  /// Searching object for x-direction
791 
792  /// Searching object for y-direction
794 
795  private:
796 
800  (const interp2_direct<vec_t,mat_t,mat_row_t,mat_column_t>&);
801 
802 #endif
803 
804  };
805 
806 #ifndef DOXYGEN_NO_O2NS
807 }
808 #endif
809 
810 #endif
811 
812 
813 
Interpolation class for pre-specified vector.
Definition: interp.h:1795
virtual double deriv_y(double x, double y) const
Compute the partial derivative in the y-direction.
ubmatrix zx
Partial derivative with respect to x.
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
virtual double deriv_x(double x, double y) const
Compute the partial derivative in the x-direction.
size_t ny
The number of y grid points.
Definition: interp2.h:128
invalid argument supplied by user
Definition: err_hnd.h:59
bool data_set
True if the data has been specified by the user.
virtual double eval_gen(int ix, int iy, double x0, double x1, double y0, double y1) const
Compute a general interpolation result.
virtual double integ_x(double x0, double x1, double y) const
Compute the integral in the x-direction between x=x0 and x=x1.
virtual double deriv_xx(double x, double y) const
Compute the partial second derivative in the x-direction.
requested feature not (yet) implemented
Definition: err_hnd.h:99
size_t itype
Interpolation type.
virtual double integ_y(double x, double y0, double y1) const
Compute the integral in the y-direction between y=y0 and y=y1.
Cubic spline for natural boundary conditions.
Definition: interp.h:73
ubmatrix zxy
Mixed partial derivative.
mat_t * datap
The data.
Definition: interp2.h:137
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
search_vec< vec_t > svy
Searching object for y-direction.
virtual double eval(double x, double y) const
Perform the 2-d interpolation.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
size_t nx
The number of x grid points.
Definition: interp2.h:125
void set_data(size_t n_x, size_t n_y, vec_t &x_grid, vec_t &y_grid, mat_t &data, size_t interp_type=itp_cspline)
Initialize the data for the 2-dimensional interpolation.
vec_t * xfun
The x grid.
Definition: interp2.h:131
virtual double deriv_xy(double x, double y) const
Compute the mixed partial derivative .
vec_t * yfun
The y grid.
Definition: interp2.h:134
Two-dimensional interpolation base class [abstract].
Definition: interp2.h:40
Cubic spline for periodic boundary conditions.
Definition: interp.h:75
Searching class for monotonic data with caching.
Definition: search_vec.h:74
Bilinear or bicubic two-dimensional interpolation.
ubmatrix zy
Partial derivative with respect to y.
virtual double deriv_yy(double x, double y) const
Compute the partial second derivative in the y-direction.
static const double x1[5]
Definition: inte_qng_gsl.h:48
search_vec< vec_t > svx
Searching object for x-direction.
Linear.
Definition: interp.h:71

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).