GammaLib  2.1.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GMath.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GMath.cpp - Mathematical functions *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2012-2018 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GMath.cpp
23  * @brief Mathematical function implementations
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cstdlib>
32 #include <cmath>
33 #include "GMath.hpp"
34 #include "GTools.hpp"
35 
36 /* __ Method name definitions ____________________________________________ */
37 
38 /* __ Macros _____________________________________________________________ */
39 
40 /* __ Coding definitions _________________________________________________ */
41 //#define G_USE_ASIN_FOR_ACOS //!< Use asin for acos computations
42 
43 /* __ Debug definitions __________________________________________________ */
44 
45 
46 /*==========================================================================
47  = =
48  = Functions =
49  = =
50  ==========================================================================*/
51 
52 /***********************************************************************//**
53  * @brief Computes acos by avoiding NaN due to rounding errors
54  *
55  * @param[in] arg Argument.
56  * @return Arc cosine of argument.
57  *
58  * Returns the arc cosine by restricting the argument to [-1,1].
59  *
60  * If the compile option G_USE_ASIN_FOR_ACOS is defined, the function will
61  * compute
62  *
63  * \f[
64  * acos(x) = \frac{\pi}{2} - asin(x)
65  * \f]
66  *
67  * which happens to be faster on most systems.
68  ***************************************************************************/
69 double gammalib::acos(const double& arg)
70 {
71  // Allocate result
72  double arccos;
73 
74  // Compute acos
75  if (arg >= 1) {
76  arccos = 0.0;
77  }
78  else if (arg <= -1.0) {
79  arccos = gammalib::pi;
80  }
81  else {
82  #if defined(G_USE_ASIN_FOR_ACOS)
83  arccos = gammalib::pihalf - std::asin(arg);
84  #else
85  arccos = std::acos(arg);
86  #endif
87  }
88 
89  // Return result
90  return arccos;
91 }
92 
93 
94 /***********************************************************************//**
95  * @brief Compute arc tangens in radians
96  *
97  * @param[in] y Nominator
98  * @param[in] x Denominator
99  *
100  * This code has been adapted from the WCSLIB function wcstrig.c::atan2d().
101  ***************************************************************************/
102 double gammalib::atan2(const double& y, const double& x)
103 {
104  // Check for rounding errors
105  if (y == 0.0) {
106  if (x >= 0.0) {
107  return 0.0;
108  }
109  else if (x < 0.0) {
110  return pi;
111  }
112  }
113  else if (x == 0.0) {
114  if (y > 0.0) {
115  return pihalf;
116  }
117  else if (y < 0.0) {
118  return -pihalf;
119  }
120  }
121 
122  // Return arc sine
123  return (std::atan2(y,x));
124 }
125 
126 
127 /***********************************************************************//**
128  * @brief Compute cosine of angle in degrees
129  *
130  * @param[in] angle Angle in degrees
131  *
132  * This code has been adapted from the WCSLIB function wcstrig.c::cosd().
133  ***************************************************************************/
134 double gammalib::cosd(const double& angle)
135 {
136  // Check for rounding errors
137  if (fmod(angle, 90.0) == 0.0) {
138  int i = std::abs((int)std::floor(angle/90.0 + 0.5)) % 4;
139  switch (i) {
140  case 0:
141  return 1.0;
142  case 1:
143  return 0.0;
144  case 2:
145  return -1.0;
146  case 3:
147  return 0.0;
148  }
149  }
150 
151  // Return cosine
152  return std::cos(angle * gammalib::deg2rad);
153 }
154 
155 
156 /***********************************************************************//**
157  * @brief Compute sine of angle in degrees
158  *
159  * @param[in] angle Angle in degrees
160  *
161  * This code has been adapted from the WCSLIB function wcstrig.c::sind().
162  ***************************************************************************/
163 double gammalib::sind(const double& angle)
164 {
165  // Check for rounding errors
166  if (fmod(angle, 90.0) == 0.0) {
167  int i = std::abs((int)std::floor(angle/90.0 - 0.5))%4;
168  switch (i) {
169  case 0:
170  return 1.0;
171  case 1:
172  return 0.0;
173  case 2:
174  return -1.0;
175  case 3:
176  return 0.0;
177  }
178  }
179 
180  // Return sine
181  return std::sin(angle * gammalib::deg2rad);
182 }
183 
184 
185 /***********************************************************************//**
186  * @brief Compute tangens of angle in degrees
187  *
188  * @param[in] angle Angle in degrees
189  *
190  * This code has been adapted from the WCSLIB function wcstrig.c::tand().
191  ***************************************************************************/
192 double gammalib::tand(const double& angle)
193 {
194  // Check for rounding errors
195  double resid = fmod(angle, 360.0);
196  if (resid == 0.0 || std::abs(resid) == 180.0) {
197  return 0.0;
198  }
199  else if (resid == 45.0 || resid == 225.0) {
200  return 1.0;
201  }
202  else if (resid == -135.0 || resid == -315.0) {
203  return -1.0;
204  }
205 
206  // Return tangens
207  return std::tan(angle * gammalib::deg2rad);
208 }
209 
210 
211 /***********************************************************************//**
212  * @brief Compute arc cosine in degrees
213  *
214  * @param[in] value Value
215  *
216  * This code has been adapted from the WCSLIB function wcstrig.c::acosd().
217  ***************************************************************************/
218 double gammalib::acosd(const double& value)
219 {
220  // Domain tolerance
221  const double wcstrig_tol = 1.0e-10;
222 
223  // Check for rounding errors
224  if (value >= 1.0) {
225  if (value-1.0 < wcstrig_tol) {
226  return 0.0;
227  }
228  }
229  else if (value == 0.0) {
230  return 90.0;
231  }
232  else if (value <= -1.0) {
233  if (value+1.0 > -wcstrig_tol) {
234  return 180.0;
235  }
236  }
237 
238  // Return arc cosine
239  return std::acos(value) * gammalib::rad2deg;
240 }
241 
242 
243 /***********************************************************************//**
244  * @brief Compute arc sine in degrees
245  *
246  * @param[in] value Value
247  *
248  * This code has been adapted from the WCSLIB function wcstrig.c::asind().
249  ***************************************************************************/
250 double gammalib::asind(const double& value)
251 {
252  // Domain tolerance
253  const double wcstrig_tol = 1.0e-10;
254 
255  // Check for rounding errors
256  if (value <= -1.0) {
257  if (value+1.0 > -wcstrig_tol) {
258  return -90.0;
259  }
260  }
261  else if (value == 0.0) {
262  return 0.0;
263  }
264  else if (value >= 1.0) {
265  if (value-1.0 < wcstrig_tol) {
266  return 90.0;
267  }
268  }
269 
270  // Return arc sine
271  return std::asin(value) * gammalib::rad2deg;
272 }
273 
274 
275 /***********************************************************************//**
276  * @brief Compute arc tangens in degrees
277  *
278  * @param[in] value Value
279  *
280  * This code has been adapted from the WCSLIB function wcstrig.c::atand().
281  ***************************************************************************/
282 double gammalib::atand(const double& value)
283 {
284  // Check for rounding errors
285  if (value == -1.0) {
286  return -45.0;
287  }
288  else if (value == 0.0) {
289  return 0.0;
290  }
291  else if (value == 1.0) {
292  return 45.0;
293  }
294 
295  // Return arc sine
296  return std::atan(value) * gammalib::rad2deg;
297 }
298 
299 
300 /***********************************************************************//**
301  * @brief Compute arc tangens in degrees
302  *
303  * @param[in] y Nominator
304  * @param[in] x Denominator
305  *
306  * This code has been adapted from the WCSLIB function wcstrig.c::atan2d().
307  ***************************************************************************/
308 double gammalib::atan2d(const double& y, const double& x)
309 {
310  // Check for rounding errors
311  if (y == 0.0) {
312  if (x >= 0.0) {
313  return 0.0;
314  }
315  else if (x < 0.0) {
316  return 180.0;
317  }
318  }
319  else if (x == 0.0) {
320  if (y > 0.0) {
321  return 90.0;
322  }
323  else if (y < 0.0) {
324  return -90.0;
325  }
326  }
327 
328  // Return arc sine
329  return std::atan2(y,x) * gammalib::rad2deg;
330 }
331 
332 
333 /***********************************************************************//**
334  * @brief Compute sine and cosine of angle in degrees
335  *
336  * @param[in] angle Angle [degrees].
337  * @param[out] s Sine of angle.
338  * @param[out] c Cosine of angle.
339  *
340  * This code has been adapted from the WCSLIB function wcstrig.c::sincosd().
341  ***************************************************************************/
342 void gammalib::sincosd(const double& angle, double *s, double *c)
343 
344 {
345  // Check for rounding errors
346  if (fmod(angle, 90.0) == 0.0) {
347  int i = std::abs((int)std::floor(angle/90.0 + 0.5)) % 4;
348  switch (i) {
349  case 0:
350  *s = 0.0;
351  *c = 1.0;
352  return;
353  case 1:
354  *s = (angle > 0.0) ? 1.0 : -1.0;
355  *c = 0.0;
356  return;
357  case 2:
358  *s = 0.0;
359  *c = -1.0;
360  return;
361  case 3:
362  *s = (angle > 0.0) ? -1.0 : 1.0;
363  *c = 0.0;
364  return;
365  }
366  }
367 
368  // Compute sine and cosine
369  *s = std::sin(angle * gammalib::deg2rad);
370  *c = std::cos(angle * gammalib::deg2rad);
371 
372  // Return
373  return;
374 }
375 
376 
377 /***********************************************************************//**
378  * @brief Computes logarithm of gamma function
379  *
380  * @param[in] arg Argument.
381  * @return Logarithm of gamma function.
382  ***************************************************************************/
383 double gammalib::gammln(const double& arg) {
384 
385  // Define static constants
386  static const double cof[6] = { 76.18009172947146,
387  -86.50532032941677,
388  24.01409824083091,
389  -1.231739572450155,
390  0.1208650973866179e-2,
391  -0.5395239384953e-5};
392 
393  // Evaluate logarithm of gamma function
394  double a = arg;
395  double b = arg;
396  double c = a + 5.5;
397  c -= (a + 0.5) * std::log(c);
398  double d = 1.000000000190015;
399  for (int i = 0; i < 6; ++i) {
400  d += cof[i]/++b;
401  }
402  double result = std::log(2.5066282746310005 * d/a) - c;
403 
404  // Return result
405  return result;
406 }
407 
408 
409 /***********************************************************************//**
410  * @brief Computes error function
411  *
412  * @param[in] arg Argument.
413  * @return Error function.
414  *
415  * Reference: http://en.wikipedia.org/wiki/Complementary_error_function
416  ***************************************************************************/
417 double gammalib::erf(const double& arg)
418 {
419  // Initialise
420  double z = std::abs(arg);
421  double t = 1.0/(1.0+0.5*z);
422 
423  // Compute tau
424  double tau = t*std::exp(-z*z - 1.26551223 + t *
425  ( 1.00002368 + t *
426  ( 0.37409196 + t *
427  ( 0.09678418 + t *
428  (-0.18628806 + t *
429  ( 0.27886807 + t *
430  (-1.13520398 + t *
431  ( 1.48851587 + t *
432  (-0.82215223 + t * 0.17087277)))))))));
433 
434  // Compute result
435  double result = (arg >= 0.0) ? 1.0 - tau : tau - 1.0;
436 
437  // Return result
438  return result;
439 }
440 
441 
442 /***********************************************************************//**
443  * @brief Computes complementary error function
444  *
445  * @param[in] arg Argument.
446  * @return Complementary error function.
447  *
448  * Reference: http://en.wikipedia.org/wiki/Complementary_error_function
449  ***************************************************************************/
450 double gammalib::erfc(const double& arg)
451 {
452  // Return result
453  return (1.0 - erf(arg));
454 }
455 
456 
457 /***********************************************************************//**
458  * @brief Computes inverse error function
459  *
460  * @param[in] arg Argument.
461  * @return Inverse error function.
462  ***************************************************************************/
463 double gammalib::erfinv(const double& arg)
464 {
465  // Define constants
466  static const double a[4] = { 0.886226899, -1.645349621, 0.914624893,
467  -0.140543331};
468  static const double b[4] = {-2.118377725, 1.442710462, -0.329097515,
469  0.012229801};
470  static const double c[4] = {-1.970840454, -1.624906493, 3.429567803,
471  1.641345311};
472  static const double d[2] = { 3.543889200, 1.637067800};
473  static const double e = 2.0/std::sqrt(gammalib::pi);
474 
475  // Allocate result
476  double result;
477 
478  // Return NAN if out of range
479  if (std::abs(arg) > 1.0) {
480  result = std::atof("NAN");
481  }
482 
483  // Return maximum double if at range limit
484  else if (std::abs(arg) == 1.0) {
485  result = copysign(1.0, arg) * DBL_MAX;
486  }
487 
488  // Otherwise compute inverse of error function
489  else {
490 
491  // Compute a rational approximation of the inverse error function
492  if (std::abs(arg) <= 0.7) {
493  double z = arg * arg;
494  double num = (((a[3]*z + a[2])*z + a[1])*z + a[0]);
495  double dem = ((((b[3]*z + b[2])*z + b[1])*z +b[0])*z + 1.0);
496  result = arg * num/dem;
497  }
498  else {
499  double z = std::sqrt(-std::log((1.0-std::abs(arg))/2.0));
500  double num = ((c[3]*z + c[2])*z + c[1])*z + c[0];
501  double dem = (d[1]*z + d[0])*z + 1.0;
502  result = (copysign(1.0, arg)) * num/dem;
503  }
504 
505  // Now do two steps of Newton-Raphson correction
506  result -= (erf(result) - arg) / (e * std::exp(-result*result));
507  result -= (erf(result) - arg) / (e * std::exp(-result*result));
508 
509  }
510 
511  // Return result
512  return result;
513 }
514 
515 
516 /***********************************************************************//**
517  * @brief Returns the remainder of the division
518  *
519  * @param[in] v1 Nominator.
520  * @param[in] v2 Denominator.
521  * @return Remainder of division
522  *
523  * Returns the remainder of the division @a v1/v2. The result is non-negative.
524  * @a v1 can be positive or negative; @a v2 must be positive.
525  ***************************************************************************/
526 double gammalib::modulo(const double& v1, const double& v2)
527 {
528  // Declare result
529  double result;
530 
531  //
532  if (v1 >= 0.0) {
533  result = (v1 < v2) ? v1 : std::fmod(v1,v2);
534  }
535  else {
536  double tmp = std::fmod(v1,v2) + v2;
537  result = (tmp == v2) ? 0.0 : tmp;
538  }
539 
540  // Return result
541  return result;
542 }
543 
544 
545 /***********************************************************************//**
546  * @brief Returns the integral of a power law
547  *
548  * @param[in] x1 First x value.
549  * @param[in] f1 Power law value at first x value.
550  * @param[in] x2 Second x value.
551  * @param[in] f2 Power law value at second x value.
552  * @return Integral of power law
553  *
554  * Analytically computes
555  *
556  * \f[\int_{x_1}^{x_2} F_1 \left( \frac{x}{x_1} \right)^m dx\f]
557  *
558  * where
559  *
560  * \f[m = \frac{\ln (F_2 / F_1)}{\ln (x_2 / x_1)}\f]
561  *
562  * and
563  * \f$F_1\f$ is the power law value at point \f$x_1\f$ and
564  * \f$F_2\f$ is the power law value at point \f$x_2\f$.
565  ***************************************************************************/
566 double gammalib::plaw_integral(const double& x1,
567  const double& f1,
568  const double& x2,
569  const double& f2)
570 {
571  // Compute power law slope
572  double x2x1 = x2/x1;
573  double fratio = std::log(f2/f1);
574  double xratio = std::log(x2x1);
575  double slope = fratio / xratio;
576  double slopep1 = slope + 1.0;
577 
578  // Compute integral. Computations dependend on the slope. We add here a
579  // kluge to assure numerical accuracy.
580  double integral;
581  if (std::abs(slopep1) > 1.0e-11) {
582  integral = f1 / slopep1 * (x2 * std::pow(x2x1, slope) - x1);
583  }
584  else {
585  integral = f1 * x1 * xratio;
586  }
587 
588  // Return integral
589  return integral;
590 }
591 
592 
593 /***********************************************************************//**
594  * @brief Returns the integral of a Gaussian function
595  *
596  * @param[in] x1 Lower x boundary (in units of Gaussian sigma).
597  * @param[in] x2 Upper x boundary (in units of Gaussian sigma).
598  * @return Integral of Gaussian
599  *
600  * Analytically computes
601  *
602  * \f[\frac{1}{\sqrt{\pi}}\int_{x_1}^{x_2} e^{-\frac{1}{2}x^2} dx\f]
603  *
604  ***************************************************************************/
605 double gammalib::gauss_integral(const double& x1,
606  const double& x2)
607 {
608  // Get integral
609  const double norm = 1.0 / gammalib::sqrt_two;
610  double xmin = x1 * norm;
611  double xmax = x2 * norm;
612  double integral = 0.5 * (gammalib::erfc(xmin) - gammalib::erfc(xmax));
613 
614  // Return integral
615  return integral;
616 }
double tand(const double &angle)
Compute tangens of angle in degrees.
Definition: GMath.cpp:192
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:864
const double pi
Definition: GMath.hpp:35
double erf(const double &arg)
Computes error function.
Definition: GMath.cpp:417
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1253
double gauss_integral(const double &x1, const double &x2)
Returns the integral of a Gaussian function.
Definition: GMath.cpp:605
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1190
double modulo(const double &v1, const double &v2)
Returns the remainder of the division.
Definition: GMath.cpp:526
double acos(const double &arg)
Computes acos by avoiding NaN due to rounding errors.
Definition: GMath.cpp:69
Gammalib tools definition.
const double sqrt_two
Definition: GMath.hpp:54
double atand(const double &value)
Compute arc tangens in degrees.
Definition: GMath.cpp:282
const double pihalf
Definition: GMath.hpp:38
const double deg2rad
Definition: GMath.hpp:43
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1358
#define copysign(X, Y)
Definition: GWcsMOL.cpp:44
double acosd(const double &value)
Compute arc cosine in degrees.
Definition: GMath.cpp:218
GVector tan(const GVector &vector)
Computes tangens of vector elements.
Definition: GVector.cpp:1379
double erfinv(const double &arg)
Computes inverse error function.
Definition: GMath.cpp:463
void sincosd(const double &angle, double *s, double *c)
Compute sine and cosine of angle in degrees.
Definition: GMath.cpp:342
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1274
double angle(const GVector &a, const GVector &b)
Computes angle between vectors.
Definition: GVector.cpp:974
double erfc(const double &arg)
Computes complementary error function.
Definition: GMath.cpp:450
GVector asin(const GVector &vector)
Computes arcsin of vector elements.
Definition: GVector.cpp:1106
double plaw_integral(const double &x1, const double &f1, const double &x2, const double &f2)
Returns the integral of a power law.
Definition: GMath.cpp:566
double gammln(const double &arg)
Computes logarithm of gamma function.
Definition: GMath.cpp:383
double atan2d(const double &y, const double &x)
Compute arc tangens in degrees.
Definition: GMath.cpp:308
double cosd(const double &angle)
Compute cosine of angle in degrees.
Definition: GMath.cpp:134
GVector pow(const GVector &vector, const double &power)
Computes tanh of vector elements.
Definition: GVector.cpp:1422
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1316
double asind(const double &value)
Compute arc sine in degrees.
Definition: GMath.cpp:250
double atan2(const double &y, const double &x)
Compute arc tangens in radians.
Definition: GMath.cpp:102
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1232
const double rad2deg
Definition: GMath.hpp:44
GVector atan(const GVector &vector)
Computes arctan of vector elements.
Definition: GVector.cpp:1148
double sind(const double &angle)
Compute sine of angle in degrees.
Definition: GMath.cpp:163
Mathematical function definitions.