PLUTO
Functions
math_misc.c File Reference

Miscellaneous math functions. More...

#include "pluto.h"
Include dependency graph for math_misc.c:

Go to the source code of this file.

Functions

double BesselJ0 (double x)
 
double BesselJ1 (double x)
 
double BesselI0 (double x)
 
double BesselI1 (double x)
 
double BesselK0 (double x)
 
double BesselK1 (double x)
 
double BesselKn (int n, double x)
 
double RandomNumber (double rmin, double rmax)
 
void VectorCartesianComponents (double *v, double x1, double x2, double x3)
 

Detailed Description

Miscellaneous math functions.

Author
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t) B. Vaidya
Date
June 28, 2015

Definition in file math_misc.c.

Function Documentation

double BesselI0 ( double  x)

Returns the modified Bessel function I0(x) for any real x. (Adapted from Numerical Recipes)

Definition at line 80 of file math_misc.c.

86 {
87  double ax,ans;
88  double y; /* Accumulate polynomials in double precision.*/
89  if ((ax=fabs(x)) < 3.75) { /* Polynomial fit. */
90  y = x/3.75;
91  y *= y;
92  ans = 1.0 + y*(3.5156229+y*(3.0899424
93  + y*(1.2067492 + y*(0.2659732 + y*(0.360768e-1
94  + y*0.45813e-2)))));
95  } else {
96  y = 3.75/ax;
97  ans = (exp(ax)/sqrt(ax))*(0.39894228 + y*(0.1328592e-1
98  + y*(0.225319e-2 + y*(-0.157565e-2 + y*(0.916281e-2
99  + y*(-0.2057706e-1 + y*(0.2635537e-1
100  + y*(-0.1647633e-1 + y*0.392377e-2))))))));
101  }
102  return ans;
103 }

Here is the caller graph for this function:

double BesselI1 ( double  x)

Returns the modified Bessel function I1(x) for any real x. (Adapted from Numerical Recipes)

Definition at line 107 of file math_misc.c.

113 {
114  double ax, ans1, ans2, ans;
115  double y; /* Accumulate polynomials in double precision.*/
116  if ((ax=fabs(x)) < 3.75) { /* Polynomial fit. */
117  y = x/3.75;
118  y *= y;
119  ans = ax*(0.5+y*(0.87890594+y*(0.51498869+y*(0.15084934
120  +y*(0.2658733e-1+y*(0.301532e-2+y*0.32411e-3))))));
121  } else {
122  y = 3.75/ax;
123  ans1 = 0.2282967e-1+y*(-0.2895312e-1+y*(0.1787654e-1
124  -y*0.420059e-2));
125  ans2 = 0.39894228+y*(-0.3988024e-1 + y*(-0.362018e-2
126  + y*(0.163801e-2+y*(-0.1031555e-1+y*ans1))));
127  ans = ans2*(exp(ax)/sqrt(ax));
128  }
129  return x < 0.0 ? -ans : ans;
130 }

Here is the caller graph for this function:

double BesselJ0 ( double  x)

Returns the Bessel function J0(x) for any real x. (Adapted from Numerical Recipes)

Definition at line 13 of file math_misc.c.

19 {
20  double ax,z;
21  double xx,y,ans,ans1,ans2; /* Accumulate pol.in double prec. */
22 
23  if ((ax=fabs(x)) < 8.0) { /* Direct rational function fit. */
24  y = x*x;
25  ans1 = 57568490574.0+y*(-13362590354.0+y*(651619640.7
26  + y*(-11214424.18+y*(77392.33017+y*(-184.9052456)))));
27  ans2 = 57568490411.0+y*(1029532985.0+y*(9494680.718
28  + y*(59272.64853+y*(267.8532712+y*1.0))));
29  ans = ans1/ans2;
30  } else { /*Fitting function (6.5.9).*/
31  z = 8.0/ax;
32  y = z*z;
33  xx = ax-0.785398164;
34  ans1 = 1.0 + y*(-0.1098628627e-2+y*(0.2734510407e-4
35  + y*(-0.2073370639e-5+y*0.2093887211e-6)));
36  ans2 = - 0.1562499995e-1+y*(0.1430488765e-3
37  + y*(-0.6911147651e-5+y*(0.7621095161e-6
38  - y*0.934945152e-7)));
39  ans = sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
40  }
41  return ans;
42 }

