[706f466] | 1 | r""" |
---|
[597878a] | 2 | Special Functions |
---|
| 3 | ................. |
---|
| 4 | |
---|
[d0dc9a3] | 5 | This following standard C99 math functions are available: |
---|
[597878a] | 6 | |
---|
| 7 | M_PI, M_PI_2, M_PI_4, M_SQRT1_2, M_E: |
---|
| 8 | $\pi$, $\pi/2$, $\pi/4$, $1/\sqrt{2}$ and Euler's constant $e$ |
---|
[706f466] | 9 | |
---|
[d0dc9a3] | 10 | exp, log, pow(x,y), expm1, log1p, sqrt, cbrt: |
---|
| 11 | Power functions $e^x$, $\ln x$, $x^y$, $e^x - 1$, $\ln 1 + x$, |
---|
| 12 | $\sqrt{x}$, $\sqrt[3]{x}$. The functions expm1(x) and log1p(x) |
---|
| 13 | are accurate across all $x$, including $x$ very close to zero. |
---|
[706f466] | 14 | |
---|
[597878a] | 15 | sin, cos, tan, asin, acos, atan: |
---|
| 16 | Trigonometry functions and inverses, operating on radians. |
---|
[706f466] | 17 | |
---|
[597878a] | 18 | sinh, cosh, tanh, asinh, acosh, atanh: |
---|
| 19 | Hyperbolic trigonometry functions. |
---|
[706f466] | 20 | |
---|
[597878a] | 21 | atan2(y,x): |
---|
| 22 | Angle from the $x$\ -axis to the point $(x,y)$, which is equal to |
---|
| 23 | $\tan^{-1}(y/x)$ corrected for quadrant. That is, if $x$ and $y$ are |
---|
| 24 | both negative, then atan2(y,x) returns a value in quadrant III where |
---|
| 25 | atan(y/x) would return a value in quadrant I. Similarly for |
---|
| 26 | quadrants II and IV when $x$ and $y$ have opposite sign. |
---|
[706f466] | 27 | |
---|
[d0dc9a3] | 28 | fabs(x), fmin(x,y), fmax(x,y), trunc, rint: |
---|
[597878a] | 29 | Floating point functions. rint(x) returns the nearest integer. |
---|
[706f466] | 30 | |
---|
[597878a] | 31 | NAN: |
---|
| 32 | NaN, Not a Number, $0/0$. Use isnan(x) to test for NaN. Note that |
---|
| 33 | you cannot use :code:`x == NAN` to test for NaN values since that |
---|
[d0dc9a3] | 34 | will always return false. NAN does not equal NAN! The alternative, |
---|
| 35 | :code:`x != x` may fail if the compiler optimizes the test away. |
---|
[706f466] | 36 | |
---|
[597878a] | 37 | INFINITY: |
---|
| 38 | $\infty, 1/0$. Use isinf(x) to test for infinity, or isfinite(x) |
---|
| 39 | to test for finite and not NaN. |
---|
[706f466] | 40 | |
---|
[597878a] | 41 | erf, erfc, tgamma, lgamma: **do not use** |
---|
| 42 | Special functions that should be part of the standard, but are missing |
---|
| 43 | or inaccurate on some platforms. Use sas_erf, sas_erfc and sas_gamma |
---|
| 44 | instead (see below). Note: lgamma(x) has not yet been tested. |
---|
| 45 | |
---|
| 46 | Some non-standard constants and functions are also provided: |
---|
| 47 | |
---|
| 48 | M_PI_180, M_4PI_3: |
---|
| 49 | $\frac{\pi}{180}$, $\frac{4\pi}{3}$ |
---|
[706f466] | 50 | |
---|
[597878a] | 51 | SINCOS(x, s, c): |
---|
| 52 | Macro which sets s=sin(x) and c=cos(x). The variables *c* and *s* |
---|
| 53 | must be declared first. |
---|
[706f466] | 54 | |
---|
[597878a] | 55 | square(x): |
---|
| 56 | $x^2$ |
---|
[706f466] | 57 | |
---|
[597878a] | 58 | cube(x): |
---|
| 59 | $x^3$ |
---|
[706f466] | 60 | |
---|
[597878a] | 61 | sas_sinx_x(x): |
---|
| 62 | $\sin(x)/x$, with limit $\sin(0)/0 = 1$. |
---|
[706f466] | 63 | |
---|
[597878a] | 64 | powr(x, y): |
---|
| 65 | $x^y$ for $x \ge 0$; this is faster than general $x^y$ on some GPUs. |
---|
[706f466] | 66 | |
---|
[597878a] | 67 | pown(x, n): |
---|
| 68 | $x^n$ for $n$ integer; this is faster than general $x^n$ on some GPUs. |
---|
[706f466] | 69 | |
---|
[597878a] | 70 | FLOAT_SIZE: |
---|
| 71 | The number of bytes in a floating point value. Even though all |
---|
| 72 | variables are declared double, they may be converted to single |
---|
| 73 | precision float before running. If your algorithm depends on |
---|
| 74 | precision (which is not uncommon for numerical algorithms), use |
---|
| 75 | the following:: |
---|
| 76 | |
---|
| 77 | #if FLOAT_SIZE>4 |
---|
| 78 | ... code for double precision ... |
---|
| 79 | #else |
---|
| 80 | ... code for single precision ... |
---|
| 81 | #endif |
---|
[706f466] | 82 | |
---|
[597878a] | 83 | SAS_DOUBLE: |
---|
| 84 | A replacement for :code:`double` so that the declared variable will |
---|
| 85 | stay double precision; this should generally not be used since some |
---|
| 86 | graphics cards do not support double precision. There is no provision |
---|
| 87 | for forcing a constant to stay double precision. |
---|
| 88 | |
---|
[d0dc9a3] | 89 | The following special functions and scattering calculations are defined. |
---|
[597878a] | 90 | These functions have been tuned to be fast and numerically stable down |
---|
| 91 | to $q=0$ even in single precision. In some cases they work around bugs |
---|
| 92 | which appear on some platforms but not others, so use them where needed. |
---|
| 93 | Add the files listed in :code:`source = ["lib/file.c", ...]` to your *model.py* |
---|
| 94 | file in the order given, otherwise these functions will not be available. |
---|
| 95 | |
---|
| 96 | polevl(x, c, n): |
---|
| 97 | Polynomial evaluation $p(x) = \sum_{i=0}^n c_i x^i$ using Horner's |
---|
| 98 | method so it is faster and more accurate. |
---|
| 99 | |
---|
| 100 | $c = \{c_n, c_{n-1}, \ldots, c_0 \}$ is the table of coefficients, |
---|
| 101 | sorted from highest to lowest. |
---|
| 102 | |
---|
| 103 | p1evl(x, c, n): |
---|
[110f69c] | 104 | Evaluate normalized polynomial $p(x) = x^n + \sum_{i=0}^{n-1} c_i x^i$ |
---|
[597878a] | 105 | using Horner's method so it is faster and more accurate. |
---|
| 106 | |
---|
| 107 | $c = \{c_{n-1}, c_{n-2} \ldots, c_0 \}$ is the table of coefficients, |
---|
| 108 | sorted from highest to lowest. |
---|
| 109 | |
---|
| 110 | sas_gamma(x): |
---|
| 111 | Gamma function $\text{sas_gamma}(x) = \Gamma(x)$. |
---|
| 112 | |
---|
| 113 | The standard math function, tgamma(x) is unstable for $x < 1$ |
---|
| 114 | on some platforms. |
---|
| 115 | |
---|
[fba9ca0] | 116 | sas_gammaln(x): |
---|
| 117 | log gamma function sas_gammaln\ $(x) = \log \Gamma(|x|)$. |
---|
| 118 | |
---|
| 119 | The standard math function, lgamma(x), is incorrect for single |
---|
| 120 | precision on some platforms. |
---|
| 121 | |
---|
| 122 | sas_gammainc(a, x), sas_gammaincc(a, x): |
---|
| 123 | Incomplete gamma function |
---|
| 124 | sas_gammainc\ $(a, x) = \int_0^x t^{a-1}e^{-t}\,dt / \Gamma(a)$ |
---|
| 125 | and complementary incomplete gamma function |
---|
| 126 | sas_gammaincc\ $(a, x) = \int_x^\infty t^{a-1}e^{-t}\,dt / \Gamma(a)$ |
---|
| 127 | |
---|
[597878a] | 128 | sas_erf(x), sas_erfc(x): |
---|
| 129 | Error function |
---|
| 130 | $\text{sas_erf}(x) = \frac{2}{\sqrt\pi}\int_0^x e^{-t^2}\,dt$ |
---|
| 131 | and complementary error function |
---|
| 132 | $\text{sas_erfc}(x) = \frac{2}{\sqrt\pi}\int_x^{\infty} e^{-t^2}\,dt$. |
---|
| 133 | |
---|
| 134 | The standard math functions erf(x) and erfc(x) are slower and broken |
---|
| 135 | on some platforms. |
---|
| 136 | |
---|
| 137 | sas_J0(x): |
---|
| 138 | Bessel function of the first kind $\text{sas_J0}(x)=J_0(x)$ where |
---|
| 139 | $J_0(x) = \frac{1}{\pi}\int_0^\pi \cos(x\sin(\tau))\,d\tau$. |
---|
| 140 | |
---|
| 141 | The standard math function j0(x) is not available on all platforms. |
---|
| 142 | |
---|
| 143 | sas_J1(x): |
---|
| 144 | Bessel function of the first kind $\text{sas_J1}(x)=J_1(x)$ where |
---|
| 145 | $J_1(x) = \frac{1}{\pi}\int_0^\pi \cos(\tau - x\sin(\tau))\,d\tau$. |
---|
| 146 | |
---|
| 147 | The standard math function j1(x) is not available on all platforms. |
---|
| 148 | |
---|
| 149 | sas_JN(n, x): |
---|
| 150 | Bessel function of the first kind and integer order $n$: |
---|
| 151 | $\text{sas_JN}(n, x)=J_n(x)$ where |
---|
| 152 | $J_n(x) = \frac{1}{\pi}\int_0^\pi \cos(n\tau - x\sin(\tau))\,d\tau$. |
---|
| 153 | If $n$ = 0 or 1, it uses sas_J0(x) or sas_J1(x), respectively. |
---|
| 154 | |
---|
| 155 | The standard math function jn(n, x) is not available on all platforms. |
---|
| 156 | |
---|
| 157 | sas_Si(x): |
---|
| 158 | Sine integral $\text{Si}(x) = \int_0^x \tfrac{\sin t}{t}\,dt$. |
---|
| 159 | |
---|
| 160 | This function uses Taylor series for small and large arguments: |
---|
| 161 | |
---|
| 162 | For large arguments, |
---|
| 163 | |
---|
| 164 | .. math:: |
---|
| 165 | |
---|
[110f69c] | 166 | \text{Si}(x) \sim \frac{\pi}{2} |
---|
| 167 | - \frac{\cos(x)}{x} |
---|
| 168 | \left(1 - \frac{2!}{x^2} + \frac{4!}{x^4} - \frac{6!}{x^6} \right) |
---|
| 169 | - \frac{\sin(x)}{x} |
---|
| 170 | \left(\frac{1}{x} - \frac{3!}{x^3} + \frac{5!}{x^5} - \frac{7!}{x^7}\right) |
---|
[597878a] | 171 | |
---|
| 172 | For small arguments, |
---|
| 173 | |
---|
| 174 | .. math:: |
---|
| 175 | |
---|
| 176 | \text{Si}(x) \sim x |
---|
| 177 | - \frac{x^3}{3\times 3!} + \frac{x^5}{5 \times 5!} - \frac{x^7}{7 \times 7!} |
---|
| 178 | + \frac{x^9}{9\times 9!} - \frac{x^{11}}{11\times 11!} |
---|
| 179 | |
---|
| 180 | sas_3j1x_x(x): |
---|
| 181 | Spherical Bessel form |
---|
| 182 | $\text{sph_j1c}(x) = 3 j_1(x)/x = 3 (\sin(x) - x \cos(x))/x^3$, |
---|
| 183 | with a limiting value of 1 at $x=0$, where $j_1(x)$ is the spherical |
---|
| 184 | Bessel function of the first kind and first order. |
---|
| 185 | |
---|
| 186 | This function uses a Taylor series for small $x$ for numerical accuracy. |
---|
| 187 | |
---|
| 188 | |
---|
| 189 | sas_2J1x_x(x): |
---|
| 190 | Bessel form $\text{sas_J1c}(x) = 2 J_1(x)/x$, with a limiting value |
---|
| 191 | of 1 at $x=0$, where $J_1(x)$ is the Bessel function of first kind |
---|
| 192 | and first order. |
---|
| 193 | |
---|
| 194 | |
---|
[d0dc9a3] | 195 | gauss76.n, gauss76.z[i], gauss76.w[i]: |
---|
[597878a] | 196 | Points $z_i$ and weights $w_i$ for 76-point Gaussian quadrature, respectively, |
---|
| 197 | computing $\int_{-1}^1 f(z)\,dz \approx \sum_{i=1}^{76} w_i\,f(z_i)$. |
---|
[d0dc9a3] | 198 | When translating the model to C, include 'lib/gauss76.c' in the source |
---|
| 199 | and use :code:`GAUSS_N`, :code:`GAUSS_Z`, and :code:`GAUSS_W`. |
---|
[597878a] | 200 | |
---|
[d0dc9a3] | 201 | Similar arrays are available in :code:`gauss20` for 20-point quadrature |
---|
| 202 | and :code:`gauss150.c` for 150-point quadrature. By using |
---|
| 203 | :code:`import gauss76 as gauss` it is easy to change the number of |
---|
| 204 | points in the integration. |
---|
[597878a] | 205 | """ |
---|
[e65c3ba] | 206 | # pylint: disable=unused-import |
---|
| 207 | |
---|
[597878a] | 208 | import numpy as np |
---|
| 209 | |
---|
| 210 | # Functions to add to our standard set |
---|
| 211 | from numpy import degrees, radians |
---|
| 212 | |
---|
| 213 | # C99 standard math library functions |
---|
[df69efa] | 214 | from numpy import exp, log, power as pow, expm1, log1p, sqrt, cbrt |
---|
[597878a] | 215 | from numpy import sin, cos, tan, arcsin as asin, arccos as acos, arctan as atan |
---|
| 216 | from numpy import sinh, cosh, tanh, arcsinh as asinh, arccosh as acosh, arctanh as atanh |
---|
| 217 | from numpy import arctan2 as atan2 |
---|
[d0dc9a3] | 218 | from numpy import fabs, fmin, fmax, trunc, rint |
---|
| 219 | from numpy import pi, nan, inf |
---|
[e65c3ba] | 220 | from scipy.special import gamma as sas_gamma |
---|
[fba9ca0] | 221 | from scipy.special import gammaln as sas_gammaln |
---|
| 222 | from scipy.special import gammainc as sas_gammainc |
---|
| 223 | from scipy.special import gammaincc as sas_gammaincc |
---|
[e65c3ba] | 224 | from scipy.special import erf as sas_erf |
---|
| 225 | from scipy.special import erfc as sas_erfc |
---|
| 226 | from scipy.special import j0 as sas_J0 |
---|
| 227 | from scipy.special import j1 as sas_J1 |
---|
| 228 | from scipy.special import jn as sas_JN |
---|
| 229 | |
---|
[597878a] | 230 | # erf, erfc, tgamma, lgamma **do not use** |
---|
| 231 | |
---|
[e65c3ba] | 232 | # C99 standard math constants |
---|
| 233 | M_PI, M_PI_2, M_PI_4, M_SQRT1_2, M_E = np.pi, np.pi/2, np.pi/4, np.sqrt(0.5), np.e |
---|
[d0dc9a3] | 234 | NAN = nan |
---|
| 235 | INFINITY = inf |
---|
[e65c3ba] | 236 | |
---|
| 237 | # non-standard constants |
---|
[597878a] | 238 | M_PI_180, M_4PI_3 = M_PI/180, 4*M_PI/3 |
---|
| 239 | |
---|
| 240 | # can't do SINCOS in python; use "s, c = SINCOS(x)" instead |
---|
[e65c3ba] | 241 | def SINCOS(x): |
---|
| 242 | """return sin(x), cos(x)""" |
---|
| 243 | return sin(x), cos(x) |
---|
[d0dc9a3] | 244 | sincos = SINCOS |
---|
[e65c3ba] | 245 | |
---|
| 246 | def square(x): |
---|
| 247 | """return x^2""" |
---|
| 248 | return x*x |
---|
| 249 | |
---|
| 250 | def cube(x): |
---|
| 251 | """return x^3""" |
---|
| 252 | return x*x*x |
---|
| 253 | |
---|
| 254 | def sas_sinx_x(x): |
---|
| 255 | """return sin(x)/x""" |
---|
| 256 | from numpy import sinc as _sinc |
---|
| 257 | return _sinc(x/M_PI) |
---|
| 258 | |
---|
| 259 | def powr(x, y): |
---|
| 260 | """return x^y for x>0""" |
---|
| 261 | return x**y |
---|
| 262 | def pown(x, n): |
---|
| 263 | """return x^n for n integer""" |
---|
| 264 | return x**n |
---|
[597878a] | 265 | |
---|
| 266 | FLOAT_SIZE = 8 |
---|
| 267 | |
---|
[e65c3ba] | 268 | def polevl(x, c, n): |
---|
| 269 | """return p(x) for polynomial p of degree n-1 with coefficients c""" |
---|
| 270 | return np.polyval(c[:n], x) |
---|
[597878a] | 271 | |
---|
[e65c3ba] | 272 | def p1evl(x, c, n): |
---|
| 273 | """return x^n + p(x) for polynomial p of degree n-1 with coefficients c""" |
---|
| 274 | return np.polyval(np.hstack(([1.], c))[:n], x) |
---|
[597878a] | 275 | |
---|
| 276 | def sas_Si(x): |
---|
[e65c3ba] | 277 | """return Si(x)""" |
---|
| 278 | from scipy.special import sici |
---|
| 279 | return sici(x)[0] |
---|
[597878a] | 280 | |
---|
[4f611f1] | 281 | def sas_j1(x): |
---|
[e65c3ba] | 282 | """return j1(x)""" |
---|
[4f611f1] | 283 | if np.isscalar(x): |
---|
| 284 | retvalue = (sin(x) - x*cos(x))/x**2 if x != 0. else 0. |
---|
| 285 | else: |
---|
| 286 | with np.errstate(all='ignore'): |
---|
| 287 | retvalue = (sin(x) - x*cos(x))/x**2 |
---|
| 288 | retvalue[x == 0.] = 0. |
---|
[e65c3ba] | 289 | return retvalue |
---|
[706f466] | 290 | |
---|
[597878a] | 291 | def sas_3j1x_x(x): |
---|
[e65c3ba] | 292 | """return 3*j1(x)/x""" |
---|
[4f611f1] | 293 | if np.isscalar(x): |
---|
| 294 | retvalue = 3*(sin(x) - x*cos(x))/x**3 if x != 0. else 1. |
---|
| 295 | else: |
---|
| 296 | with np.errstate(all='ignore'): |
---|
| 297 | retvalue = 3*(sin(x) - x*cos(x))/x**3 |
---|
| 298 | retvalue[x == 0.] = 1. |
---|
| 299 | return retvalue |
---|
[859567e] | 300 | |
---|
[597878a] | 301 | def sas_2J1x_x(x): |
---|
[e65c3ba] | 302 | """return 2*J1(x)/x""" |
---|
[4f611f1] | 303 | if np.isscalar(x): |
---|
| 304 | retvalue = 2*sas_J1(x)/x if x != 0 else 1. |
---|
| 305 | else: |
---|
| 306 | with np.errstate(all='ignore'): |
---|
| 307 | retvalue = 2*sas_J1(x)/x |
---|
[110f69c] | 308 | retvalue[x == 0] = 1. |
---|
[4f611f1] | 309 | return retvalue |
---|
[597878a] | 310 | |
---|
| 311 | |
---|
| 312 | # Gaussians |
---|
[b297ba9] | 313 | class Gauss(object): |
---|
| 314 | """Gauss-Legendre integration weights""" |
---|
[d0dc9a3] | 315 | def __init__(self, w, z): |
---|
| 316 | self.n = len(w) |
---|
| 317 | self.w = w |
---|
| 318 | self.z = z |
---|
| 319 | |
---|
| 320 | gauss20 = Gauss( |
---|
| 321 | w=np.array([ |
---|
| 322 | .0176140071391521, |
---|
| 323 | .0406014298003869, |
---|
| 324 | .0626720483341091, |
---|
| 325 | .0832767415767047, |
---|
| 326 | .10193011981724, |
---|
| 327 | .118194531961518, |
---|
| 328 | .131688638449177, |
---|
| 329 | .142096109318382, |
---|
| 330 | .149172986472604, |
---|
| 331 | .152753387130726, |
---|
| 332 | .152753387130726, |
---|
| 333 | .149172986472604, |
---|
| 334 | .142096109318382, |
---|
| 335 | .131688638449177, |
---|
| 336 | .118194531961518, |
---|
| 337 | .10193011981724, |
---|
| 338 | .0832767415767047, |
---|
| 339 | .0626720483341091, |
---|
| 340 | .0406014298003869, |
---|
| 341 | .0176140071391521 |
---|
| 342 | ]), |
---|
| 343 | z=np.array([ |
---|
| 344 | -.993128599185095, |
---|
| 345 | -.963971927277914, |
---|
| 346 | -.912234428251326, |
---|
| 347 | -.839116971822219, |
---|
| 348 | -.746331906460151, |
---|
| 349 | -.636053680726515, |
---|
| 350 | -.510867001950827, |
---|
| 351 | -.37370608871542, |
---|
| 352 | -.227785851141645, |
---|
| 353 | -.076526521133497, |
---|
| 354 | .0765265211334973, |
---|
| 355 | .227785851141645, |
---|
| 356 | .37370608871542, |
---|
| 357 | .510867001950827, |
---|
| 358 | .636053680726515, |
---|
| 359 | .746331906460151, |
---|
| 360 | .839116971822219, |
---|
| 361 | .912234428251326, |
---|
| 362 | .963971927277914, |
---|
| 363 | .993128599185095 |
---|
| 364 | ]) |
---|
| 365 | ) |
---|
| 366 | |
---|
| 367 | gauss76 = Gauss( |
---|
| 368 | w=np.array([ |
---|
| 369 | .00126779163408536, #0 |
---|
| 370 | .00294910295364247, |
---|
| 371 | .00462793522803742, |
---|
| 372 | .00629918049732845, |
---|
| 373 | .00795984747723973, |
---|
| 374 | .00960710541471375, |
---|
| 375 | .0112381685696677, |
---|
| 376 | .0128502838475101, |
---|
| 377 | .0144407317482767, |
---|
| 378 | .0160068299122486, |
---|
| 379 | .0175459372914742, #10 |
---|
| 380 | .0190554584671906, |
---|
| 381 | .020532847967908, |
---|
| 382 | .0219756145344162, |
---|
| 383 | .0233813253070112, |
---|
| 384 | .0247476099206597, |
---|
| 385 | .026072164497986, |
---|
| 386 | .0273527555318275, |
---|
| 387 | .028587223650054, |
---|
| 388 | .029773487255905, |
---|
| 389 | .0309095460374916, #20 |
---|
| 390 | .0319934843404216, |
---|
| 391 | .0330234743977917, |
---|
| 392 | .0339977794120564, |
---|
| 393 | .0349147564835508, |
---|
| 394 | .0357728593807139, |
---|
| 395 | .0365706411473296, |
---|
| 396 | .0373067565423816, |
---|
| 397 | .0379799643084053, |
---|
| 398 | .0385891292645067, |
---|
| 399 | .0391332242205184, #30 |
---|
| 400 | .0396113317090621, |
---|
| 401 | .0400226455325968, |
---|
| 402 | .040366472122844, |
---|
| 403 | .0406422317102947, |
---|
| 404 | .0408494593018285, |
---|
| 405 | .040987805464794, |
---|
| 406 | .0410570369162294, |
---|
| 407 | .0410570369162294, |
---|
| 408 | .040987805464794, |
---|
| 409 | .0408494593018285, #40 |
---|
| 410 | .0406422317102947, |
---|
| 411 | .040366472122844, |
---|
| 412 | .0400226455325968, |
---|
| 413 | .0396113317090621, |
---|
| 414 | .0391332242205184, |
---|
| 415 | .0385891292645067, |
---|
| 416 | .0379799643084053, |
---|
| 417 | .0373067565423816, |
---|
| 418 | .0365706411473296, |
---|
| 419 | .0357728593807139, #50 |
---|
| 420 | .0349147564835508, |
---|
| 421 | .0339977794120564, |
---|
| 422 | .0330234743977917, |
---|
| 423 | .0319934843404216, |
---|
| 424 | .0309095460374916, |
---|
| 425 | .029773487255905, |
---|
| 426 | .028587223650054, |
---|
| 427 | .0273527555318275, |
---|
| 428 | .026072164497986, |
---|
| 429 | .0247476099206597, #60 |
---|
| 430 | .0233813253070112, |
---|
| 431 | .0219756145344162, |
---|
| 432 | .020532847967908, |
---|
| 433 | .0190554584671906, |
---|
| 434 | .0175459372914742, |
---|
| 435 | .0160068299122486, |
---|
| 436 | .0144407317482767, |
---|
| 437 | .0128502838475101, |
---|
| 438 | .0112381685696677, |
---|
| 439 | .00960710541471375, #70 |
---|
| 440 | .00795984747723973, |
---|
| 441 | .00629918049732845, |
---|
| 442 | .00462793522803742, |
---|
| 443 | .00294910295364247, |
---|
| 444 | .00126779163408536 #75 (indexed from 0) |
---|
| 445 | ]), |
---|
| 446 | z=np.array([ |
---|
| 447 | -.999505948362153, #0 |
---|
| 448 | -.997397786355355, |
---|
| 449 | -.993608772723527, |
---|
| 450 | -.988144453359837, |
---|
| 451 | -.981013938975656, |
---|
| 452 | -.972229228520377, |
---|
| 453 | -.961805126758768, |
---|
| 454 | -.949759207710896, |
---|
| 455 | -.936111781934811, |
---|
| 456 | -.92088586125215, |
---|
| 457 | -.904107119545567, #10 |
---|
| 458 | -.885803849292083, |
---|
| 459 | -.866006913771982, |
---|
| 460 | -.844749694983342, |
---|
| 461 | -.822068037328975, |
---|
| 462 | -.7980001871612, |
---|
| 463 | -.77258672828181, |
---|
| 464 | -.74587051350361, |
---|
| 465 | -.717896592387704, |
---|
| 466 | -.688712135277641, |
---|
| 467 | -.658366353758143, #20 |
---|
| 468 | -.626910417672267, |
---|
| 469 | -.594397368836793, |
---|
| 470 | -.560882031601237, |
---|
| 471 | -.526420920401243, |
---|
| 472 | -.491072144462194, |
---|
| 473 | -.454895309813726, |
---|
| 474 | -.417951418780327, |
---|
| 475 | -.380302767117504, |
---|
| 476 | -.342012838966962, |
---|
| 477 | -.303146199807908, #30 |
---|
| 478 | -.263768387584994, |
---|
| 479 | -.223945802196474, |
---|
| 480 | -.183745593528914, |
---|
| 481 | -.143235548227268, |
---|
| 482 | -.102483975391227, |
---|
| 483 | -.0615595913906112, |
---|
| 484 | -.0205314039939986, |
---|
| 485 | .0205314039939986, |
---|
| 486 | .0615595913906112, |
---|
| 487 | .102483975391227, #40 |
---|
| 488 | .143235548227268, |
---|
| 489 | .183745593528914, |
---|
| 490 | .223945802196474, |
---|
| 491 | .263768387584994, |
---|
| 492 | .303146199807908, |
---|
| 493 | .342012838966962, |
---|
| 494 | .380302767117504, |
---|
| 495 | .417951418780327, |
---|
| 496 | .454895309813726, |
---|
| 497 | .491072144462194, #50 |
---|
| 498 | .526420920401243, |
---|
| 499 | .560882031601237, |
---|
| 500 | .594397368836793, |
---|
| 501 | .626910417672267, |
---|
| 502 | .658366353758143, |
---|
| 503 | .688712135277641, |
---|
| 504 | .717896592387704, |
---|
| 505 | .74587051350361, |
---|
| 506 | .77258672828181, |
---|
| 507 | .7980001871612, #60 |
---|
| 508 | .822068037328975, |
---|
| 509 | .844749694983342, |
---|
| 510 | .866006913771982, |
---|
| 511 | .885803849292083, |
---|
| 512 | .904107119545567, |
---|
| 513 | .92088586125215, |
---|
| 514 | .936111781934811, |
---|
| 515 | .949759207710896, |
---|
| 516 | .961805126758768, |
---|
| 517 | .972229228520377, #70 |
---|
| 518 | .981013938975656, |
---|
| 519 | .988144453359837, |
---|
| 520 | .993608772723527, |
---|
| 521 | .997397786355355, |
---|
| 522 | .999505948362153 #75 |
---|
| 523 | ]) |
---|
| 524 | ) |
---|
| 525 | |
---|
| 526 | gauss150 = Gauss( |
---|
| 527 | z=np.array([ |
---|
| 528 | -0.9998723404457334, |
---|
| 529 | -0.9993274305065947, |
---|
| 530 | -0.9983473449340834, |
---|
| 531 | -0.9969322929775997, |
---|
| 532 | -0.9950828645255290, |
---|
| 533 | -0.9927998590434373, |
---|
| 534 | -0.9900842691660192, |
---|
| 535 | -0.9869372772712794, |
---|
| 536 | -0.9833602541697529, |
---|
| 537 | -0.9793547582425894, |
---|
| 538 | -0.9749225346595943, |
---|
| 539 | -0.9700655145738374, |
---|
| 540 | -0.9647858142586956, |
---|
| 541 | -0.9590857341746905, |
---|
| 542 | -0.9529677579610971, |
---|
| 543 | -0.9464345513503147, |
---|
| 544 | -0.9394889610042837, |
---|
| 545 | -0.9321340132728527, |
---|
| 546 | -0.9243729128743136, |
---|
| 547 | -0.9162090414984952, |
---|
| 548 | -0.9076459563329236, |
---|
| 549 | -0.8986873885126239, |
---|
| 550 | -0.8893372414942055, |
---|
| 551 | -0.8795995893549102, |
---|
| 552 | -0.8694786750173527, |
---|
| 553 | -0.8589789084007133, |
---|
| 554 | -0.8481048644991847, |
---|
| 555 | -0.8368612813885015, |
---|
| 556 | -0.8252530581614230, |
---|
| 557 | -0.8132852527930605, |
---|
| 558 | -0.8009630799369827, |
---|
| 559 | -0.7882919086530552, |
---|
| 560 | -0.7752772600680049, |
---|
| 561 | -0.7619248049697269, |
---|
| 562 | -0.7482403613363824, |
---|
| 563 | -0.7342298918013638, |
---|
| 564 | -0.7198995010552305, |
---|
| 565 | -0.7052554331857488, |
---|
| 566 | -0.6903040689571928, |
---|
| 567 | -0.6750519230300931, |
---|
| 568 | -0.6595056411226444, |
---|
| 569 | -0.6436719971150083, |
---|
| 570 | -0.6275578900977726, |
---|
| 571 | -0.6111703413658551, |
---|
| 572 | -0.5945164913591590, |
---|
| 573 | -0.5776035965513142, |
---|
| 574 | -0.5604390262878617, |
---|
| 575 | -0.5430302595752546, |
---|
| 576 | -0.5253848818220803, |
---|
| 577 | -0.5075105815339176, |
---|
| 578 | -0.4894151469632753, |
---|
| 579 | -0.4711064627160663, |
---|
| 580 | -0.4525925063160997, |
---|
| 581 | -0.4338813447290861, |
---|
| 582 | -0.4149811308476706, |
---|
| 583 | -0.3959000999390257, |
---|
| 584 | -0.3766465660565522, |
---|
| 585 | -0.3572289184172501, |
---|
| 586 | -0.3376556177463400, |
---|
| 587 | -0.3179351925907259, |
---|
| 588 | -0.2980762356029071, |
---|
| 589 | -0.2780873997969574, |
---|
| 590 | -0.2579773947782034, |
---|
| 591 | -0.2377549829482451, |
---|
| 592 | -0.2174289756869712, |
---|
| 593 | -0.1970082295132342, |
---|
| 594 | -0.1765016422258567, |
---|
| 595 | -0.1559181490266516, |
---|
| 596 | -0.1352667186271445, |
---|
| 597 | -0.1145563493406956, |
---|
| 598 | -0.0937960651617229, |
---|
| 599 | -0.0729949118337358, |
---|
| 600 | -0.0521619529078925, |
---|
| 601 | -0.0313062657937972, |
---|
| 602 | -0.0104369378042598, |
---|
| 603 | 0.0104369378042598, |
---|
| 604 | 0.0313062657937972, |
---|
| 605 | 0.0521619529078925, |
---|
| 606 | 0.0729949118337358, |
---|
| 607 | 0.0937960651617229, |
---|
| 608 | 0.1145563493406956, |
---|
| 609 | 0.1352667186271445, |
---|
| 610 | 0.1559181490266516, |
---|
| 611 | 0.1765016422258567, |
---|
| 612 | 0.1970082295132342, |
---|
| 613 | 0.2174289756869712, |
---|
| 614 | 0.2377549829482451, |
---|
| 615 | 0.2579773947782034, |
---|
| 616 | 0.2780873997969574, |
---|
| 617 | 0.2980762356029071, |
---|
| 618 | 0.3179351925907259, |
---|
| 619 | 0.3376556177463400, |
---|
| 620 | 0.3572289184172501, |
---|
| 621 | 0.3766465660565522, |
---|
| 622 | 0.3959000999390257, |
---|
| 623 | 0.4149811308476706, |
---|
| 624 | 0.4338813447290861, |
---|
| 625 | 0.4525925063160997, |
---|
| 626 | 0.4711064627160663, |
---|
| 627 | 0.4894151469632753, |
---|
| 628 | 0.5075105815339176, |
---|
| 629 | 0.5253848818220803, |
---|
| 630 | 0.5430302595752546, |
---|
| 631 | 0.5604390262878617, |
---|
| 632 | 0.5776035965513142, |
---|
| 633 | 0.5945164913591590, |
---|
| 634 | 0.6111703413658551, |
---|
| 635 | 0.6275578900977726, |
---|
| 636 | 0.6436719971150083, |
---|
| 637 | 0.6595056411226444, |
---|
| 638 | 0.6750519230300931, |
---|
| 639 | 0.6903040689571928, |
---|
| 640 | 0.7052554331857488, |
---|
| 641 | 0.7198995010552305, |
---|
| 642 | 0.7342298918013638, |
---|
| 643 | 0.7482403613363824, |
---|
| 644 | 0.7619248049697269, |
---|
| 645 | 0.7752772600680049, |
---|
| 646 | 0.7882919086530552, |
---|
| 647 | 0.8009630799369827, |
---|
| 648 | 0.8132852527930605, |
---|
| 649 | 0.8252530581614230, |
---|
| 650 | 0.8368612813885015, |
---|
| 651 | 0.8481048644991847, |
---|
| 652 | 0.8589789084007133, |
---|
| 653 | 0.8694786750173527, |
---|
| 654 | 0.8795995893549102, |
---|
| 655 | 0.8893372414942055, |
---|
| 656 | 0.8986873885126239, |
---|
| 657 | 0.9076459563329236, |
---|
| 658 | 0.9162090414984952, |
---|
| 659 | 0.9243729128743136, |
---|
| 660 | 0.9321340132728527, |
---|
| 661 | 0.9394889610042837, |
---|
| 662 | 0.9464345513503147, |
---|
| 663 | 0.9529677579610971, |
---|
| 664 | 0.9590857341746905, |
---|
| 665 | 0.9647858142586956, |
---|
| 666 | 0.9700655145738374, |
---|
| 667 | 0.9749225346595943, |
---|
| 668 | 0.9793547582425894, |
---|
| 669 | 0.9833602541697529, |
---|
| 670 | 0.9869372772712794, |
---|
| 671 | 0.9900842691660192, |
---|
| 672 | 0.9927998590434373, |
---|
| 673 | 0.9950828645255290, |
---|
| 674 | 0.9969322929775997, |
---|
| 675 | 0.9983473449340834, |
---|
| 676 | 0.9993274305065947, |
---|
| 677 | 0.9998723404457334 |
---|
| 678 | ]), |
---|
| 679 | w=np.array([ |
---|
| 680 | 0.0003276086705538, |
---|
| 681 | 0.0007624720924706, |
---|
| 682 | 0.0011976474864367, |
---|
| 683 | 0.0016323569986067, |
---|
| 684 | 0.0020663664924131, |
---|
| 685 | 0.0024994789888943, |
---|
| 686 | 0.0029315036836558, |
---|
| 687 | 0.0033622516236779, |
---|
| 688 | 0.0037915348363451, |
---|
| 689 | 0.0042191661429919, |
---|
| 690 | 0.0046449591497966, |
---|
| 691 | 0.0050687282939456, |
---|
| 692 | 0.0054902889094487, |
---|
| 693 | 0.0059094573005900, |
---|
| 694 | 0.0063260508184704, |
---|
| 695 | 0.0067398879387430, |
---|
| 696 | 0.0071507883396855, |
---|
| 697 | 0.0075585729801782, |
---|
| 698 | 0.0079630641773633, |
---|
| 699 | 0.0083640856838475, |
---|
| 700 | 0.0087614627643580, |
---|
| 701 | 0.0091550222717888, |
---|
| 702 | 0.0095445927225849, |
---|
| 703 | 0.0099300043714212, |
---|
| 704 | 0.0103110892851360, |
---|
| 705 | 0.0106876814158841, |
---|
| 706 | 0.0110596166734735, |
---|
| 707 | 0.0114267329968529, |
---|
| 708 | 0.0117888704247183, |
---|
| 709 | 0.0121458711652067, |
---|
| 710 | 0.0124975796646449, |
---|
| 711 | 0.0128438426753249, |
---|
| 712 | 0.0131845093222756, |
---|
| 713 | 0.0135194311690004, |
---|
| 714 | 0.0138484622795371, |
---|
| 715 | 0.0141714592928592, |
---|
| 716 | 0.0144882814685445, |
---|
| 717 | 0.0147987907597169, |
---|
| 718 | 0.0151028518701744, |
---|
| 719 | 0.0154003323133401, |
---|
| 720 | 0.0156911024699895, |
---|
| 721 | 0.0159750356447283, |
---|
| 722 | 0.0162520081211971, |
---|
| 723 | 0.0165218992159766, |
---|
| 724 | 0.0167845913311726, |
---|
| 725 | 0.0170399700056559, |
---|
| 726 | 0.0172879239649355, |
---|
| 727 | 0.0175283451696437, |
---|
| 728 | 0.0177611288626114, |
---|
| 729 | 0.0179861736145128, |
---|
| 730 | 0.0182033813680609, |
---|
| 731 | 0.0184126574807331, |
---|
| 732 | 0.0186139107660094, |
---|
| 733 | 0.0188070535331042, |
---|
| 734 | 0.0189920016251754, |
---|
| 735 | 0.0191686744559934, |
---|
| 736 | 0.0193369950450545, |
---|
| 737 | 0.0194968900511231, |
---|
| 738 | 0.0196482898041878, |
---|
| 739 | 0.0197911283358190, |
---|
| 740 | 0.0199253434079123, |
---|
| 741 | 0.0200508765398072, |
---|
| 742 | 0.0201676730337687, |
---|
| 743 | 0.0202756819988200, |
---|
| 744 | 0.0203748563729175, |
---|
| 745 | 0.0204651529434560, |
---|
| 746 | 0.0205465323660984, |
---|
| 747 | 0.0206189591819181, |
---|
| 748 | 0.0206824018328499, |
---|
| 749 | 0.0207368326754401, |
---|
| 750 | 0.0207822279928917, |
---|
| 751 | 0.0208185680053983, |
---|
| 752 | 0.0208458368787627, |
---|
| 753 | 0.0208640227312962, |
---|
| 754 | 0.0208731176389954, |
---|
| 755 | 0.0208731176389954, |
---|
| 756 | 0.0208640227312962, |
---|
| 757 | 0.0208458368787627, |
---|
| 758 | 0.0208185680053983, |
---|
| 759 | 0.0207822279928917, |
---|
| 760 | 0.0207368326754401, |
---|
| 761 | 0.0206824018328499, |
---|
| 762 | 0.0206189591819181, |
---|
| 763 | 0.0205465323660984, |
---|
| 764 | 0.0204651529434560, |
---|
| 765 | 0.0203748563729175, |
---|
| 766 | 0.0202756819988200, |
---|
| 767 | 0.0201676730337687, |
---|
| 768 | 0.0200508765398072, |
---|
| 769 | 0.0199253434079123, |
---|
| 770 | 0.0197911283358190, |
---|
| 771 | 0.0196482898041878, |
---|
| 772 | 0.0194968900511231, |
---|
| 773 | 0.0193369950450545, |
---|
| 774 | 0.0191686744559934, |
---|
| 775 | 0.0189920016251754, |
---|
| 776 | 0.0188070535331042, |
---|
| 777 | 0.0186139107660094, |
---|
| 778 | 0.0184126574807331, |
---|
| 779 | 0.0182033813680609, |
---|
| 780 | 0.0179861736145128, |
---|
| 781 | 0.0177611288626114, |
---|
| 782 | 0.0175283451696437, |
---|
| 783 | 0.0172879239649355, |
---|
| 784 | 0.0170399700056559, |
---|
| 785 | 0.0167845913311726, |
---|
| 786 | 0.0165218992159766, |
---|
| 787 | 0.0162520081211971, |
---|
| 788 | 0.0159750356447283, |
---|
| 789 | 0.0156911024699895, |
---|
| 790 | 0.0154003323133401, |
---|
| 791 | 0.0151028518701744, |
---|
| 792 | 0.0147987907597169, |
---|
| 793 | 0.0144882814685445, |
---|
| 794 | 0.0141714592928592, |
---|
| 795 | 0.0138484622795371, |
---|
| 796 | 0.0135194311690004, |
---|
| 797 | 0.0131845093222756, |
---|
| 798 | 0.0128438426753249, |
---|
| 799 | 0.0124975796646449, |
---|
| 800 | 0.0121458711652067, |
---|
| 801 | 0.0117888704247183, |
---|
| 802 | 0.0114267329968529, |
---|
| 803 | 0.0110596166734735, |
---|
| 804 | 0.0106876814158841, |
---|
| 805 | 0.0103110892851360, |
---|
| 806 | 0.0099300043714212, |
---|
| 807 | 0.0095445927225849, |
---|
| 808 | 0.0091550222717888, |
---|
| 809 | 0.0087614627643580, |
---|
| 810 | 0.0083640856838475, |
---|
| 811 | 0.0079630641773633, |
---|
| 812 | 0.0075585729801782, |
---|
| 813 | 0.0071507883396855, |
---|
| 814 | 0.0067398879387430, |
---|
| 815 | 0.0063260508184704, |
---|
| 816 | 0.0059094573005900, |
---|
| 817 | 0.0054902889094487, |
---|
| 818 | 0.0050687282939456, |
---|
| 819 | 0.0046449591497966, |
---|
| 820 | 0.0042191661429919, |
---|
| 821 | 0.0037915348363451, |
---|
| 822 | 0.0033622516236779, |
---|
| 823 | 0.0029315036836558, |
---|
| 824 | 0.0024994789888943, |
---|
| 825 | 0.0020663664924131, |
---|
| 826 | 0.0016323569986067, |
---|
| 827 | 0.0011976474864367, |
---|
| 828 | 0.0007624720924706, |
---|
| 829 | 0.0003276086705538 |
---|
| 830 | ]) |
---|
| 831 | ) |
---|