interp2_direct.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2013-2017, 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 xi=svx.find(x);
190  size_t yi=svy.find(y);
191 
192  double xmin=(*this->xfun)[xi];
193  double xmax=(*this->xfun)[xi+1];
194  double ymin=(*this->yfun)[yi];
195  double ymax=(*this->yfun)[yi+1];
196 
197  double zminmin=(*this->datap)(xi,yi);
198  double zminmax=(*this->datap)(xi,yi+1);
199  double zmaxmin=(*this->datap)(xi+1,yi);
200  double zmaxmax=(*this->datap)(xi+1,yi+1);
201 
202  double dx=xmax-xmin;
203  double dy=ymax-ymin;
204 
205  double t=(x-xmin)/dx;
206  double u=(y-ymin)/dy;
207  double dt=1.0/dx;
208  double du=1.0/dy;
209 
210  if (itype==itp_linear) {
211  return (1.0-t)*(1.0-u)*zminmin+t*(1.0-u)*zmaxmin+
212  (1.0-t)*u*zminmax+t*u*zmaxmax;
213  }
214 
215  double zxminmin=zx(xi,yi)/dt;
216  double zxminmax=zx(xi,yi+1)/dt;
217  double zxmaxmin=zx(xi+1,yi)/dt;
218  double zxmaxmax=zx(xi+1,yi+1)/dt;
219 
220  double zyminmin=zy(xi,yi)/du;
221  double zyminmax=zy(xi,yi+1)/du;
222  double zymaxmin=zy(xi+1,yi)/du;
223  double zymaxmax=zy(xi+1,yi+1)/du;
224 
225  double zxyminmin=zxy(xi,yi)/du/dt;
226  double zxyminmax=zxy(xi,yi+1)/du/dt;
227  double zxymaxmin=zxy(xi+1,yi)/du/dt;
228  double zxymaxmax=zxy(xi+1,yi+1)/du/dt;
229 
230  double t0=1.0;
231  double t1=t;
232  double t2=t*t;
233  double t3=t*t2;
234  double u0=1.0;
235  double u1=u;
236  double u2=u*u;
237  double u3=u*u2;
238 
239  double z=0.0;
240  double v=zminmin;
241  z+=v*t0*u0;
242  v=zyminmin;
243  z+=v*t0*u1;
244  v=-3*zminmin+3*zminmax-2*zyminmin-zyminmax;
245  z+=v*t0*u2;
246  v=2*zminmin-2*zminmax+zyminmin+zyminmax;
247  z+=v*t0*u3;
248  v=zxminmin;
249  z+=v*t1*u0;
250  v=zxyminmin;
251  z+=v*t1*u1;
252  v=-3*zxminmin+3*zxminmax-2*zxyminmin-zxyminmax;
253  z+=v*t1*u2;
254  v=2*zxminmin-2*zxminmax+zxyminmin+zxyminmax;
255  z+=v*t1*u3;
256  v=-3*zminmin+3*zmaxmin-2*zxminmin-zxmaxmin;
257  z+=v*t2*u0;
258  v=-3*zyminmin+3*zymaxmin-2*zxyminmin-zxymaxmin;
259  z+=v*t2*u1;
260  v=9*zminmin-9*zmaxmin+9*zmaxmax-9*zminmax+6*zxminmin+3*zxmaxmin-
261  3*zxmaxmax-6*zxminmax+6*zyminmin-6*zymaxmin-3*zymaxmax+
262  3*zyminmax+4*zxyminmin+2*zxymaxmin+zxymaxmax+2*zxyminmax;
263  z+=v*t2*u2;
264  v=-6*zminmin+6*zmaxmin-6*zmaxmax+6*zminmax-4*zxminmin-2*zxmaxmin+
265  2*zxmaxmax+4*zxminmax-3*zyminmin+3*zymaxmin+3*zymaxmax-
266  3*zyminmax-2*zxyminmin-zxymaxmin-zxymaxmax-2*zxyminmax;
267  z+=v*t2*u3;
268  v=2*zminmin-2*zmaxmin+zxminmin+zxmaxmin;
269  z+=v*t3*u0;
270  v=2*zyminmin-2*zymaxmin+zxyminmin+zxymaxmin;
271  z+=v*t3*u1;
272  v=-6*zminmin+6*zmaxmin-6*zmaxmax+6*zminmax-3*zxminmin-3*zxmaxmin+
273  3*zxmaxmax+3*zxminmax-4*zyminmin+4*zymaxmin+2*zymaxmax-
274  2*zyminmax-2*zxyminmin-2*zxymaxmin-zxymaxmax-zxyminmax;
275  z+=v*t3*u2;
276  v=4*zminmin-4*zmaxmin+4*zmaxmax-4*zminmax+2*zxminmin+2*zxmaxmin-
277  2*zxmaxmax-2*zxminmax+2*zyminmin-2*zymaxmin-2*zymaxmax+
278  2*zyminmax+zxyminmin+zxymaxmin+zxymaxmax+zxyminmax;
279  z+=v*t3*u3;
280 
281  return z;
282  }
283 
284  /** \brief Compute the partial derivative in the x-direction
285  */
286  virtual double deriv_x(double x, double y) const {
287 
288  if (!data_set) {
289  O2SCL_ERR("Data not set in interp2_direct::deriv_x().",exc_einval);
290  }
291 
292  size_t xi=svx.find(x);
293  size_t yi=svy.find(y);
294 
295  double xmin=(*this->xfun)[xi];
296  double xmax=(*this->xfun)[xi+1];
297  double ymin=(*this->yfun)[yi];
298  double ymax=(*this->yfun)[yi+1];
299 
300  double zminmin=(*this->datap)(xi,yi);
301  double zminmax=(*this->datap)(xi,yi+1);
302  double zmaxmin=(*this->datap)(xi+1,yi);
303  double zmaxmax=(*this->datap)(xi+1,yi+1);
304 
305  double dx=xmax-xmin;
306  double dy=ymax-ymin;
307 
308  double t=(x-xmin)/dx;
309  double u=(y-ymin)/dy;
310  double dt=1.0/dx;
311  double du=1.0/dy;
312 
313  if (itype==itp_linear) {
314  return dt*(-(1.0-u)*zminmin+(1.0-u)*zmaxmin-u*zminmax+u*zmaxmax);
315  }
316 
317  double zxminmin=zx(xi,yi)/dt;
318  double zxminmax=zx(xi,yi+1)/dt;
319  double zxmaxmin=zx(xi+1,yi)/dt;
320  double zxmaxmax=zx(xi+1,yi+1)/dt;
321 
322  double zyminmin=zy(xi,yi)/du;
323  double zyminmax=zy(xi,yi+1)/du;
324  double zymaxmin=zy(xi+1,yi)/du;
325  double zymaxmax=zy(xi+1,yi+1)/du;
326 
327  double zxyminmin=zxy(xi,yi)/du/dt;
328  double zxyminmax=zxy(xi,yi+1)/du/dt;
329  double zxymaxmin=zxy(xi+1,yi)/du/dt;
330  double zxymaxmax=zxy(xi+1,yi+1)/du/dt;
331 
332  double t0=1.0;
333  double t1=t;
334  double t2=t*t;
335  double t3=t*t2;
336  double u0=1.0;
337  double u1=u;
338  double u2=u*u;
339  double u3=u*u2;
340 
341  double z_p=0.0;
342  double v=zxminmin;
343  z_p+=v*t0*u0;
344  v=zxyminmin;
345  z_p+=v*t0*u1;
346  v=-3*zxminmin+3*zxminmax-2*zxyminmin-zxyminmax;
347  z_p+=v*t0*u2;
348  v=2*zxminmin-2*zxminmax+zxyminmin+zxyminmax;
349  z_p+=v*t0*u3;
350  v=-3*zminmin+3*zmaxmin-2*zxminmin-zxmaxmin;
351  z_p+=2*v*t1*u0;
352  v=-3*zyminmin+3*zymaxmin-2*zxyminmin-zxymaxmin;
353  z_p+=2*v*t1*u1;
354  v=9*zminmin-9*zmaxmin+9*zmaxmax-9*zminmax+6*zxminmin+3*zxmaxmin-
355  3*zxmaxmax-6*zxminmax+6*zyminmin-6*zymaxmin-3*zymaxmax+
356  3*zyminmax+4*zxyminmin+2*zxymaxmin+zxymaxmax+2*zxyminmax;
357  z_p+=2*v*t1*u2;
358  v=-6*zminmin+6*zmaxmin-6*zmaxmax+6*zminmax-4*zxminmin-2*zxmaxmin+
359  2*zxmaxmax+4*zxminmax-3*zyminmin+3*zymaxmin+3*zymaxmax-
360  3*zyminmax-2*zxyminmin-zxymaxmin-zxymaxmax-2*zxyminmax;
361  z_p+=2*v*t1*u3;
362  v=2*zminmin-2*zmaxmin+zxminmin+zxmaxmin;
363  z_p+=3*v*t2*u0;
364  v=2*zyminmin-2*zymaxmin+zxyminmin+zxymaxmin;
365  z_p+=3*v*t2*u1;
366  v=-6*zminmin+6*zmaxmin-6*zmaxmax+6*zminmax-3*zxminmin-3*zxmaxmin+
367  3*zxmaxmax+3*zxminmax-4*zyminmin+4*zymaxmin+2*zymaxmax-
368  2*zyminmax-2*zxyminmin-2*zxymaxmin-zxymaxmax-zxyminmax;
369  z_p+=3*v*t2*u2;
370  v=4*zminmin-4*zmaxmin+4*zmaxmax-4*zminmax+2*zxminmin+2*zxmaxmin-
371  2*zxmaxmax-2*zxminmax+2*zyminmin-2*zymaxmin-2*zymaxmax+
372  2*zyminmax+zxyminmin+zxymaxmin+zxymaxmax+zxyminmax;
373  z_p+=3*v*t2*u3;
374  z_p*=dt;
375 
376  return z_p;
377  }
378 
379  /** \brief Compute the partial second derivative in the x-direction
380  */
381  virtual double deriv_xx(double x, double y) const {
382 
383  if (!data_set) {
384  O2SCL_ERR("Data not set in interp2_direct::deriv_xx().",exc_einval);
385  }
386 
387  if (itype==itp_linear) {
388  return 0.0;
389  }
390 
391  size_t xi=svx.find(x);
392  size_t yi=svy.find(y);
393 
394  double xmin=(*this->xfun)[xi];
395  double xmax=(*this->xfun)[xi+1];
396  double ymin=(*this->yfun)[yi];
397  double ymax=(*this->yfun)[yi+1];
398 
399  double zminmin=(*this->datap)(xi,yi);
400  double zminmax=(*this->datap)(xi,yi+1);
401  double zmaxmin=(*this->datap)(xi+1,yi);
402  double zmaxmax=(*this->datap)(xi+1,yi+1);
403 
404  double dx=xmax-xmin;
405  double dy=ymax-ymin;
406 
407  double t=(x-xmin)/dx;
408  double u=(y-ymin)/dy;
409  double dt=1.0/dx;
410  double du=1.0/dy;
411 
412  double zxminmin=zx(xi,yi)/dt;
413  double zxminmax=zx(xi,yi+1)/dt;
414  double zxmaxmin=zx(xi+1,yi)/dt;
415  double zxmaxmax=zx(xi+1,yi+1)/dt;
416 
417  double zyminmin=zy(xi,yi)/du;
418  double zyminmax=zy(xi,yi+1)/du;
419  double zymaxmin=zy(xi+1,yi)/du;
420  double zymaxmax=zy(xi+1,yi+1)/du;
421 
422  double zxyminmin=zxy(xi,yi)/du/dt;
423  double zxyminmax=zxy(xi,yi+1)/du/dt;
424  double zxymaxmin=zxy(xi+1,yi)/du/dt;
425  double zxymaxmax=zxy(xi+1,yi+1)/du/dt;
426 
427  double t0=1.0;
428  double t1=t;
429  double t2=t*t;
430  double t3=t*t2;
431  double u0=1.0;
432  double u1=u;
433  double u2=u*u;
434  double u3=u*u2;
435 
436  double z_pp=0.0;
437  double v=-3*zminmin+3*zmaxmin-2*zxminmin-zxmaxmin;
438  z_pp+=2*v*t0*u0;
439  v=-3*zyminmin+3*zymaxmin-2*zxyminmin-zxymaxmin;
440  z_pp+=2*v*t0*u1;
441  v=9*zminmin-9*zmaxmin+9*zmaxmax-9*zminmax+6*zxminmin+3*zxmaxmin-
442  3*zxmaxmax-6*zxminmax+6*zyminmin-6*zymaxmin-3*zymaxmax+
443  3*zyminmax+4*zxyminmin+2*zxymaxmin+zxymaxmax+2*zxyminmax;
444  z_pp+=2*v*t0*u2;
445  v=-6*zminmin+6*zmaxmin-6*zmaxmax+6*zminmax-4*zxminmin-2*zxmaxmin+
446  2*zxmaxmax+4*zxminmax-3*zyminmin+3*zymaxmin+3*zymaxmax-
447  3*zyminmax-2*zxyminmin-zxymaxmin-zxymaxmax-2*zxyminmax;
448  z_pp+=2*v*t0*u3;
449  v=2*zminmin-2*zmaxmin+zxminmin+zxmaxmin;
450  z_pp+=6*v*t1*u0;
451  v=2*zyminmin-2*zymaxmin+zxyminmin+zxymaxmin;
452  z_pp+=6*v*t1*u1;
453  v=-6*zminmin+6*zmaxmin-6*zmaxmax+6*zminmax-3*zxminmin-3*zxmaxmin+
454  3*zxmaxmax+3*zxminmax-4*zyminmin+4*zymaxmin+2*zymaxmax-
455  2*zyminmax-2*zxyminmin-2*zxymaxmin-zxymaxmax-zxyminmax;
456  z_pp+=6*v*t1*u2;
457  v=4*zminmin-4*zmaxmin+4*zmaxmax-4*zminmax+2*zxminmin+2*zxmaxmin-
458  2*zxmaxmax-2*zxminmax+2*zyminmin-2*zymaxmin-2*zymaxmax+
459  2*zyminmax+zxyminmin+zxymaxmin+zxymaxmax+zxyminmax;
460  z_pp+=6*v*t1*u3;
461  z_pp*=dt*dt;
462 
463  return z_pp;
464  }
465 
466  /** \brief Compute the partial derivative in the y-direction
467  */
468  virtual double deriv_y(double x, double y) const {
469 
470  if (!data_set) {
471  O2SCL_ERR("Data not set in interp2_direct::deriv_y().",exc_einval);
472  }
473 
474  size_t xi=svx.find(x);
475  size_t yi=svy.find(y);
476 
477  double xmin=(*this->xfun)[xi];
478  double xmax=(*this->xfun)[xi+1];
479  double ymin=(*this->yfun)[yi];
480  double ymax=(*this->yfun)[yi+1];
481 
482  double zminmin=(*this->datap)(xi,yi);
483  double zminmax=(*this->datap)(xi,yi+1);
484  double zmaxmin=(*this->datap)(xi+1,yi);
485  double zmaxmax=(*this->datap)(xi+1,yi+1);
486 
487  double dx=xmax-xmin;
488  double dy=ymax-ymin;
489 
490  double t=(x-xmin)/dx;
491  double u=(y-ymin)/dy;
492  double dt=1.0/dx;
493  double du=1.0/dy;
494 
495  if (itype==itp_linear) {
496  return du*(-(1.0-t)*zminmin-t*zmaxmin+(1.0-t)*zminmax+t*zmaxmax);
497  }
498 
499  double zxminmin=zx(xi,yi)/dt;
500  double zxminmax=zx(xi,yi+1)/dt;
501  double zxmaxmin=zx(xi+1,yi)/dt;
502  double zxmaxmax=zx(xi+1,yi+1)/dt;
503 
504  double zyminmin=zy(xi,yi)/du;
505  double zyminmax=zy(xi,yi+1)/du;
506  double zymaxmin=zy(xi+1,yi)/du;
507  double zymaxmax=zy(xi+1,yi+1)/du;
508 
509  double zxyminmin=zxy(xi,yi)/du/dt;
510  double zxyminmax=zxy(xi,yi+1)/du/dt;
511  double zxymaxmin=zxy(xi+1,yi)/du/dt;
512  double zxymaxmax=zxy(xi+1,yi+1)/du/dt;
513 
514  double t0=1.0;
515  double t1=t;
516  double t2=t*t;
517  double t3=t*t2;
518  double u0=1.0;
519  double u1=u;
520  double u2=u*u;
521  double u3=u*u2;
522 
523  double z_p=0.0;
524  double v=zyminmin;
525  z_p+=v*t0*u0;
526  v=-3*zminmin+3*zminmax-2*zyminmin-zyminmax;
527  z_p+=2*v*t0*u1;
528  v=2*zminmin-2*zminmax+zyminmin+zyminmax;
529  z_p+=3*v*t0*u2;
530  v=zxyminmin;
531  z_p+=v*t1*u0;
532  v=-3*zxminmin+3*zxminmax-2*zxyminmin-zxyminmax;
533  z_p+=2*v*t1*u1;
534  v=2*zxminmin-2*zxminmax+zxyminmin+zxyminmax;
535  z_p+=3*v*t1*u2;
536  v=-3*zyminmin+3*zymaxmin-2*zxyminmin-zxymaxmin;
537  z_p+=v*t2*u0;
538  v=9*zminmin-9*zmaxmin+9*zmaxmax-9*zminmax+6*zxminmin+3*zxmaxmin-
539  3*zxmaxmax-6*zxminmax+6*zyminmin-6*zymaxmin-3*zymaxmax+
540  3*zyminmax+4*zxyminmin+2*zxymaxmin+zxymaxmax+2*zxyminmax;
541  z_p+=2*v*t2*u1;
542  v=-6*zminmin+6*zmaxmin-6*zmaxmax+6*zminmax-4*zxminmin-2*zxmaxmin+
543  2*zxmaxmax+4*zxminmax-3*zyminmin+3*zymaxmin+3*zymaxmax-
544  3*zyminmax-2*zxyminmin-zxymaxmin-zxymaxmax-2*zxyminmax;
545  z_p+=3*v*t2*u2;
546  v=2*zyminmin-2*zymaxmin+zxyminmin+zxymaxmin;
547  z_p+=v*t3*u0;
548  v=-6*zminmin+6*zmaxmin-6*zmaxmax+6*zminmax-3*zxminmin-3*zxmaxmin+
549  3*zxmaxmax+3*zxminmax-4*zyminmin+4*zymaxmin+2*zymaxmax-
550  2*zyminmax-2*zxyminmin-2*zxymaxmin-zxymaxmax-zxyminmax;
551  z_p+=2*v*t3*u1;
552  v=4*zminmin-4*zmaxmin+4*zmaxmax-4*zminmax+2*zxminmin+2*zxmaxmin-
553  2*zxmaxmax-2*zxminmax+2*zyminmin-2*zymaxmin-2*zymaxmax+
554  2*zyminmax+zxyminmin+zxymaxmin+zxymaxmax+zxyminmax;
555  z_p+=3*v*t3*u2;
556  z_p*=du;
557 
558  return z_p;
559  }
560 
561  /** \brief Compute the partial second derivative in the y-direction
562  */
563  virtual double deriv_yy(double x, double y) const {
564 
565  if (!data_set) {
566  O2SCL_ERR("Data not set in interp2_direct::deriv_yy().",exc_einval);
567  }
568 
569  if (itype==itp_linear) {
570  return 0.0;
571  }
572 
573  size_t xi=svx.find(x);
574  size_t yi=svy.find(y);
575 
576  double xmin=(*this->xfun)[xi];
577  double xmax=(*this->xfun)[xi+1];
578  double ymin=(*this->yfun)[yi];
579  double ymax=(*this->yfun)[yi+1];
580 
581  double zminmin=(*this->datap)(xi,yi);
582  double zminmax=(*this->datap)(xi,yi+1);
583  double zmaxmin=(*this->datap)(xi+1,yi);
584  double zmaxmax=(*this->datap)(xi+1,yi+1);
585 
586  double dx=xmax-xmin;
587  double dy=ymax-ymin;
588 
589  double t=(x-xmin)/dx;
590  double u=(y-ymin)/dy;
591  double dt=1.0/dx;
592  double du=1.0/dy;
593 
594  double zxminmin=zx(xi,yi)/dt;
595  double zxminmax=zx(xi,yi+1)/dt;
596  double zxmaxmin=zx(xi+1,yi)/dt;
597  double zxmaxmax=zx(xi+1,yi+1)/dt;
598 
599  double zyminmin=zy(xi,yi)/du;
600  double zyminmax=zy(xi,yi+1)/du;
601  double zymaxmin=zy(xi+1,yi)/du;
602  double zymaxmax=zy(xi+1,yi+1)/du;
603 
604  double zxyminmin=zxy(xi,yi)/du/dt;
605  double zxyminmax=zxy(xi,yi+1)/du/dt;
606  double zxymaxmin=zxy(xi+1,yi)/du/dt;
607  double zxymaxmax=zxy(xi+1,yi+1)/du/dt;
608 
609  double t0=1.0;
610  double t1=t;
611  double t2=t*t;
612  double t3=t*t2;
613  double u0=1.0;
614  double u1=u;
615  double u2=u*u;
616  double u3=u*u2;
617 
618  double z_pp=0.0;
619  double v=-3*zminmin+3*zminmax-2*zyminmin-zyminmax;
620  z_pp+=2*v*t0*u0;
621  v=2*zminmin-2*zminmax+zyminmin+zyminmax;
622  z_pp+=6*v*t0*u1;
623  v=-3*zxminmin+3*zxminmax-2*zxyminmin-zxyminmax;
624  z_pp+=2*v*t1*u0;
625  v=2*zxminmin-2*zxminmax+zxyminmin+zxyminmax;
626  z_pp+=6*v*t1*u1;
627  v=9*zminmin-9*zmaxmin+9*zmaxmax-9*zminmax+6*zxminmin+3*zxmaxmin-
628  3*zxmaxmax-6*zxminmax+6*zyminmin-6*zymaxmin-3*zymaxmax+
629  3*zyminmax+4*zxyminmin+2*zxymaxmin+zxymaxmax+2*zxyminmax;
630  z_pp+=2*v*t2*u0;
631  v=-6*zminmin+6*zmaxmin-6*zmaxmax+6*zminmax-4*zxminmin-2*zxmaxmin+
632  2*zxmaxmax+4*zxminmax-3*zyminmin+3*zymaxmin+3*zymaxmax-
633  3*zyminmax-2*zxyminmin-zxymaxmin-zxymaxmax-2*zxyminmax;
634  z_pp+=6*v*t2*u1;
635  v=-6*zminmin+6*zmaxmin-6*zmaxmax+6*zminmax-3*zxminmin-3*zxmaxmin+
636  3*zxmaxmax+3*zxminmax-4*zyminmin+4*zymaxmin+2*zymaxmax-
637  2*zyminmax-2*zxyminmin-2*zxymaxmin-zxymaxmax-zxyminmax;
638  z_pp+=2*v*t3*u0;
639  v=4*zminmin-4*zmaxmin+4*zmaxmax-4*zminmax+2*zxminmin+2*zxmaxmin-
640  2*zxmaxmax-2*zxminmax+2*zyminmin-2*zymaxmin-2*zymaxmax+
641  2*zyminmax+zxyminmin+zxymaxmin+zxymaxmax+zxyminmax;
642  z_pp+=6*v*t3*u1;
643  z_pp*=du*du;
644 
645  return z_pp;
646  }
647 
648  /** \brief Compute the mixed partial derivative
649  \f$ \frac{\partial^2 f}{\partial x \partial y} \f$
650  */
651  virtual double deriv_xy(double x, double y) const {
652 
653  if (!data_set) {
654  O2SCL_ERR("Data not set in interp2_direct::deriv_xy().",exc_einval);
655  }
656 
657  size_t xi=svx.find(x);
658  size_t yi=svy.find(y);
659 
660  double xmin=(*this->xfun)[xi];
661  double xmax=(*this->xfun)[xi+1];
662  double ymin=(*this->yfun)[yi];
663  double ymax=(*this->yfun)[yi+1];
664 
665  double zminmin=(*this->datap)(xi,yi);
666  double zminmax=(*this->datap)(xi,yi+1);
667  double zmaxmin=(*this->datap)(xi+1,yi);
668  double zmaxmax=(*this->datap)(xi+1,yi+1);
669 
670  double dx=xmax-xmin;
671  double dy=ymax-ymin;
672 
673  double t=(x-xmin)/dx;
674  double u=(y-ymin)/dy;
675  double dt=1.0/dx;
676  double du=1.0/dy;
677 
678  if (itype==itp_linear) {
679  return dt*du*(zminmin-zmaxmin-zminmax+zmaxmax);
680  }
681 
682  double zxminmin=zx(xi,yi)/dt;
683  double zxminmax=zx(xi,yi+1)/dt;
684  double zxmaxmin=zx(xi+1,yi)/dt;
685  double zxmaxmax=zx(xi+1,yi+1)/dt;
686 
687  double zyminmin=zy(xi,yi)/du;
688  double zyminmax=zy(xi,yi+1)/du;
689  double zymaxmin=zy(xi+1,yi)/du;
690  double zymaxmax=zy(xi+1,yi+1)/du;
691 
692  double zxyminmin=zxy(xi,yi)/du/dt;
693  double zxyminmax=zxy(xi,yi+1)/du/dt;
694  double zxymaxmin=zxy(xi+1,yi)/du/dt;
695  double zxymaxmax=zxy(xi+1,yi+1)/du/dt;
696 
697  double t0=1.0;
698  double t1=t;
699  double t2=t*t;
700  double t3=t*t2;
701  double u0=1.0;
702  double u1=u;
703  double u2=u*u;
704  double u3=u*u2;
705 
706  double z_pp=0.0;
707  double v=zxyminmin;
708  z_pp+=v*t0*u0;
709  v=-3*zxminmin+3*zxminmax-2*zxyminmin-zxyminmax;
710  z_pp+=2*v*t0*u1;
711  v=2*zxminmin-2*zxminmax+zxyminmin+zxyminmax;
712  z_pp+=3*v*t0*u2;
713  v=-3*zyminmin+3*zymaxmin-2*zxyminmin-zxymaxmin;
714  z_pp+=2*v*t1*u0;
715  v=9*zminmin-9*zmaxmin+9*zmaxmax-9*zminmax+6*zxminmin+3*zxmaxmin-
716  3*zxmaxmax-6*zxminmax+6*zyminmin-6*zymaxmin-3*zymaxmax+
717  3*zyminmax+4*zxyminmin+2*zxymaxmin+zxymaxmax+2*zxyminmax;
718  z_pp+=4*v*t1*u1;
719  v=-6*zminmin+6*zmaxmin-6*zmaxmax+6*zminmax-4*zxminmin-2*zxmaxmin+
720  2*zxmaxmax+4*zxminmax-3*zyminmin+3*zymaxmin+3*zymaxmax-
721  3*zyminmax-2*zxyminmin-zxymaxmin-zxymaxmax-2*zxyminmax;
722  z_pp+=6*v*t1*u2;
723  v=2*zyminmin-2*zymaxmin+zxyminmin+zxymaxmin;
724  z_pp+=3*v*t2*u0;
725  v=-6*zminmin+6*zmaxmin-6*zmaxmax+6*zminmax-3*zxminmin-3*zxmaxmin+
726  3*zxmaxmax+3*zxminmax-4*zyminmin+4*zymaxmin+2*zymaxmax-
727  2*zyminmax-2*zxyminmin-2*zxymaxmin-zxymaxmax-zxyminmax;
728  z_pp+=6*v*t2*u1;
729  v=4*zminmin-4*zmaxmin+4*zmaxmax-4*zminmax+2*zxminmin+2*zxmaxmin-
730  2*zxmaxmax-2*zxminmax+2*zyminmin-2*zymaxmin-2*zymaxmax+
731  2*zyminmax+zxyminmin+zxymaxmin+zxymaxmax+zxyminmax;
732  z_pp+=9*v*t2*u2;
733  z_pp*=dt*du;
734 
735  return z_pp;
736  }
737 
738  virtual double integ_x(double x0, double x1, double y) const {
739  O2SCL_ERR("Integration unimplemented in interp2_direct.",
740  exc_eunimpl);
741  return 0.0;
742  }
743 
744  virtual double integ_y(double x, double y0, double y1) const {
745  O2SCL_ERR("Integration unimplemented in interp2_direct.",
746  exc_eunimpl);
747  return 0.0;
748  }
749 
750  virtual double eval_gen(int ix, int iy, double x0, double x1,
751  double y0, double y1) const {
752  O2SCL_ERR("Function eval_gen() unimplemented in interp2_direct.",
753  exc_eunimpl);
754  return 0.0;
755  }
756 
757 
758 #ifndef DOXYGEN_NO_O2NS
759 
760  protected:
761 
762  /// True if the data has been specified by the user
763  bool data_set;
764 
765  /// Interpolation type
766  size_t itype;
767 
768  /// Partial derivative with respect to x
769  ubmatrix zx;
770 
771  /// Partial derivative with respect to y
772  ubmatrix zy;
773 
774  /// Mixed partial derivative
775  ubmatrix zxy;
776 
777  /// Searching object for x-direction
779 
780  /// Searching object for y-direction
782 
783  private:
784 
788  (const interp2_direct<vec_t,mat_t,mat_row_t,mat_column_t>&);
789 
790 #endif
791 
792  };
793 
794 #ifndef DOXYGEN_NO_O2NS
795 }
796 #endif
797 
798 #endif
799 
800 
801 
Interpolation class for pre-specified vector.
Definition: interp.h:1685
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:72
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:74
Searching class for monotonic data with caching.
Definition: search_vec.h:79
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:70

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