* Bug 16015 fixed: intg(a,b,f) called f(x) with x outside domain
[scilab.git] / scilab / modules / differential_equations / sci_gateway / cpp / sci_intg.cpp
1 /*
2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2013 - Scilab Enterprises - Cedric DELAMARRE
4 *
5  * Copyright (C) 2012 - 2016 - Scilab Enterprises
6  *
7  * This file is hereby licensed under the terms of the GNU GPL v2.0,
8  * pursuant to article 5.3.4 of the CeCILL v.2.1.
9  * This file was originally licensed under the terms of the CeCILL v2.1,
10  * and continues to be available under such terms.
11  * For more information, see the COPYING file which you should have received
12  * along with this program.
13 *
14 */
15 /*--------------------------------------------------------------------------*/
16 #include <cmath>
17
18 #include "differential_equations_gw.hxx"
19 #include "function.hxx"
20 #include "double.hxx"
21 #include "string.hxx"
22 #include "list.hxx"
23 #include "callable.hxx"
24 #include "differentialequationfunctions.hxx"
25
26 extern "C"
27 {
28 #include "sci_malloc.h"
29 #include "localization.h"
30 #include "Scierror.h"
31 #include "scifunctions.h"
32 #include "configvariable_interface.h"
33 #include "sciprint.h"
34 #include "matrix_division.h"
35 #include "vfinite.h"
36 }
37
38 /*--------------------------------------------------------------------------*/
39 types::Function::ReturnValue sci_intg(types::typed_list &in, int _iRetCount, types::typed_list &out)
40 {
41     double pdA    = 0;
42     double pdB    = 0;
43     double pdEpsR = 1.0e-8;
44     double pdEpsA = 1.0e-13;
45
46     double result = 0;
47     double abserr = 0;
48
49     int iOne = 1;
50
51     // error message catched
52     std::wostringstream os;
53     bool bCatch = false;
54
55     // *** check the minimal number of input args. ***
56     if (in.size() < 3 || in.size() > 5)
57     {
58         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "intg", 3);
59         return types::Function::Error;
60     }
61
62     // *** check number of output args ***
63     if (_iRetCount > 3)
64     {
65         Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "intg", 3);
66         return types::Function::Error;
67     }
68
69     // *** check type of input args and get it. ***
70     // A
71     if (in[0]->isDouble() == false)
72     {
73         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "intg", 1);
74         return types::Function::Error;
75     }
76
77     types::Double* pDblA = in[0]->getAs<types::Double>();
78
79     if (pDblA->isScalar() == false)
80     {
81         Scierror(999, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), "intg", 1);
82         return types::Function::Error;
83     }
84
85     pdA = pDblA->get(0);
86
87     if (ISNAN(pdA) || C2F(vfinite)(&iOne , &pdA) == false)
88     {
89         Scierror(264, _("%s: Wrong type for input argument #%d: Must not contain NaN or Inf.\n"), "intg", 1);
90         return types::Function::Error;
91     }
92
93     // B
94     if (in[1]->isDouble() == false)
95     {
96         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "intg", 2);
97         return types::Function::Error;
98     }
99
100     types::Double* pDblB = in[1]->getAs<types::Double>();
101
102     if (pDblB->isScalar() == false)
103     {
104         Scierror(999, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), "intg", 2);
105         return types::Function::Error;
106     }
107
108     pdB = pDblB->get(0);
109
110     if (ISNAN(pdB) || C2F(vfinite)(&iOne , &pdB) == false)
111     {
112         Scierror(264, _("%s: Wrong type for input argument #%d: Must not contain NaN or Inf.\n"), "intg", 1);
113         return types::Function::Error;
114     }
115
116     // function
117     DifferentialEquationFunctions deFunctionsManager(L"intg");
118     DifferentialEquation::addDifferentialEquationFunctions(&deFunctionsManager);
119
120     if (in[2]->isCallable())
121     {
122         types::Callable* pCall = in[2]->getAs<types::Callable>();
123         deFunctionsManager.setFFunction(pCall);
124
125         // check function
126         double t = pdA;
127         double ret = intg_f(&t);
128         /* if (ret == 0)
129         {
130             Scierror(50, _("%s: Argument #%d: Variable returned by scilab argument function is incorrect.\n"), "intg", 3);
131             DifferentialEquation::removeDifferentialEquationFunctions();
132             return types::Function::Error;
133         }*/
134     }
135     else if (in[2]->isString())
136     {
137         bool bOK = false;
138         types::String* pStr = in[2]->getAs<types::String>();
139         bOK = deFunctionsManager.setFFunction(pStr);
140
141         if (bOK == false)
142         {
143             char* pst = wide_string_to_UTF8(pStr->get(0));
144             Scierror(50, _("%s: Subroutine not found: %s\n"), "intg", pst);
145             FREE(pst);
146             DifferentialEquation::removeDifferentialEquationFunctions();
147             return types::Function::Error;
148         }
149     }
150     else if (in[2]->isList())
151     {
152         types::List* pList = in[2]->getAs<types::List>();
153
154         if (pList->getSize() == 0)
155         {
156             Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "intg", 3, "(string empty)");
157             DifferentialEquation::removeDifferentialEquationFunctions();
158             return types::Function::Error;
159         }
160
161         if (pList->get(0)->isCallable())
162         {
163             deFunctionsManager.setFFunction(pList->get(0)->getAs<types::Callable>());
164             for (int iter = 1; iter < pList->getSize(); iter++)
165             {
166                 deFunctionsManager.setFArgs(pList->get(iter)->getAs<types::InternalType>());
167             }
168         }
169         else
170         {
171             Scierror(999, _("%s: Wrong type for input argument #%d: The first argument in the list must be a Scilab function.\n"), "intg", 3);
172             DifferentialEquation::removeDifferentialEquationFunctions();
173             return types::Function::Error;
174         }
175     }
176     else
177     {
178         Scierror(999, _("%s: Wrong type for input argument #%d: A function expected.\n"), "intg", 3);
179         DifferentialEquation::removeDifferentialEquationFunctions();
180         return types::Function::Error;
181     }
182
183     if (in.size() > 3)
184     {
185         if (in[3]->isDouble() == false)
186         {
187             Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "intg", 4);
188             DifferentialEquation::removeDifferentialEquationFunctions();
189             return types::Function::Error;
190         }
191
192         types::Double* pDblEpsA = in[3]->getAs<types::Double>();
193
194         if (pDblEpsA->isScalar() == false)
195         {
196             Scierror(999, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), "intg", 4);
197             DifferentialEquation::removeDifferentialEquationFunctions();
198             return types::Function::Error;
199         }
200
201         pdEpsA = pDblEpsA->get(0);
202     }
203
204     if (in.size() == 5)
205     {
206         if (in[4]->isDouble() == false)
207         {
208             Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "intg", 5);
209             DifferentialEquation::removeDifferentialEquationFunctions();
210             return types::Function::Error;
211         }
212
213         types::Double* pDblEpsR = in[4]->getAs<types::Double>();
214
215         if (pDblEpsR->isScalar() == false)
216         {
217             Scierror(999, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), "intg", 5);
218             DifferentialEquation::removeDifferentialEquationFunctions();
219             return types::Function::Error;
220         }
221         pdEpsR = pDblEpsR->get(0);
222     }
223
224     // *** Create working table. ***
225     int limit   = 750;
226     int neval   = 0;
227     int last    = 0;
228     int lenw    = 4 * limit;
229
230     double* dwork   = (double*)MALLOC(lenw * sizeof(double));
231     int* iwork      = (int*)MALLOC(limit * sizeof(int));
232
233     double epsabs   = fabs(pdEpsA);
234     double epsrel   = fabs(pdEpsR);
235
236     // *** Perform operation. ***
237     int ier = 0;
238     try
239     {
240         C2F(dqags)(intg_f, &pdA, &pdB, &epsabs, &epsrel,
241                    &result, &abserr, &neval, &ier,
242                    &limit, &lenw, &last, iwork, dwork);
243     }
244     catch (ast::InternalError &ie)
245     {
246         os << ie.GetErrorMessage();
247         bCatch = true;
248     }
249
250     FREE(dwork);
251     FREE(iwork);
252     DifferentialEquation::removeDifferentialEquationFunctions();
253
254     if (bCatch)
255     {
256         wchar_t szError[bsiz];
257         os_swprintf(szError, bsiz, _W("%ls: An error occurred in '%ls' subroutine.\n").c_str(), L"intg", L"dqags");
258         os << szError;
259         throw ast::InternalError(os.str());
260     }
261
262     if (ier)
263     {
264         char* msg = NULL;
265         switch (ier)
266         {
267             case 1 :
268             {
269                 msg = _("%s: Maximum number of subdivisions achieved. Splitting the interval might help.\n");
270                 break;
271             }
272             case 2 :
273             {
274                 msg = _("%s: Round-off error detected, the requested tolerance (or default) cannot be achieved. Try using bigger tolerances.\n");
275                 break;
276             }
277             case 3 :
278             {
279                 msg = _("%s: Bad integrand behavior occurs at some points of the integration interval.\n");
280                 break;
281             }
282             case 4 :
283             {
284                 msg = _("%s: Convergence problem, round-off error detected. Try using bigger tolerances.\n");
285                 break;
286             }
287             case 5 :
288             {
289                 msg = _("%s: The integral is probably divergent, or slowly convergent.\n");
290                 break;
291             }
292             case 6 :
293             {
294                 msg = _("%s: Invalid input, absolute tolerance <= 0 and relative tolerance < 2.e-14.\n");
295                 break;
296             }
297             default :
298                 msg = _("%s: Convergence problem...\n");
299         }
300
301         if (_iRetCount == 3)
302         {
303             if (getWarningMode())
304             {
305                 sciprint(msg, "intg: Warning");
306             }
307         }
308         else
309         {
310             Scierror(999, msg, "intg: Error");
311             return types::Function::Error;
312         }
313     }
314
315     // *** Return result in Scilab. ***
316     types::Double* pDblOut = new types::Double(result);
317     out.push_back(pDblOut);
318
319     if (_iRetCount > 1)
320     {
321         out.push_back(new types::Double(abserr));
322     }
323
324     if (_iRetCount == 3)
325     {
326         out.push_back(new types::Double((double)ier));
327     }
328
329     return types::Function::OK;
330 }
331 /*--------------------------------------------------------------------------*/
332