1 /***************************************************************************
2 File : Integration.cpp
3 Project : SciDAVis
4 --------------------------------------------------------------------
5 Copyright : (C) 2007 by Ion Vasilief
6 Email (use @ for *) : ion_vasilief*yahoo.fr
7 Description : Numerical integration of data sets
8
9 ***************************************************************************/
10
11 /***************************************************************************
12 * *
13 * This program is free software; you can redistribute it and/or modify *
14 * it under the terms of the GNU General Public License as published by *
15 * the Free Software Foundation; either version 2 of the License, or *
16 * (at your option) any later version. *
17 * *
18 * This program is distributed in the hope that it will be useful, *
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
21 * GNU General Public License for more details. *
22 * *
23 * You should have received a copy of the GNU General Public License *
24 * along with this program; if not, write to the Free Software *
25 * Foundation, Inc., 51 Franklin Street, Fifth Floor, *
26 * Boston, MA 02110-1301 USA *
27 * *
28 ***************************************************************************/
29 #include "Integration.h"
30 #include "MultiLayer.h"
31 #include "Legend.h"
32
33 #include <QMessageBox>
34 #include <QDateTime>
35 #include <QLocale>
36
37 #include <gsl/gsl_spline.h>
38 #include <gsl/gsl_interp.h>
39 #include <gsl/gsl_vector.h>
40
41 #include <stdexcept>
42 using namespace std;
43
Integration(ApplicationWindow * parent,Graph * g)44 Integration::Integration(ApplicationWindow *parent, Graph *g) : Filter(parent, g)
45 {
46 init();
47 }
48
Integration(ApplicationWindow * parent,Graph * g,const QString & curveTitle)49 Integration::Integration(ApplicationWindow *parent, Graph *g, const QString &curveTitle)
50 : Filter(parent, g)
51 {
52 init();
53 setDataFromCurve(curveTitle);
54 }
55
Integration(ApplicationWindow * parent,Graph * g,const QString & curveTitle,double start,double end)56 Integration::Integration(ApplicationWindow *parent, Graph *g, const QString &curveTitle,
57 double start, double end)
58 : Filter(parent, g)
59 {
60 init();
61 setDataFromCurve(curveTitle, start, end);
62 }
63
init()64 void Integration::init()
65 {
66 setObjectName(tr("Integration"));
67 d_method = Linear;
68 d_sort_data = true;
69 d_result = NAN;
70 }
71
trapezoid()72 double Integration::trapezoid()
73 {
74 double sum = 0.0;
75 vector<double> result;
76 result.reserve(d_n);
77 int size = d_n - 1;
78 for (int i = 0; i < size; i++) {
79 int j = i + 1;
80 result.push_back(sum);
81 sum += 0.5 * (d_y[j] + d_y[i]) * (d_x[j] - d_x[i]);
82 }
83
84 result.push_back(sum);
85 d_points = d_n;
86 addResultCurve(d_x, &result[0]);
87 return sum;
88 }
89
isDataAcceptable()90 bool Integration::isDataAcceptable()
91 {
92 switch (d_method) {
93 case Linear:
94 break;
95 case Cubic:
96 break;
97 case Akima:
98 break;
99 default:
100 QMessageBox::critical((ApplicationWindow *)parent(), tr("SciDAVis") + " - " + tr("Error"),
101 tr("Unknown interpolation method. Valid values are: 0 - Linear, 1 - "
102 "Cubic, 2 - Akima."));
103 d_init_err = true;
104 return true;
105 }
106 // GSL interpolation routines fail with division by zero on such data
107 for (unsigned i = 1; i < d_n; i++)
108 if (d_x[i - 1] == d_x[i]) {
109 QMessageBox::critical((ApplicationWindow *)parent(),
110 tr("SciDAVis") + " - " + tr("Error"),
111 tr("Several data points have the same x value causing divisions "
112 "by zero, operation aborted!"));
113 return false;
114 }
115
116 return Filter::isDataAcceptable();
117 }
118
logInfo()119 QString Integration::logInfo()
120 {
121 const gsl_interp_type *method_t;
122 QString method_name;
123 switch (d_method) {
124 case Linear:
125 method_t = gsl_interp_linear;
126 method_name = tr("Linear");
127 break;
128 case Cubic:
129 method_t = gsl_interp_cspline;
130 method_name = tr("Cubic");
131 break;
132 case Akima:
133 method_t = gsl_interp_akima;
134 method_name = tr("Akima");
135 break;
136 default:
137 throw runtime_error("invalid method");
138 }
139
140 if (d_n < gsl_interp_type_min_size(method_t)) {
141 QMessageBox::critical((ApplicationWindow *)parent(), tr("SciDAVis") + " - " + tr("Error"),
142 tr("You need at least %1 points in order to perform this operation!")
143 .arg(gsl_interp_type_min_size(method_t)));
144 d_init_err = true;
145 return QString("");
146 }
147
148 gsl_interp *interpolation = gsl_interp_alloc(method_t, d_n);
149 gsl_interp_init(interpolation, d_x, d_y, d_n);
150
151 QString logInfo = "[" + QLocale().toString(QDateTime::currentDateTime()) + "\t" + tr("Plot")
152 + ": ''" + d_graph->parentPlotName() + "'']\n";
153 logInfo += "\n" + tr("Numerical integration of") + ": " + d_curve->title().text()
154 + tr(" using ") + method_name + tr("Interpolation") + "\n";
155
156 ApplicationWindow *app = (ApplicationWindow *)parent();
157 int prec = app->d_decimal_digits;
158 logInfo += tr("Points") + ": " + QString::number(d_n) + " " + tr("from")
159 + " x = " + QLocale().toString(d_from, 'g', prec) + " ";
160 logInfo += tr("to") + " x = " + QLocale().toString(d_to, 'g', prec) + "\n";
161
162 // using GSL to find maximum value of data set
163 gsl_vector *aux = gsl_vector_alloc(d_n);
164 for (unsigned i = 0; i < d_n; i++)
165 gsl_vector_set(aux, i, d_y[i]);
166 int maxID = static_cast<int>(gsl_vector_max_index(aux));
167 gsl_vector_free(aux);
168
169 // calculate result
170 d_result = gsl_interp_eval_integ(interpolation, d_x, d_y, d_from, d_to, 0);
171
172 logInfo += tr("Peak at") + " x = " + QLocale().toString(d_x[maxID], 'g', prec) + "\t";
173 logInfo += "y = " + QLocale().toString(d_y[maxID], 'g', prec) + "\n";
174
175 logInfo += tr("Area") + "=";
176 logInfo += QLocale().toString(d_result, 'g', prec);
177 logInfo += "\n-------------------------------------------------------------\n";
178 // use trapezoid rule to get the integral curve and plot it
179 trapezoid();
180
181 gsl_interp_free(interpolation);
182
183 return logInfo;
184 }
185