Homec4science
Diffusion lammps (efficient_neuronet) Browse

lammps/lib/awpmd/ivutils/include/cerf.h
7025ec968fc7efficient_neuronet

Last Change
Raw File

cerf.h
View Options

# ifndef CERF_H
# define CERF_H
# include <complex>
# include <cmath>
/*
Copyright (C) 1998, 1999 John Smith
This file is part of Octave.
Or it might be one day....
Octave is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2, or (at your option) any
later version.
Octave is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with Octave; see the file COPYING. If not, write to the Free
Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/
// Put together by John Smith john at arrows dot demon dot co dot uk,
// using ideas by others.
//
// Calculate erf(z) for complex z.
// Three methods are implemented; which one is used depends on z.
//
// The code includes some hard coded constants that are intended to
// give about 14 decimal places of accuracy. This is appropriate for
// 64-bit floating point numbers.
//
// Oct 1999: Fixed a typo that in
// const Complex cerf_continued_fraction( const Complex z )
// that caused erroneous answers for erf(z) where real(z) negative
//
//#include <math.h>
//#include <octave/oct.h>
//#include "f77-fcn.h"
//
// Abramowitz and Stegun: (eqn: 7.1.14) gives this continued
// fraction for erfc(z)
//
// erfc(z) = sqrt(pi).exp(-z^2). 1 1/2 1 3/2 2 5/2
// --- --- --- --- --- --- ...
// z + z + z + z + z + z +
//
// This is evaluated using Lentz's method, as described in the narative
// of Numerical Recipes in C.
//
// The continued fraction is true providing real(z)>0. In practice we
// like real(z) to be significantly greater than 0, say greater than 0.5.
//
template< class Complex>
const Complex cerfc_continued_fraction( const Complex z )
{
double tiny = 1e-20 ; // a small number, large enough to calculate 1/tiny
double eps = 1e-15 ; // large enough so that 1.0+eps > 1.0, when using
// the floating point arithmetic
//
// first calculate z+ 1/2 1
// --- --- ...
// z + z +
Complex f(z) ;
Complex C(f) ;
Complex D(0.0) ;
Complex delta ;
double a ;
a = 0.0 ;
do
{
a = a + 0.5 ;
D = z + a*D ;
C = z + a/C ;
if (D.real() == 0.0 && D.imag() == 0.0)
D = tiny ;
D = 1.0 / D ;
delta = (C * D) ;
f = f * delta ;
} while (abs(1.0-delta) > eps ) ;
//
// Do the first term of the continued fraction
//
f = 1.0 / f ;
//
// and do the final scaling
//
f = f * exp(-z*z)/ sqrt(M_PI) ;
return f ;
}
template< class Complex>
const Complex cerf_continued_fraction( const Complex z )
{
// warning("cerf_continued_fraction:");
if (z.real() > 0)
return 1.0 - cerfc_continued_fraction( z ) ;
else
return -1.0 + cerfc_continued_fraction( -z ) ;
}
//
// Abramawitz and Stegun, Eqn. 7.1.5 gives a series for erf(z)
// good for all z, but converges faster for smallish abs(z), say abs(z)<2.
//
template< class Complex>
const Complex cerf_series( const Complex z )
{
double tiny = 1e-20 ; // a small number compared with 1.
// warning("cerf_series:");
Complex sum(0.0) ;
Complex term(z) ;
Complex z2(z*z) ;
for (int n=0; n<3 || abs(term) > abs(sum)*tiny; n++)
{
sum = sum + term / (2*n+1) ;
term = -term * z2 / (n+1) ;
}
return sum * 2.0 / sqrt(M_PI) ;
}
//
// Numerical Recipes quotes a formula due to Rybicki for evaluating
// Dawson's Integral:
//
// exp(-x^2) integral exp(t^2).dt = 1/sqrt(pi) lim sum exp(-(z-n.h)^2) / n
// 0 to x h->0 n odd
//
// This can be adapted to erf(z).
//
template< class Complex>
const Complex cerf_rybicki( const Complex z )
{
// warning("cerf_rybicki:");
double h = 0.2 ; // numerical experiment suggests this is small enough
//
// choose an even n0, and then shift z->z-n0.h and n->n-h.
// n0 is chosen so that real((z-n0.h)^2) is as small as possible.
//
int n0 = 2*(int) (floor( z.imag()/(2*h) + 0.5 )) ;
Complex z0( 0.0, n0*h ) ;
Complex zp(z-z0) ;
Complex sum(0.0,0.0) ;
//
// limits of sum chosen so that the end sums of the sum are
// fairly small. In this case exp(-(35.h)^2)=5e-22
//
//
for (int np=-35; np<=35; np+=2)
{
Complex t( zp.real(), zp.imag()-np*h) ;
Complex b( exp(t*t) / (np+n0) ) ;
sum += b ;
}
sum = sum * 2 * exp(-z*z) / M_PI ;
return Complex(-sum.imag(), sum.real()) ;
}
template< class Complex>
const Complex cerf( const Complex z )
{
//
// Use the method appropriate to size of z -
// there probably ought to be an extra option for NaN z, or infinite z
//
//
if (abs(z) < 2.0)
return cerf_series( z ) ;
else if (abs(z.real()) < 0.5)
return cerf_rybicki( z ) ;
else
return cerf_continued_fraction( z ) ;
}
//
// Footnote:
//
// Using the definitions from Abramowitz and Stegun (7.3.1, 7.3.2)
// The fresnel intgerals defined as:
//
// / t=x
// C(x) = | cos(pi/2 t^2) dt
// /
// t=0
//
// and
// / t=x
// S(x) = | sin(pi/2 t^2) dt
// /
// t=0
//
// These can be derived from erf(x) using 7.3.22
//
// C(z) +iS(z) = (1+i) erf( sqrt(pi)/2 (1-i) z )
// -----
// 2
//
// --------------------------------------------------------------------------
// Some test examples -
// comparative data taken from Abramowitz and Stegun table 7.9.
// Table 7.9 tabulates w(z), where w(z) = exp(-z*z) erfc(iz)
// I have copied twelve values of w(z) from the table, and separately
// calculated them using this code. The results are identical.
//
// x y Abramowitz & Stegun | Octave Calculations
// w(x+iy) | w(x+iy) cerf ( i.(x+iy))
// 0.2 0.2 0.783538+0.157403i | 0.783538 +0.157403 0.23154672 -0.219516
// 0.2 0.7 0.515991+0.077275i | 0.515991 +0.077275 0.69741968 -0.138277
// 0.2 1.7 0.289309+0.027154i | 0.289309 +0.027154 0.98797507 -0.011744
// 0.2 2.7 0.196050+0.013002i | 0.196050 +0.013002 0.99994252 -0.000127
// 1.2 0.2 0.270928+0.469488i | 0.270928 +0.469488 0.90465623 -2.196064
// 1.2 0.7 0.280740+0.291851i | 0.280740 +0.291851 1.82926135 -0.639343
// 1.2 1.7 0.222436+0.129684i | 0.222436 +0.129684 1.00630308 +0.060067
// 1.2 2.7 0.170538+0.068617i | 0.170538 +0.068617 0.99955699 -0.000290
// 2.2 0.2 0.041927+0.287771i | 0.041927 +0.287771 24.70460755-26.205981
// 2.2 0.7 0.099943+0.242947i | 0.099943 +0.242947 9.88734713+18.310797
// 2.2 1.7 0.135021+0.153161i | 0.135021 +0.153161 1.65541359 -1.276707
// 2.2 2.7 0.127900+0.096330i | 0.127900 +0.096330 0.98619434 +0.000564
# endif
c4science · Help