Here is the caller graph for this function:

double BesselJ1 ( double  x)

Returns the Bessel function J1(x) for any real x. (Adapted from Numerical Recipes)

Definition at line 45 of file math_misc.c.

51 {
52  double ax,z;
53  double xx,y,ans,ans1,ans2; /* Accumulate polynomials in double precision.*/
54  if ((ax=fabs(x)) < 8.0) { /* Direct rational approximation.*/
55 
56  y=x*x;
57  ans1 = x*(72362614232.0+y*(-7895059235.0+y*(242396853.1
58  + y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
59  ans2 = 144725228442.0+y*(2300535178.0+y*(18583304.74
60  + y*(99447.43394+y*(376.9991397+y*1.0))));
61  ans = ans1/ans2;
62  } else { /*Fitting function (6.5.9).*/
63  z = 8.0/ax;
64  y = z*z;
65  xx = ax-2.356194491;
66  ans1 = 1.0 +y*(0.183105e-2+y*(-0.3516396496e-4
67  + y*(0.2457520174e-5+y*(-0.240337019e-6))));
68  ans2 = 0.04687499995+y*(-0.2002690873e-3
69  + y*(0.8449199096e-5+y*(-0.88228987e-6
70  + y*0.105787412e-6)));
71  ans = sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
72  if (x < 0.0) ans = -ans;
73  }
74  return ans;
75 }

Here is the caller graph for this function:

double BesselK0 ( double  x)

Returns the modified Bessel function K0(x) for any real x. (Adapted from Numerical Recipes)

Definition at line 134 of file math_misc.c.

140 {
141  double I0, ans;
142  double y; /* Accumulate polynomials in double precision.*/
143  if (x <= 2.0) { /* Polynomial fit. */
144  y = x*x/4.0;
145  I0 = BesselI0(x);
146  ans = (-log(x/2.0)*I0)+(-0.57721566+y*(0.42278420
147  + y*(0.23069756+y*(0.3488590e-1
148  + y*(0.262698e-2 +y*(0.10750e-3+y*0.74e-5))))));
149  } else {
150  y = 2.0/x;
151  ans = (exp(-x)/sqrt(x))*(1.25331414+y*(-0.7832358e-1
152  + y*(0.2189568e-1+y*(-0.1062446e-1+y*(0.587872e-2
153  + y*(-0.251540e-2+y*0.53208e-3))))));
154  }
155  return ans;
156 }
BesselI0
double BesselI0(double x)
Definition: math_misc.c:80

Here is the call graph for this function:

Here is the caller graph for this function:

double BesselK1 ( double  x)

Returns the modified Bessel function K1(x) for any real x. (Adapted from Numerical Recipes)

Definition at line 159 of file math_misc.c.

165 {
166  double I1, ans;
167  double y; /* Accumulate polynomials in double precision.*/
168  if (x <= 2.0) { /* Polynomial fit. */
169  y = x*x/4.0;
170  I1 = BesselI1(x);
171  ans = (log(x/2.0)*I1)+(1.0/x)*(1.0+y*(0.15443144
172  + y*(-0.67278579+y*(-0.18156897+y*(-0.1919402e-1
173  + y*(-0.110404e-2+y*(-0.4686e-4)))))));
174  }else{
175  y = 2.0/x;
176  ans = (exp(-x)/sqrt(x))*(1.25331414+y*(0.23498619 +y*(-0.3655620e-1
177  + y*(0.1504268e-1+y*(-0.780353e-2
178  + y*(0.325614e-2+y*(-0.68245e-3)))))));
179  }
180  return ans;
181 }
BesselI1
double BesselI1(double x)
Definition: math_misc.c:107

Here is the call graph for this function:

Here is the caller graph for this function:

double BesselKn ( int  n,
double  x 
)

Returns the modified Bessel function Kn(x) for positive x and n >= 2. (Adapted from Numerical Recipes)

Definition at line 184 of file math_misc.c.

190 {
191  int j;
192  double bk,bkm,bkp,tox;
193  if (n < 2){
194  print1("!BesselKn: Index n less than 2");
195  QUIT_PLUTO(1);
196  }
197  tox = 2.0/x;
198  bkm = BesselK0(x);
199  bk = BesselK1(x);
200  for (j=1;j<n;j++) {
201  bkp = bkm + j*tox*bk;
202  bkm = bk;
203  bk = bkp;
204  }
205  return bk;
206 }
BesselK1
double BesselK1(double x)
Definition: math_misc.c:159
n
static int n
Definition: analysis.c:3
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
BesselK0
double BesselK0(double x)
Definition: math_misc.c:134
j
int j
Definition: analysis.c:2
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125

Here is the call graph for this function:

double RandomNumber ( double  rmin,
double  rmax 
)

Generate and return a random number between [rmin, rmax]

Definition at line 209 of file math_misc.c.

213 {
214  static int first_call = 1;
215  double eps, rnd;
216 
217  if (first_call == 1){
218  srand(time(NULL) + prank);
219  first_call = 0;
220  }
221  eps = (double)(rand())/((double)RAND_MAX + 1.0); /* 0 < eps < 1 */
222  rnd = rmin + eps*(rmax-rmin); /* Map random number between [rmin,rmax] */
223  return rnd;
224 }
prank
int prank
Processor rank.
Definition: globals.h:33
eps
#define eps

Here is the caller graph for this function:

void VectorCartesianComponents ( double *  v,
double  x1,
double  x2,
double  x3 
)

Transform the vector $ v = (v_{x1}, v_{x2}, v_{x3})$ from the chosen coordinate system (= GEOMETRY) to Cartesian components.

Parameters
[in,out]*vthe original vector. On output v is replaced by the three Cartesian components.
[in]x1,x2,x3the grid coordinates

Definition at line 227 of file math_misc.c.

238 {
239  double vx1, vx2, vx3;
240 
241  vx1 = v[0];
242  vx2 = v[1];
243  vx3 = v[2];
244 
245  #if GEOMETRY == CARTESIAN
246 
247  /* nothing to do */
248 
249  #elif GEOMETRY == POLAR
250 
251  EXPAND(v[0] = vx1*cos(x2) - vx2*sin(x2); ,
252  v[1] = vx1*sin(x2) + vx2*cos(x2); ,
253  v[2] = vx3;)
254 
255  #elif GEOMETRY == SPHERICAL
256 
257  #if DIMENSIONS == 2
258  EXPAND(v[0] = vx1*sin(x2) + vx2*cos(x2); ,
259  v[1] = vx1*cos(x2) - vx2*sin(x2); ,
260  v[2] = vx3;)
261  #elif DIMENSIONS == 3
262  v[0] = (vx1*sin(x2) + vx2*cos(x2))*cos(x3) - vx3*sin(x3);
263  v[1] = (vx1*sin(x2) + vx2*cos(x2))*sin(x3) + vx3*cos(x3);
264  v[2] = (vx1*cos(x2) - vx2*sin(x2));
265  #endif
266  #endif
267 }
func_param::v
double v[NVAR]
Definition: eos.h:106
SPHERICAL
#define SPHERICAL
Definition: pluto.h:36
GEOMETRY
#define GEOMETRY
Definition: definitions_01.h:4
DIMENSIONS
#define DIMENSIONS
Definition: definitions_01.h:2

Here is the caller graph for this function: