From f61865f4638158170082dca249903e6f846da6be Mon Sep 17 00:00:00 2001 From: daniau Date: Fri, 29 Jun 2007 13:36:06 +0000 Subject: [PATCH] git-svn-id: https://lxsd.femto-st.fr/svn/fvn@14 b657c933-2333-4658-acf2-d3c7c2708721 --- fvnlib.f90 | 196 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 195 insertions(+), 1 deletion(-) diff --git a/fvnlib.f90 b/fvnlib.f90 index c4eceac..e426022 100644 --- a/fvnlib.f90 +++ b/fvnlib.f90 @@ -22,7 +22,8 @@ module fvn ! if you find it usefull it would be kind to give credits ;-) Nevertheless, you ! may give credits to quadpack authors. ! -! Version 1.1 +! svn version +! June 2007 : added some complex trigonometric functions ! ! TO DO LIST : ! + Order eigenvalues and vectors in decreasing eigenvalue's modulus order -> atm @@ -36,6 +37,11 @@ module fvn !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! implicit none +! We define pi and i for the module +real(kind=8),parameter :: fvn_pi = 3.141592653589793_8 +complex(kind=8),parameter :: fvn_i = (0._8,1._8) + + ! All quadpack routines are private to the module private :: d1mach,dqag,dqag_2d_inner,dqag_2d_outer,dqage,dqage_2d_inner, & dqage_2d_outer,dqk15,dqk15_2d_inner,dqk15_2d_outer,dqk21,dqk21_2d_inner,dqk21_2d_outer, & @@ -1772,6 +1778,194 @@ include "fvn_quadpack/dqk61_2d_inner.f" include "fvn_quadpack/dqk21_2d_inner.f" include "fvn_quadpack/dqk15_2d_outer.f" +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! +! Trigonometric functions +! +! fvn_z_acos, fvn_z_asin : complex arc cosine and sine +! fvn_d_acosh : arc cosinus hyperbolic +! fvn_d_asinh : arc sinus hyperbolic +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +function fvn_z_acos(z) + ! double complex arccos function adapted from + ! the c gsl library + ! http://www.gnu.org/software/gsl/ + implicit none + complex(kind=8) :: fvn_z_acos + complex(kind=8) :: z + real(kind=8) :: rz,iz,x,y,a,b,y2,r,s,d,apx,am1 + real(kind=8),parameter :: a_crossover=1.5_8,b_crossover = 0.6417_8 + complex(kind=8),parameter :: i=(0._8,1._8) + real(kind=8) :: r_res,i_res + + rz=dble(z) + iz=aimag(z) + if ( iz == 0._8 ) then + fvn_z_acos=fvn_z_acos_real(rz) + return + end if + + x=dabs(rz) + y=dabs(iz) + r=fvn_d_hypot(x+1.,y) + s=fvn_d_hypot(x-1.,y) + a=0.5*(r + s) + b=x/a + y2=y*y + + if (b <= b_crossover) then + r_res=dacos(b) + else + if (x <= 1.) then + d=0.5*(a+x)*(y2/(r+x+1)+(s + (1 - x))) + r_res=datan(dsqrt(d)/x) + else + apx=a+x + d=0.5*(apx/(r+x+1)+apx/(s + (x - 1))) + r_res=datan((y*dsqrt(d))/x); + end if + end if + + if (a <= a_crossover) then + if (x < 1.) then + am1=0.5*(y2 / (r + (x + 1)) + y2 / (s + (1 - x))) + else + am1=0.5*(y2 / (r + (x + 1)) + (s + (x - 1))) + end if + i_res = dlog(1.+(am1 + sqrt (am1 * (a + 1)))); + else + i_res = dlog (a + dsqrt (a*a - 1.)); + end if + if (rz <0.) then + r_res=fvn_pi-r_res + end if + i_res=-sign(1._8,iz)*i_res + fvn_z_acos=dcmplx(r_res)+fvn_i*dcmplx(i_res) + +end function fvn_z_acos + +function fvn_z_acos_real(r) + ! return the double complex arc cosinus for a + ! double precision argument + implicit none + real(kind=8) :: r + complex(kind=8) :: fvn_z_acos_real + + if (dabs(r)<=1._8) then + fvn_z_acos_real=dcmplx(dacos(r)) + return + end if + if (r < 0._8) then + fvn_z_acos_real=dcmplx(fvn_pi)-fvn_i*dcmplx(fvn_d_acosh(-r)) + else + fvn_z_acos_real=fvn_i*dcmplx(fvn_d_acosh(r)) + end if +end function + + +function fvn_z_asin(z) + ! double complex arcsin function derived from + ! the c gsl library + ! http://www.gnu.org/software/gsl/ + implicit none + complex(kind=8) :: fvn_z_asin + complex(kind=8) :: z + real(kind=8) :: rz,iz,x,y,a,b,y2,r,s,d,apx,am1 + real(kind=8),parameter :: a_crossover=1.5_8,b_crossover = 0.6417_8 + real(kind=8) :: r_res,i_res + + rz=dble(z) + iz=aimag(z) + if ( iz == 0._8 ) then + ! z is real + fvn_z_asin=fvn_z_asin_real(rz) + return + end if + + x=dabs(rz) + y=dabs(iz) + r=fvn_d_hypot(x+1.,y) + s=fvn_d_hypot(x-1.,y) + a=0.5*(r + s) + b=x/a + y2=y*y + + if (b <= b_crossover) then + r_res=dasin(b) + else + if (x <= 1.) then + d=0.5*(a+x)*(y2/(r+x+1)+(s + (1 - x))) + r_res=datan(x/dsqrt(d)) + else + apx=a+x + d=0.5*(apx/(r+x+1)+apx/(s + (x - 1))) + r_res=datan(x/(y*dsqrt(d))); + end if + end if + + if (a <= a_crossover) then + if (x < 1.) then + am1=0.5*(y2 / (r + (x + 1)) + y2 / (s + (1 - x))) + else + am1=0.5*(y2 / (r + (x + 1)) + (s + (x - 1))) + end if + i_res = dlog(1.+(am1 + sqrt (am1 * (a + 1)))); + else + i_res = dlog (a + dsqrt (a*a - 1.)); + end if + r_res=sign(1._8,rz)*r_res + i_res=sign(1._8,iz)*i_res + fvn_z_asin=dcmplx(r_res)+fvn_i*dcmplx(i_res) + +end function fvn_z_asin + +function fvn_z_asin_real(r) + ! return the double complex arc sinus for a + ! double precision argument + implicit none + real(kind=8) :: r + complex(kind=8) :: fvn_z_asin_real + + if (dabs(r)<=1._8) then + fvn_z_asin_real=dcmplx(dasin(r)) + return + end if + if (r < 0._8) then + fvn_z_asin_real=dcmplx(-fvn_pi/2._8)+fvn_i*dcmplx(fvn_d_acosh(-r)) + else + fvn_z_asin_real=dcmplx(fvn_pi/2._8)-fvn_i*dcmplx(fvn_d_acosh(r)) + end if +end function fvn_z_asin_real + +function fvn_d_acosh(r) + ! return the arc hyperbolic cosine + implicit none + real(kind=8) :: r + real(kind=8) :: fvn_d_acosh + if (r >=1) then + fvn_d_acosh=dlog(r+dsqrt(r*r-1)) + else + !! TODO : Better error handling!!!!!! + stop "Argument to fvn_d_acosh lesser than 1" + end if +end function fvn_d_acosh + +function fvn_d_asinh(r) + ! return the arc hyperbolic sine + implicit none + real(kind=8) :: r + real(kind=8) :: fvn_d_asinh + fvn_d_asinh=dlog(r+dsqrt(r*r+1)) +end function fvn_d_asinh + +function fvn_d_hypot(a,b) + implicit none + ! return the euclidian norm of vector(a,b) + real(kind=8) :: a,b + real(kind=8) :: fvn_d_hypot + fvn_d_hypot=dsqrt(a*a+b*b) +end function -- 2.16.4