f4a38c8ac6e6dec23fa47eecb5d7213c39b81317
[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) 2011 - DIGITEO - Cedric DELAMARRE
4 *
5 * This file must be used under the terms of the CeCILL.
6 * This source file is licensed as described in the file COPYING, which
7 * you should have received as part of this distribution.  The terms
8 * are also available at
9 * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
10 *
11 */
12 /*--------------------------------------------------------------------------*/
13
14 #include "differential_equations_gw.hxx"
15 #include "function.hxx"
16 #include "double.hxx"
17 #include "string.hxx"
18 #include "list.hxx"
19 #include "callable.hxx"
20 #include "differentialequationfunctions.hxx"
21
22 extern "C"
23 {
24 #include "MALLOC.h"
25 #include "localization.h"
26 #include "Scierror.h"
27 #include "scifunctions.h"
28 #include "sci_warning.h"
29 #include "sciprint.h"
30 }
31
32 /*--------------------------------------------------------------------------*/
33 types::Function::ReturnValue sci_intg(types::typed_list &in, int _iRetCount, types::typed_list &out)
34 {
35     double pdA    = 0;
36     double pdB    = 0;
37     double pdEpsR = 1.0e-8;
38     double pdEpsA = 1.0e-14;
39
40     double result = 0;
41     double abserr = 0;
42
43 // *** check the minimal number of input args. ***
44     if(in.size() < 3 || in.size() > 5)
45     {
46         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "intg", 3);
47         return types::Function::Error;
48     }
49
50 // *** check number of output args ***
51     if(_iRetCount > 2)
52     {
53         Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "intg", 2);
54         return types::Function::Error;
55     }
56
57 // *** check type of input args and get it. ***
58     // A
59     if(in[0]->isDouble() == false)
60     {
61         Scierror(999, _("%s: Wrong type for input argument #%d : A matrix expected.\n"), "intg", 1);
62         return types::Function::Error;
63     }
64
65     types::Double* pDblA = in[0]->getAs<types::Double>();
66
67     if(pDblA->isScalar() == false)
68     {
69         Scierror(999, _("%s: Wrong size for input argument #%d : A scalar expected.\n"), "intg", 1);
70         return types::Function::Error;
71     }
72
73     pdA = pDblA->get(0);
74
75     // B
76     if(in[1]->isDouble() == false)
77     {
78         Scierror(999, _("%s: Wrong type for input argument #%d : A matrix expected.\n"), "intg", 2);
79         return types::Function::Error;
80     }
81
82     types::Double* pDblB = in[1]->getAs<types::Double>();
83
84     if(pDblB->isScalar() == false)
85     {
86         Scierror(999, _("%s: Wrong size for input argument #%d : A scalar expected.\n"), "intg", 2);
87         return types::Function::Error;
88     }
89
90     pdB = pDblB->get(0);
91
92     // function
93     DifferentialEquationFunctions* deFunctionsManager = new DifferentialEquationFunctions(L"intg");
94     DifferentialEquation::addDifferentialEquationFunctions(deFunctionsManager);
95
96     if(in[2]->isCallable())
97     {
98         types::Callable* pCall = in[2]->getAs<types::Callable>();
99         deFunctionsManager->setFFunction(pCall);
100
101         // check function
102         double t = 1;
103         double ret = intg_f(&t);
104         if(ret == 0)
105         {
106             Scierror(50, _("%s: Argument #%d : Variable returned by scilab argument function is incorrect.\n"), "intg", 3);
107             DifferentialEquation::removeDifferentialEquationFunctions();
108             return types::Function::Error;
109         }
110     }
111     else if(in[2]->isString())
112     {
113         bool bOK = false;
114         types::String* pStr = in[2]->getAs<types::String>();
115         bOK = deFunctionsManager->setFFunction(pStr);
116
117         if(bOK == false)
118         {
119             char* pst = wide_string_to_UTF8(pStr->get(0));
120             Scierror(50, _("%s: Subroutine not found: %s\n"), "intg", pst);
121             FREE(pst);
122             DifferentialEquation::removeDifferentialEquationFunctions();
123             return types::Function::Error;
124         }
125     }
126     else if(in[2]->isList())
127     {
128         types::List* pList = in[2]->getAs<types::List>();
129
130         if(pList->getSize() == 0)
131         {
132             Scierror(50, _("%s: Argument #%d : Subroutine not found in list: %s\n"), "intg", 3, "(string empty)");
133             DifferentialEquation::removeDifferentialEquationFunctions();
134             return types::Function::Error;
135         }
136
137         if(pList->get(0)->isCallable())
138         {
139             deFunctionsManager->setFFunction(pList->get(0)->getAs<types::Callable>());
140             for(int iter = 1; iter < pList->getSize(); iter++)
141             {
142                 deFunctionsManager->setFArgs(pList->get(iter)->getAs<types::InternalType>());
143             }
144         }
145         else
146         {
147             Scierror(999, _("%s: Wrong type for input argument #%d : The first argument in the list must be a Scilab function.\n"), "intg", 3);
148             DifferentialEquation::removeDifferentialEquationFunctions();
149             return types::Function::Error;
150         }
151     }
152     else
153     {
154         Scierror(999, _("%s: Wrong type for input argument #%d : A function expected.\n"), "intg", 3);
155         DifferentialEquation::removeDifferentialEquationFunctions();
156         return types::Function::Error;
157     }
158
159     if(in.size() > 3)
160     {
161         if(in[3]->isDouble() == false)
162         {
163             Scierror(999, _("%s: Wrong type for input argument #%d : A matrix expected.\n"), "intg", 4);
164             DifferentialEquation::removeDifferentialEquationFunctions();
165             return types::Function::Error;
166         }
167
168         types::Double* pDblEpsR = in[3]->getAs<types::Double>();
169
170         if(pDblEpsR->isScalar() == false)
171         {
172             Scierror(999, _("%s: Wrong size for input argument #%d : A scalar expected.\n"), "intg", 4);
173             DifferentialEquation::removeDifferentialEquationFunctions();
174             return types::Function::Error;
175         }
176
177         pdEpsR = pDblEpsR->get(0);
178     }
179
180     if(in.size() == 5)
181     {
182         if(in[4]->isDouble() == false)
183         {
184             Scierror(999, _("%s: Wrong type for input argument #%d : A matrix expected.\n"), "intg", 5);
185             DifferentialEquation::removeDifferentialEquationFunctions();
186             return types::Function::Error;
187         }
188
189         types::Double* pDblEpsA = in[4]->getAs<types::Double>();
190
191         if(pDblEpsA->isScalar() == false)
192         {
193             Scierror(999, _("%s: Wrong size for input argument #%d : A scalar expected.\n"), "intg", 5);
194             DifferentialEquation::removeDifferentialEquationFunctions();
195             return types::Function::Error;
196         }
197         pdEpsA = pDblEpsA->get(0);
198     }
199
200 // *** Create working table. ***
201     int limit = 750;
202
203     // rwork
204     double* alist   = (double*)malloc(limit * sizeof(double));
205     double* blist   = (double*)malloc(limit * sizeof(double));
206     double* elist   = (double*)malloc(limit * sizeof(double));
207     double* rlist   = (double*)malloc(limit * sizeof(double));
208
209     int* iwork      = (int*)malloc(limit * sizeof(int));
210
211     double epsabs   = fabs(pdEpsA);
212     double epsrel   = fabs(pdEpsR);
213
214 // *** Perform operation. ***
215     int ier = 0;
216     C2F(dqags)(intg_f, &pdA, &pdB, &epsabs, &epsrel, alist, blist, elist, rlist, &limit, iwork, &limit, &result, &abserr, &ier);
217
218     free(alist);
219     free(blist);
220     free(elist);
221     free(rlist);
222     free(iwork);
223     DifferentialEquation::removeDifferentialEquationFunctions();
224
225     if(ier)
226     {
227         switch(ier)
228         {
229             case 1 :
230             {
231                 Scierror(999, _("%s: Maximum number of subdivisions allowed has been achieved.\n"), "intg");
232                 return types::Function::Error;
233             }
234             case 2 :
235             {
236                 if(getWarningMode())
237                 {
238                     sciprint(_("%ls: Warning : The occurrence of roundoff error is detected, which prevents the requested tolerance from being achieved. The error may be under-estimated.\n"),L"intg");
239                 }
240                 break;
241             }
242             case 3 :
243             {
244                 Scierror(999, _("%s: Extremely bad integrand behaviour occurs at some points of the integration interval.\n"), "intg");
245                 return types::Function::Error;
246             }
247             case 4 :
248             {
249                 if(getWarningMode())
250                 {
251                     sciprint(_("%ls: Warning : The algorithm does not converge. Roundoff error is detected in the extrapolation table. It is presumed that the requested tolerance cannot be achieved, and that the returned result is the best which can be obtained.\n"),L"intg");
252                 }
253                 break;
254             }
255             case 5 :
256             {
257                 Scierror(999, _("%s: The integral is probably divergent, or slowly convergent.\n"), "intg");
258                 return types::Function::Error;
259             }
260         }
261     }
262
263 // *** Return result in Scilab. ***
264     types::Double* pDblOut = new types::Double(result);
265     out.push_back(pDblOut);
266
267     if(_iRetCount == 2)
268     {
269         types::Double* pDblErrOut = new types::Double(abserr);
270         out.push_back(pDblErrOut);
271     }
272
273     return types::Function::OK;
274 }
275 /*--------------------------------------------------------------------------*/
276