[cig-commits] r8473 - seismo/2D/SPECFEM2D/trunk
walter at geodynamics.org
walter at geodynamics.org
Fri Dec 7 15:48:59 PST 2007
Author: walter
Date: 2007-12-07 15:48:59 -0800 (Fri, 07 Dec 2007)
New Revision: 8473
Added:
seismo/2D/SPECFEM2D/trunk/numerical_recipes.f90
Log:
added SPECFEM2D/SEM_2D_Dimitri/numerical_recipes.f90
Added: seismo/2D/SPECFEM2D/trunk/numerical_recipes.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/numerical_recipes.f90 2005-02-24 11:24:59 UTC (rev 8472)
+++ seismo/2D/SPECFEM2D/trunk/numerical_recipes.f90 2007-12-07 23:48:59 UTC (rev 8473)
@@ -0,0 +1,169 @@
+!=====================================================================
+!
+! Routines from "Numerical Recipes: the Art of Scientific Computing"
+! W. H. Press et al., Cambridge University Press
+!
+!=====================================================================
+
+! double precision routines
+
+ double precision function erf(x)
+
+ implicit none
+
+ double precision x
+
+! this routine uses routine gammp
+ double precision gammp
+
+ if(x<0.)then
+ erf=-gammp(0.5d0,x**2)
+ else
+ erf=gammp(0.5d0,x**2)
+ endif
+
+ end function erf
+
+! ---------------------------------
+
+ double precision function gammp(a,x)
+
+ implicit none
+
+ double precision a,x
+
+! this routine uses routines gcf and gser
+ double precision gammcf,gamser,gln
+
+ if(x<0.d0 .or. a <= 0.d0) stop 'bad arguments in gammp'
+
+ if(x<a+1.d0)then
+ call gser(gamser,a,x,gln)
+ gammp=gamser
+ else
+ call gcf(gammcf,a,x,gln)
+ gammp=1.d0-gammcf
+ endif
+
+ end function gammp
+
+! ---------------------------------
+
+ subroutine gcf(gammcf,a,x,gln)
+
+ implicit none
+
+ double precision a,gammcf,gln,x
+
+ double precision, parameter :: EPS=3.d-7,FPMIN=1.d-30
+ integer, parameter :: ITMAX=100
+
+! this routine uses routine gammln
+
+ integer i
+ double precision an,b,c,d,del,h
+
+ double precision, external :: gammln
+
+ gln=gammln(a)
+ b=x+1.d0-a
+ c=1.d0/FPMIN
+ d=1.d0/b
+ h=d
+ do i=1,ITMAX
+ an=-i*(i-a)
+ b=b+2.d0
+ d=an*d+b
+ if(dabs(d)<FPMIN)d=FPMIN
+ c=b+an/c
+ if(dabs(c)<FPMIN)c=FPMIN
+ d=1.d0/d
+ del=d*c
+ h=h*del
+ if(dabs(del-1.d0)<EPS) then
+ gammcf=exp(-x+a*log(x)-gln)*h
+ return
+ endif
+ enddo
+
+ stop 'a too large, ITMAX too small in gcf'
+
+ end subroutine gcf
+
+! ---------------------------------
+
+ subroutine gser(gamser,a,x,gln)
+
+ implicit none
+
+ double precision a,gamser,gln,x
+
+ integer, parameter :: ITMAX=100
+ double precision, parameter :: EPS=3.d-7
+
+! this routine uses routine gammln
+
+ integer n
+ double precision ap,del,sumval
+
+ double precision, external :: gammln
+
+ gln=gammln(a)
+
+ if(x <= 0.d0)then
+ if(x<0.d0) stop 'x < 0 in gser'
+ gamser=0.d0
+ return
+ endif
+
+ ap=a
+ sumval=1.d0/a
+ del=sumval
+
+ do n=1,ITMAX
+ ap=ap+1.d0
+ del=del*x/ap
+ sumval=sumval+del
+ if(dabs(del)<dabs(sumval)*EPS) then
+ gamser=sumval*exp(-x+a*log(x)-gln)
+ return
+ endif
+ enddo
+
+ stop 'a too large, ITMAX too small in gser'
+
+ end subroutine gser
+
+! ---------------------------------
+
+ double precision function gammln(xx)
+
+ implicit none
+
+ double precision xx
+
+ integer j
+ double precision ser,stp,tmp,x,y,cof(6)
+
+ cof(1) = 76.18009172947146d0
+ cof(2) = -86.50532032941677d0
+ cof(3) = 24.01409824083091d0
+ cof(4) = -1.231739572450155d0
+ cof(5) = 0.1208650973866179d-2
+ cof(6) = -0.5395239384953d-5
+
+ stp = 2.5066282746310005d0
+
+ x=xx
+ y=x
+ tmp=x+5.5d0
+ tmp=(x+0.5d0)*log(tmp)-tmp
+ ser=1.000000000190015d0
+ do j=1,6
+ y=y+1.d0
+ ser=ser+cof(j)/y
+ enddo
+ gammln=tmp+log(stp*ser/x)
+
+ end function gammln
+
More information about the cig-commits
mailing list