source: trunk/yat/regression/GSLInterpolation.cc @ 1724

Last change on this file since 1724 was 1724, checked in by Jari Häkkinen, 12 years ago

Fixes #466. GSL_error is now thrown on GSL errors. Last thing in this ticket.

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