← Back to branch summary

~valavanisalex/ubuntu/maverick/scidavis/scidavis-fix-565206

« back to all changes in this revision

Viewing changes to scidavis/src/Integration.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Ruben Molina
  • Date: 2009-09-06 11:34:04 UTC
  • Revision ID: james.westby@ubuntu.com-20090906113404-4awaey82l3686w4q
Tags: upstream-0.2.3
Import upstream version 0.2.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
scidavis/src/Integration.cpp
 
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 "nrutil.h"
 
31
#include "MultiLayer.h"
 
32
#include "Legend.h"
 
33
 
 
34
#include <QMessageBox>
 
35
#include <QDateTime>
 
36
#include <QLocale>
 
37
 
 
38
#include <gsl/gsl_spline.h>
 
39
#include <gsl/gsl_interp.h>
 
40
#include <gsl/gsl_vector.h>
 
41
 
 
42
Integration::Integration(ApplicationWindow *parent, Graph *g)
 
43
: Filter(parent, g)
 
44
{
 
45
        init();
 
46
}
 
47
 
 
48
Integration::Integration(ApplicationWindow *parent, Graph *g, const QString& curveTitle)
 
49
: Filter(parent, g)
 
50
{
 
51
        init();
 
52
        setDataFromCurve(curveTitle);
 
53
}
 
54
 
 
55
Integration::Integration(ApplicationWindow *parent, Graph *g, const QString& curveTitle, double start, double end)
 
56
: Filter(parent, g)
 
57
{
 
58
        init();
 
59
        setDataFromCurve(curveTitle, start, end);
 
60
}
 
61
 
 
62
void Integration::init()
 
63
{
 
64
        setName(tr("Integration"));
 
65
        d_method = 1;
 
66
    d_max_iterations = 40;
 
67
    d_sort_data = true;
 
68
}
 
69
 
 
70
QString Integration::logInfo()
 
