GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Interpolator2D.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  Steve Dennis <s.r.dennis \at liverpool.ac.uk >
7  University of Liverpool
8 */
9 //____________________________________________________________________________
10 
11 #include <cassert>
12 #include <cstdlib>
13 
15 //#include "gsl/gsl_version.h"
16 
17 using namespace genie;
18 
19 // Implemented primarily with the GSL v2 spline2d class.
20 // However, not all sites have this, so there's a backup implementation using TGraph2D
21 #if GSL_MAJOR_VERSION >= 2
22 #include "gsl/gsl_spline2d.h"
23 
24 // define our Pimpl structs
26 {
27  spline2d_container() : spl(NULL) {};
28  ~spline2d_container() { if (spl) delete spl; };
29  gsl_spline2d * spl;
30 };
31 //____________________________________________________________________________
33 {
34  interp_accel_container() : acc(NULL) {};
35  ~interp_accel_container() { if (acc) delete acc; };
36  gsl_interp_accel * acc;
37 };
38 //____________________________________________________________________________
39 
40 // Now define our actual code (GSL)
41 //____________________________________________________________________________
43  const size_t & size_x, const double * grid_x,
44  const size_t & size_y, const double * grid_y,
45  const double * knots) :
46  fSpline (new Interpolator2D::spline2d_container() ),
47  fAcc_x (new Interpolator2D::interp_accel_container() ),
48  fAcc_y (new Interpolator2D::interp_accel_container() )
49 {
50  fSpline->spl = gsl_spline2d_alloc(gsl_interp2d_bilinear,size_x,size_y);
51  gsl_spline2d_init(fSpline->spl,grid_x,grid_y,knots,size_x,size_y);
52  fAcc_x->acc = gsl_interp_accel_alloc();
53  fAcc_y->acc = gsl_interp_accel_alloc();
54 }
55 //____________________________________________________________________________
57 {
58  if (fSpline) delete fSpline;
59  if (fAcc_x ) delete fAcc_x;
60  if (fAcc_y ) delete fAcc_y;
61 }
62 //____________________________________________________________________________
63 double Interpolator2D::Eval(const double & x, const double & y) const
64 {
65  return gsl_spline2d_eval(
66  fSpline->spl, x, y,
67  fAcc_x->acc,
68  fAcc_y->acc);
69 }
70 //____________________________________________________________________________
71 double Interpolator2D::DerivX(const double & x, const double & y) const
72 {
73  return gsl_spline2d_eval_deriv_x(
74  fSpline->spl, x, y,
75  fAcc_x->acc,
76  fAcc_y->acc);
77 }
78 //____________________________________________________________________________
79 double Interpolator2D::DerivY(const double & x, const double & y) const
80 {
81  return gsl_spline2d_eval_deriv_y(
82  fSpline->spl, x, y,
83  fAcc_x->acc,
84  fAcc_y->acc);
85 }
86 //____________________________________________________________________________
87 double Interpolator2D::DerivXX(const double & x, const double & y) const
88 {
89  return gsl_spline2d_eval_deriv_xx(
90  fSpline->spl, x, y,
91  fAcc_x->acc,
92  fAcc_y->acc);
93 }
94 //____________________________________________________________________________
95 double Interpolator2D::DerivXY(const double & x, const double & y) const
96 {
97  return gsl_spline2d_eval_deriv_xy(
98  fSpline->spl, x, y,
99  fAcc_x->acc,
100  fAcc_y->acc);
101 }
102 //____________________________________________________________________________
103 double Interpolator2D::DerivYY(const double & x, const double & y) const
104 {
105  return gsl_spline2d_eval_deriv_yy(
106  fSpline->spl, x, y,
107  fAcc_x->acc,
108  fAcc_y->acc);
109 }
110 //____________________________________________________________________________
111 //____________________________________________________________________________
112 //____________________________________________________________________________
113 // And now the TGraph2D version
114 #else
115 #include "TGraph2D.h"
116 // define our Pimpl structs for TGraph
118 {
119  spline2d_container() : spl(NULL) {};
120  ~spline2d_container() { if (spl) delete spl; };
121  TGraph2D * spl;
122 };
123 //____________________________________________________________________________
125 {
128 };
129 //____________________________________________________________________________
130 // Now define our actual code (TGraph2D)
131 //____________________________________________________________________________
133  const size_t & size_x, const double * grid_x,
134  const size_t & size_y, const double * grid_y,
135  const double * knots) :
136  fSpline (new spline2d_container()),
137  fAcc_x (NULL),
138  fAcc_y (NULL)
139 {
140  fSpline->spl = new TGraph2D(size_x*size_y);
141  fSpline->spl->SetDirectory(0);
142  size_t iz = 0;
143  for (size_t iy = 0 ; iy < size_y ; iy++) {
144  for (size_t ix = 0 ; ix < size_x ; ix++) {
145  fSpline->spl->SetPoint(iz,grid_x[ix],grid_y[iy],knots[iz]);
146  iz++;
147  }
148  }
149 }
150 //____________________________________________________________________________
152 {
153  if (fSpline) delete fSpline;
154  if (fAcc_x ) delete fAcc_x;
155  if (fAcc_y ) delete fAcc_y;
156 }
157 //____________________________________________________________________________
158 double Interpolator2D::Eval(const double & x, const double & y) const
159 {
160  return fSpline->spl->Interpolate(x,y);
161 }
162 //____________________________________________________________________________
163 double Interpolator2D::DerivX(const double & x, const double & y) const
164 {
165  assert(!"Method requires GSL version 2 or higher.");
166  return -999;
167 }
168 //____________________________________________________________________________
169 double Interpolator2D::DerivY(const double & x, const double & y) const
170 {
171  assert(!"Method requires GSL version 2 or higher.");
172  return -999;
173 }
174 //____________________________________________________________________________
175 double Interpolator2D::DerivXX(const double & x, const double & y) const
176 {
177  assert(!"Method requires GSL version 2 or higher.");
178  return -999;
179 }
180 //____________________________________________________________________________
181 double Interpolator2D::DerivXY(const double & x, const double & y) const
182 {
183  assert(!"Method requires GSL version 2 or higher.");
184  return -999;
185 }
186 //____________________________________________________________________________
187 double Interpolator2D::DerivYY(const double & x, const double & y) const
188 {
189  assert(!"Method requires GSL version 2 or higher.");
190  return -999;
191 }
192 //____________________________________________________________________________
193 
194 #endif // GSL_MAJOR_VERSION
double Eval(const double &x, const double &y) const
double DerivXY(const double &x, const double &y) const
double DerivYY(const double &x, const double &y) const
double DerivXX(const double &x, const double &y) const
spline2d_container * fSpline
Interpolator2D(const size_t &size_x, const double *grid_x, const size_t &size_y, const double *grid_y, const double *knots)
double DerivX(const double &x, const double &y) const
A 2D interpolator using the GSL spline type If GSL version is not sufficient, does an inefficient ver...
interp_accel_container * fAcc_y
interp_accel_container * fAcc_x
double DerivY(const double &x, const double &y) const