TABLE OF CONTENTS
ABINIT/GAMMA_FUNCTION [ Functions ]
NAME
GAMMA_FUNCTION
FUNCTION
COPYRIGHT
Copyright (C) 2014-2018 ABINIT group (XG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
PARENTS
pspnl_hgh_rec,pspnl_operat_rec,vso_realspace_local
CHILDREN
gsl_f90_sf_gamma
SOURCE
30 #if defined HAVE_CONFIG_H 31 #include "config.h" 32 #endif 33 34 #include "abi_common.h" 35 36 37 subroutine GAMMA_FUNCTION(X,GA) 38 39 use defs_basis 40 41 !This section has been created automatically by the script Abilint (TD). 42 !Do not modify the following lines by hand. 43 #undef ABI_FUNC 44 #define ABI_FUNC 'GAMMA_FUNCTION' 45 !End of the abilint section 46 47 implicit none 48 49 #ifdef HAVE_GSL 50 ! in case we have gsl, no need to use explicit function, just wrap the 51 ! call to the GSL C function in 01_gsl_ext/ 52 53 ! arguments 54 55 real(dp),intent(in) :: x 56 real(dp),intent(out) :: ga 57 58 call gsl_f90_sf_gamma(x,ga) 59 60 #else 61 ! ==================================================== 62 ! Purpose: This program computes the gamma function 63 ! Gamma(x) using subroutine GAMMA 64 ! Examples: 65 ! x Gamma(x) 66 ! ---------------------------- 67 ! 1/3 2.678938534708 68 ! 0.5 1.772453850906 69 ! -0.5 -3.544907701811 70 ! -1.5 2.363271801207 71 ! 5.0 24.000000000000 72 ! ==================================================== 73 ! 74 ! This routine was downloaded from UIUC: 75 ! http://jin.ece.uiuc.edu/routines/routines.html 76 ! 77 ! The programs appear to accompany a book "Computation of Special 78 ! Functions" (1996) John Wiley and Sons, but are distributed online 79 ! by the authors. Exact copyright should be checked. 80 ! 81 ! Authors / copyright: 82 ! Shanjie Zhang and Jianming Jin 83 ! Proposed contact is: j-jin1@uiuc.edu 84 ! 85 ! 20 October 2008: 86 ! Incorporated into ABINIT by M. Verstraete 87 ! 88 ! 89 ! 90 ! ================================================== 91 ! Purpose: Compute the gamma function Gamma(x) 92 ! Input : x --- Argument of Gamma(x) 93 ! ( x is not equal to 0,-1,-2, etc ) 94 ! Output: GA --- Gamma(x) 95 ! ================================================== 96 ! 97 98 ! arguments 99 100 real(dp),intent(in) :: x 101 real(dp),intent(out) :: ga 102 103 ! local variables 104 integer :: k,m 105 real(dp) :: m1,z,r,gr 106 real(dp) :: G(26) 107 108 ! source code: 109 110 ! initialization of reference data 111 G=(/1.0D0,0.5772156649015329D0, & 112 & -0.6558780715202538D0, -0.420026350340952D-1, & 113 & 0.1665386113822915D0,-.421977345555443D-1, & 114 & -.96219715278770D-2, .72189432466630D-2, & 115 & -.11651675918591D-2, -.2152416741149D-3, & 116 & .1280502823882D-3, -.201348547807D-4, & 117 & -.12504934821D-5, .11330272320D-5, & 118 & -.2056338417D-6, .61160950D-8, & 119 & .50020075D-8, -.11812746D-8, & 120 & .1043427D-9, .77823D-11, & 121 & -.36968D-11, .51D-12, & 122 & -.206D-13, -.54D-14, .14D-14, .1D-15/) 123 124 125 ! for the integer case, do explicit factorial 126 if (X==int(X)) then 127 if (X > 0.0D0) then 128 GA=1.0D0 129 M1=X-1 130 do K=2,int(M1) 131 GA=GA*K 132 end do 133 else 134 GA=1.0D+300 135 end if 136 ! for the integer case, do explicit factorial 137 else 138 if (abs(X) > 1.0D0) then 139 Z=abs(X) 140 M=int(Z) 141 R=1.0D0 142 do K=1,M 143 R=R*(Z-K) 144 end do 145 Z=Z-M 146 else 147 Z=X 148 end if 149 GR=G(26) 150 do K=25,1,-1 151 GR=GR*Z+G(K) 152 end do 153 GA=1.0D0/(GR*Z) 154 if (abs(X) > 1.0D0) then 155 GA=GA*R 156 if (X < 0.0D0) GA=-PI/(X*GA*SIN(PI*X)) 157 end if 158 end if 159 return 160 161 #endif 162 ! end preproc for presence of GSL 163 164 end subroutine GAMMA_FUNCTION