71
{
 
72
        // Do the interpolation, use GSL libraries for that
 
73
        gsl_interp_accel *acc = gsl_interp_accel_alloc();
 
74
        const gsl_interp_type *method;
 
75
        // The method for interpolation is chosen based on the number of points
 
76
        if(d_n>3)
 
77
                method=gsl_interp_linear;
 
78
        else if(d_n>4)
 
79
                method=gsl_interp_cspline;
 
80
        else if(d_n>5)
 
81
                method=gsl_interp_akima;
 
82
 
 
83
        // If we have enough points use GSL libraries for interpolation, else use the polint algorithm
 
84
        gsl_spline *interp ;
 
85
        if(d_n>3)
 
86
        {
 
87
                interp = gsl_spline_alloc (method, d_n);
 
88
                gsl_spline_init (interp, d_x, d_y, d_n);
 
89
        }
 
90
 
 
91
        // Using Numerical Recipes
 
92
        // This is Romberg Integration method
 
93
        // This method uses the Nevilles' algorithm for interpollation;
 
94
        double yup, ylow;
 
95
        double xx,tnm,sum,del,ss,dss,error,tsum;
 
96
        if(d_n > 3)
 
97
        {
 
98
                yup = gsl_spline_eval (interp, d_to, acc);
 
99
                ylow = gsl_spline_eval (interp, d_from, acc);
 
100
        }
 
101
        else if (d_n<=3)
 
102
        {
 
103
                polint(d_x,d_y,d_n,d_to,&yup,&dss);
 
104
                polint(d_x,d_y,d_n,d_from,&ylow,&dss);
 
105
        }
 
106
 
 
107
        double *S = new double[d_max_iterations];
 
108
        double *h = new double[d_max_iterations];
 
109
        int j,it,l;
 
110
        bool success = false;
 
111
        h[0]=1.0;
 
112
        for(j=0; j < d_max_iterations; j++)
 
113
        {//Trapezoid Rule
 
114
                if(j==0)
 
115
                        S[0]=0.5*(d_to-d_from)*(ylow+yup);
 
116
                else
 
117
                {
 
118
                        h[j] = 0.25*h[j-1];
 
119
                        S[j] = S[j-1];
 
120
                        for(it=1,l=1;l<j-1;l++)it<<=1;
 
121
                        tnm=it;
 
122
                        del=(d_to-d_from)/tnm;
 
123
                        xx=d_from+0.5*del;
 
124
                        for(sum=0.0,l=1;l<=it;l++)
 
125
                        {
 
126
                                if(d_n>3)
 
127
                                        sum+=gsl_spline_eval (interp,xx, acc);
 
128
                                else if(d_n<=3)
 
129
                                {
 
130
                                        polint(d_x,d_y,d_n,xx,&tsum,&dss);
 
131
                                        sum+=tsum;
 
132
                                }
 
133
                                xx+=del;
 
134
                        }
 
135
                        S[j]=0.5*(S[j-1]+(d_to-d_from)*sum/tnm);
 
136
 
 
137
                }
 
138
                if(j>=d_method)
 
139
                {
 
140
                        polint(&h[j-d_method],&S[j-d_method],d_method,0,&ss,&dss);
 
141
                        S[j]=ss;
 
142
                }
 
143
                error=fabs(S[j]-S[j-1]);
 
144
                if(error<=d_tolerance) success = true;
 
145
                if(success) break;
 
146
        }
 
147
 
 
148
        QString logInfo = "[" + QDateTime::currentDateTime().toString(Qt::LocalDate) + "\t" + tr("Plot")+ ": ''" + d_graph->parentPlotName() + "'']\n";
 
149
        logInfo += "\n" + tr("Numerical integration of") + ": " + d_curve->title().text() + " " + tr("using a %1 order method").arg(d_method)+"\n";
 
150
        if(success)
 
151
                logInfo += tr("Iterations") + ": " + QString::number(j)+"\n";
 
152
        if(!success)
 
153
                logInfo += tr("Iterations") + ": " + QString::number(j-1)+"\n";
 
154
 
 
155
    ApplicationWindow *app = (ApplicationWindow *)parent();
 
156
    int prec = app->d_decimal_digits;
 
157
        logInfo += tr("Tolerance") + "(" + tr("max") + " = " + QLocale().toString(d_tolerance, 'g', prec)+"): " + QString::number(error)+ "\n";
 
158
        logInfo += tr("Points") + ": "+QString::number(d_n) + " " + tr("from") + " x = " +QLocale().toString(d_from, 'g', prec) + " ";
 
159
    logInfo += tr("to") + " x = " + QLocale().toString(d_to, 'g', prec) + "\n";
 
160
 
 
161
    // using GSL to find maximum value of data set
 
162
    gsl_vector *aux = gsl_vector_alloc(d_n);
 
163
    for(int i=0; i < d_n; i++)
 
164
        gsl_vector_set (aux, i, fabs(d_y[i]));
 
165
    int maxID=gsl_vector_max_index (aux);
 
166
    gsl_vector_free (aux);
 
167
 
 
168
    logInfo += tr("Peak at") + " x = " + QLocale().toString(d_x[maxID], 'g', prec)+"\t";
 
169
        logInfo += "y = " + QLocale().toString(d_y[maxID], 'g', prec)+"\n";
 
170
 
 
171
        logInfo += tr("Area") + "=";
 
172
        if(success)
 
173
                logInfo += QLocale().toString(S[j], 'g', prec);
 
174
        if(!success)
 
175
                logInfo += QLocale().toString(S[j-1], 'g', prec);
 
176
        logInfo += "\n-------------------------------------------------------------\n";
 
177
 
 
178
        if(d_n>3)
 
179
                gsl_spline_free (interp);
 
180
 
 
181
        gsl_interp_accel_free (acc);
 
182
    delete[] S;
 
183
    delete[] h;
 
184
    return logInfo;
 
185
}
 
186
 
 
187
void Integration::setMethodOrder(int n)
 
188
{
 
189
if (n < 1 || n > 5)
 
190
    {
 
191
        QMessageBox::critical((ApplicationWindow *)parent(), tr("Error"),
 
192
        tr("Unknown integration method. Valid values must be in the range: 1 (Trapezoidal Method) to 5."));
 
193
        return;
 
194
    }
 
195
 
 
196
d_method = n;
 
197
}