-
Notifications
You must be signed in to change notification settings - Fork 0
/
paraboliccylinder.h
122 lines (113 loc) · 3.43 KB
/
paraboliccylinder.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
#ifndef PARABOLICCYLINDER_H
#define PARABOLICCYLINDER_H
#include "constants.h"
#include <boost/math/special_functions/bessel.hpp>
/*
Adopted from Randy Read's code by Pavel Afonine, 12-MAR-2014
*/
double dvsa(double, double);
double dvla(double, double);
double vvla(double, double);
double dv(double va, double x) {
/* Compute parabolic cylinder function Dv(x)
Equivalent to Mathematica ParabolicCylinderD[va,x]
Derived from routines in program mpbdv.for by Shanjie Zhang and Jianming
Jin distributed with their book "Computation of Special Functions",
Copyright 1996 by John Wiley & Sons, Inc.
Permission has been granted for purchasers of the book to incorporate these
programs as long as the copyright is acknowledged.
va = order of the parabolic cylinder function Dv
x = argument
*/
double ax(std::abs(x));
double pd;
if (ax <= 5.8)
pd = dvsa(va, x);
else
pd = dvla(va, x);
return pd;
}
double dvsa(double va, double x) {
// Compute parabolic cylinder function Dv(x) for small values of |x| (<=5.8)
static double EPS(std::pow(10., -15));
static double SQRT2(std::sqrt(2.));
static double SQRTPI(std::sqrt(scitbx::constants::pi));
double pd;
double ep = exp(-0.25 * x * x);
double va0 = 0.5 * (1. - va);
if (va == 0.)
pd = ep;
else {
if (x == 0.) {
double ftol = std::numeric_limits<double>::epsilon();
if (va0 <= 0. && std::abs(va0 - std::floor(va0 + 0.5)) < ftol)
pd = 0.;
else
pd = SQRTPI / (boost::math::tgamma(va0) * std::pow(2, -0.5 * va));
} else {
double a0 = std::pow(2, -0.5 * va - 1.) * ep / boost::math::tgamma(-va);
double vt = -0.5 * va;
pd = boost::math::tgamma(vt);
double r(1.), r1(pd);
int m(1);
while (m <= 250 && std::abs(r1) >= std::abs(pd) * EPS) {
double vm = 0.5 * (m - va);
r = -r * SQRT2 * x / m;
r1 = boost::math::tgamma(vm) * r;
pd += r1;
m++;
}
pd *= a0;
}
}
return pd;
}
double dvla(double va, double x) {
// Compute parabolic cylinder function Dv(x) for large values of |x| (>5.8)
static double EPS(std::pow(10., -12));
double xsqr(x * x);
double ep(exp(-0.25 * xsqr));
double a0(std::pow(std::abs(x), va) * ep);
double r(1.);
double pd(1.);
int k(1);
while (k <= 16 && std::abs(r / pd) >= EPS) {
r = -0.5 * r * (2. * k - va - 1.) * (2. * k - va - 2.) / (k * xsqr);
pd += r;
k++;
}
pd *= a0;
if (x < 0.) {
double x1(-x);
double vl = vvla(va, x1);
pd = scitbx::constants::pi * vl / boost::math::tgamma(-va) +
cos(scitbx::constants::pi * va) * pd;
}
return pd;
}
double vvla(double va, double x) {
// Compute parabolic cylinder function Vv(x) for large argument
static double EPS(std::pow(10., -12));
static double SQRT2BYPI(std::sqrt(2. / scitbx::constants::pi));
double xsqr(x * x);
double qe(exp(0.25 * xsqr));
double a0 = std::pow(std::abs(x), -va - 1.) * SQRT2BYPI * qe;
double r(1.), pv(1.);
int k(1);
while (k <= 18 && std::abs(r / pv) >= EPS) {
r = 0.5 * r * (2. * k + va - 1.) * (2. * k + va) / (k * xsqr);
pv += r;
k++;
}
pv *= a0;
if (x < 0.) {
double x1(-x);
double pdl = dvla(va, x1);
double gl = boost::math::tgamma(-va);
double piva = scitbx::constants::pi * va;
double dsl = std::pow(sin(piva), 2);
pv = dsl * gl / scitbx::constants::pi * pdl - cos(piva) * pv;
}
return pv;
}
#endif // PARABOLICCYLINDER_H