source: trunk/yat/regression/GSLInterpolation.cc

Last change on this file was 2992, checked in by Peter, 8 years ago

set svndigest:ignore proprty to reflect when files were initially copied from other file (and early history should be ignored). closes #750

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
  • Property svndigest:ignore set to 1637
File size: 3.7 KB
Line 
1// $Id: GSLInterpolation.cc 2992 2013-03-03 05:03:44Z peter $
2
3/*
4  Copyright (C) 2004, 2005, 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
5  Copyright (C) 2009 Jari Häkkinen
6  Copyright (C) 2012 Peter Johansson
7
8  This file is part of the yat library, http://dev.thep.lu.se/yat
9
10  The yat library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU General Public License as
12  published by the Free Software Foundation; either version 3 of the
13  License, or (at your option) any later version.
14
15  The yat library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  General Public License for more details.
19
20  You should have received a copy of the GNU General Public License
21  along with yat. If not, see <http://www.gnu.org/licenses/>.
22*/
23
24#include <config.h>
25
26#include "GSLInterpolation.h"
27
28#include "yat/utility/Exception.h"
29#include "yat/utility/VectorBase.h"
30
31#include <cassert>
32#include <sstream>
33
34namespace theplu {
35namespace yat {
36namespace regression {
37
38
39  GSLInterpolation::GSLInterpolation(const gsl_interp_type* interpolator_type,
40                                     const utility::VectorBase& x,
41                                     const utility::VectorBase& y)
42    : evaluation_(0.0)
43  {
44    assert(x.size()==y.size());
45    size_t size=x.size();
46    x_=new double[size];
47    y_=new double[size];
48    for (size_t i=0; i<size; ++i) {
49      x_[i]=x(i);
50      y_[i]=y(i);
51    }
52    interpolator_ = gsl_interp_alloc(interpolator_type, size);
53    if (int status=gsl_interp_init(interpolator_, x_, y_, size)) {
54      std::stringstream ss;
55      ss << "GSLInterpolation::GSLInterpolation error code " << status;
56      // clean up
57      gsl_interp_free(interpolator_);
58      throw utility::GSL_error(ss.str());
59    }
60    accelerator_ = gsl_interp_accel_alloc();
61  }
62
63
64  GSLInterpolation::~GSLInterpolation(void)
65  {
66    gsl_interp_accel_free(accelerator_);
67    gsl_interp_free(interpolator_);
68    delete x_;
69    delete y_;
70  }
71
72
73  double GSLInterpolation::evaluate(double x)
74  {
75    if (int status=gsl_interp_eval_e(interpolator_, x_, y_, x, accelerator_,
76                                     &evaluation_)) {
77      std::stringstream ss;
78      ss << "GSLInterpolation::evaluate for argument " << x << " fails.\n"
79         << "   Error code " << status;
80      throw utility::GSL_error(ss.str());
81    }
82    return evaluation_;
83  }
84
85
86  double GSLInterpolation::evaluate_derivative(double x)
87  {
88    if (int status=gsl_interp_eval_deriv_e(interpolator_, x_, y_, x,
89                                           accelerator_, &evaluation_)) {
90      std::stringstream ss;
91      ss << "GSLInterpolation::evaluate_derivative for argument " << x
92         << " fails.\n   Error code " << status;
93      throw utility::GSL_error(ss.str());
94    }
95    return evaluation_;
96  }
97
98
99  double GSLInterpolation::evaluate_derivative2(double x)
100  {
101    if (int status=gsl_interp_eval_deriv2_e(interpolator_, x_, y_, x,
102                                            accelerator_, &evaluation_)) {
103      std::stringstream ss;
104      ss << "GSLInterpolation::evaluate_derivative2 for argument " << x
105         << " fails.\n   Error code " << status;
106      throw utility::GSL_error(ss.str());
107    }
108    return evaluation_;
109  }
110
111
112  double GSLInterpolation::evaluate_integral(double a, double b)
113  {
114    if (int status=gsl_interp_eval_integ_e(interpolator_, x_, y_, a, b,
115                                           accelerator_, &evaluation_)) {
116      std::stringstream ss;
117      ss << "GSLInterpolation::evaluate_integral for arguments a=" << a
118         << " and b=" << b << " fails.\n   Error code " << status;
119      throw utility::GSL_error(ss.str());
120    }
121    return evaluation_;
122  }
123
124
125  double GSLInterpolation::evaluation(void) const
126  {
127    return evaluation_;
128  }
129
130
131  unsigned int GSLInterpolation::min_size(void) const
132  {
133    return gsl_interp_min_size(interpolator_);
134  }
135
136
137}}} // of namespaces regression, yat, and theplu
Note: See TracBrowser for help on using the repository browser.