git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@6378 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
245
lib/awpmd/ivutils/include/cerf.h
Normal file
245
lib/awpmd/ivutils/include/cerf.h
Normal file
@ -0,0 +1,245 @@
|
||||
# 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
|
||||
Reference in New Issue
Block a user