1 /***************************************************************************
2     File                 : MyParser.cpp
3     Project              : QtiPlot
4     --------------------------------------------------------------------
5     Copyright            : (C) 2006 by Ion Vasilief
6     Email (use @ for *)  : ion_vasilief*yahoo.fr
7     Description          : Parser class based on muParser
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 "MyParser.h"
30 #include "muParserScripting.h"
31 #include "NonLinearFit.h"
32 
33 #include <QMessageBox>
34 #include <QApplication>
35 #include <QLocale>
36 
37 #include <gsl/gsl_const_mksa.h>
38 #include <gsl/gsl_const_num.h>
39 
MyParser()40 MyParser::MyParser()
41 :Parser()
42 {
43 	DefineConst("pi", M_PI);
44 	DefineConst("Pi", M_PI);
45 	DefineConst("PI", M_PI);
46 
47 	for (const muParserScripting::mathFunction *i=muParserScripting::math_functions; i->name; i++){
48 		if (i->numargs == 1 && i->fun1 != NULL)
49 			DefineFun(i->name, i->fun1);
50 		else if (i->numargs == 2 && i->fun2 != NULL)
51 			DefineFun(i->name, i->fun2);
52 		else if (i->numargs == 3 && i->fun3 != NULL)
53 			DefineFun(i->name, i->fun3);
54 	}
55 	gsl_set_error_handler_off();
56 
57 	setLocale(getLocale());
58 }
59 
getLocale()60 QLocale MyParser::getLocale()
61 {
62 	bool cLocale = true;
63 	foreach (QWidget *w, QApplication::allWidgets()){
64 		ApplicationWindow *app = qobject_cast<ApplicationWindow *>(w);
65 		if (app){
66 			cLocale = app->d_muparser_c_locale;
67 			break;
68 		}
69 	}
70 
71 	QLocale locale = QLocale::c();
72 	if (!cLocale)
73 		locale = QLocale();
74 
75 	return locale;
76 }
77 
setLocale(const QLocale & locale)78 void MyParser::setLocale(const QLocale& locale)
79 {
80 	const char decPoint = locale.decimalPoint().toAscii();
81 	if (decPoint != '.'){
82 		SetDecSep(decPoint);
83 		SetArgSep(';');
84 		SetThousandsSep(locale.groupSeparator().toAscii());
85 	} else
86 		ResetLocale();// reset C locale
87 }
88 
addGSLConstants()89 void MyParser::addGSLConstants()
90 {
91 	DefineConst("e", M_E);
92 	DefineConst("E", M_E);
93 
94 	//Fundamental constants provided by GSL
95 	DefineConst("c", GSL_CONST_MKSA_SPEED_OF_LIGHT);//The speed of light in vacuum
96 	DefineConst("eV", GSL_CONST_MKSA_ELECTRON_VOLT);//The energy of 1 electron volt
97 	DefineConst("g", GSL_CONST_MKSA_GRAV_ACCEL);//The standard gravitational acceleration on Earth
98 	DefineConst("G", GSL_CONST_MKSA_GRAVITATIONAL_CONSTANT);//The gravitational constant
99 	DefineConst("h", GSL_CONST_MKSA_PLANCKS_CONSTANT_H);//Planck's constant
100 	DefineConst("hbar", GSL_CONST_MKSA_PLANCKS_CONSTANT_HBAR);//Planck's constant divided by 2 pi
101 	DefineConst("k", GSL_CONST_MKSA_PLANCKS_CONSTANT_H);//The Boltzmann constant
102 	DefineConst("Na", GSL_CONST_NUM_AVOGADRO);//Avogadro's number
103 	DefineConst("R0", GSL_CONST_MKSA_MOLAR_GAS);//The molar gas constant
104 	DefineConst("V0", GSL_CONST_MKSA_STANDARD_GAS_VOLUME);//The standard gas volume
105 	DefineConst("Ry", GSL_CONST_MKSA_RYDBERG);//The Rydberg constant, Ry, in units of energy
106 }
107 
functionsList()108 const QStringList MyParser::functionsList()
109 {
110   QStringList l;
111 
112   QString argSeparator = ",";
113   if (QString(getLocale().decimalPoint()) == argSeparator)
114 	argSeparator = ";";
115 
116   for (const muParserScripting::mathFunction *i = muParserScripting::math_functions; i->name; i++){
117     if (i->numargs == 2)
118       l << QString(i->name) + "(" + argSeparator + ")";
119 	else
120       l << QString(i->name) + "()";
121   }
122   return l;
123 }
124 
functionNamesList()125 const QStringList MyParser::functionNamesList()
126 {
127 	QStringList l;
128 	for (const muParserScripting::mathFunction *i = muParserScripting::math_functions; i->name; i++)
129 		l << QString(i->name);
130 
131 	return l;
132 }
133 
explainFunction(int index)134 QString MyParser::explainFunction(int index)
135 {
136 	const muParserScripting::mathFunction i = muParserScripting::math_functions[index];
137 	QString s = QObject::tr(i.description);
138 	if (getLocale().decimalPoint() == ',')
139 		s.replace(",", ";");
140 
141 	return s;
142 }
143 
EvalRemoveSingularity(double * xvar,bool noisy) const144 double MyParser::EvalRemoveSingularity(double *xvar, bool noisy) const
145 {
146 	try {
147 		double result = Eval();
148 		if ( gsl_isinf(result) || gsl_isnan(result))
149 			throw Singularity();
150 		return result;
151 	} catch (Singularity) {
152 	    try {
153 			if (isinf(Eval()))
154 				throw Pole();
155 			double backup = *xvar;
156 			int n;
157 			frexp (*xvar, &n);
158 			double xp = *xvar + ldexp(DBL_EPSILON, n);
159 			double xm = *xvar - ldexp(DBL_EPSILON, n);
160 			*xvar = xp;
161 			double yp = Eval();
162 			if (gsl_isinf(yp) || gsl_isnan(yp))
163 				throw Pole();
164 			*xvar = xm;
165 			double ym = Eval();
166 			if (gsl_isinf(ym) || gsl_isnan(ym))
167 				throw Pole();
168 			*xvar = backup;
169 			return 0.5*(yp + ym);
170 	    } catch (Pole) {
171 	        if (noisy){
172 	        	QApplication::restoreOverrideCursor();
173 				//QMessageBox::critical(0, QObject::tr("QtiPlot - Math Error"),
174 				//QObject::tr("Found non-removable singularity at x = %1.").arg(*xvar));
175 				throw Pole();
176 			}
177 			return GSL_NAN;
178 	    }
179 	}
180 }
181 
182 //almost verbatim copy from Parser::Diff, adapted to use EvalRemoveSingularity()
183 
DiffRemoveSingularity(double * xvar,double * a_Var,double a_fPos) const184 double MyParser::DiffRemoveSingularity(double *xvar, double *a_Var, double a_fPos) const
185 {
186     double fRes(0),
187 		   fBuf(*a_Var),
188            f[4] = {0,0,0,0},
189 	       a_fEpsilon( (a_fPos == 0) ? (double)1e-10 : 1e-7 * a_fPos );
190 
191     *a_Var = a_fPos+2 * a_fEpsilon;  f[0] = EvalRemoveSingularity(xvar);
192     *a_Var = a_fPos+1 * a_fEpsilon;  f[1] = EvalRemoveSingularity(xvar);
193     *a_Var = a_fPos-1 * a_fEpsilon;  f[2] = EvalRemoveSingularity(xvar);
194     *a_Var = a_fPos-2 * a_fEpsilon;  f[3] = EvalRemoveSingularity(xvar);
195     *a_Var = fBuf; // restore variable
196 
197     fRes = (-f[0] + 8*f[1] - 8*f[2] + f[3]) / (12*a_fEpsilon);
198     return fRes;
199 }
200