./32_util/ctrap.F90

ABINIT/ctrap

ABINIT/ctrap [ Functions ]

NAME

 ctrap

FUNCTION

 Do corrected trapezoidal integral on uniform grid of spacing hh.

 Copyright (C) 1998-2017 ABINIT group (DCA, XG, FrD)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

  imax=highest index of grid=grid point number of upper limit
  ff(imax)=integrand values
  hh=spacing between x points

OUTPUT

  ans=resulting integral by corrected trapezoid

NOTES

PARENTS

      psden,psp11lo,psp11nl,psp5lo,psp5nl,psp8lo,psp8nl,vhtnzc

CHILDREN

SOURCE

 33 #if defined HAVE_CONFIG_H
 34 #include "config.h"
 35 #endif
 36 
 37 #include "abi_common.h"
 38 
 39 
 40 subroutine ctrap(imax,ff,hh,ans)
 41 
 42  use defs_basis
 43 
 44 !This section has been created automatically by the script Abilint (TD).
 45 !Do not modify the following lines by hand.
 46 #undef ABI_FUNC
 47 #define ABI_FUNC 'ctrap'
 48 !End of the abilint section
 49 
 50  implicit none
 51 
 52 !Arguments ------------------------------------
 53 !scalars
 54  integer,intent(in) :: imax
 55  real(dp),intent(in) :: hh
 56  real(dp),intent(out) :: ans
 57 !arrays
 58  real(dp),intent(in) :: ff(imax)
 59 
 60 !Local variables-------------------------------
 61 !scalars
 62  integer :: ir,ir2
 63  real(dp) :: endpt,sum
 64 
 65 ! *************************************************************************
 66 
 67  if (imax>=10)then
 68 
 69 !  endpt=end point correction terms (low and high ends)
 70    endpt  = (23.75d0*(ff(1)+ff(imax  )) &
 71 &   + 95.10d0*(ff(2)+ff(imax-1)) &
 72 &   + 55.20d0*(ff(3)+ff(imax-2)) &
 73 &   + 79.30d0*(ff(4)+ff(imax-3)) &
 74 &   + 70.65d0*(ff(5)+ff(imax-4)))/ 72.d0
 75    ir2 = imax - 5
 76    sum=0.00d0
 77    if (ir2 > 5) then
 78      do ir=6,ir2
 79        sum = sum + ff(ir)
 80      end do
 81    end if
 82    ans = (sum + endpt ) * hh
 83 
 84  else if (imax>=8)then
 85    endpt  = (17.0d0*(ff(1)+ff(imax  )) &
 86 &   + 59.0d0*(ff(2)+ff(imax-1)) &
 87 &   + 43.0d0*(ff(3)+ff(imax-2)) &
 88 &   + 49.0d0*(ff(4)+ff(imax-3)) )/ 48.d0
 89    sum=0.0d0
 90    if(imax==9)sum=ff(5)
 91    ans = (sum + endpt ) * hh
 92 
 93  else if (imax==7)then
 94    ans = (17.0d0*(ff(1)+ff(imax  )) &
 95 &   + 59.0d0*(ff(2)+ff(imax-1)) &
 96 &   + 43.0d0*(ff(3)+ff(imax-2)) &
 97 &   + 50.0d0* ff(4)                )/ 48.d0  *hh
 98 
 99  else if (imax==6)then
100    ans = (17.0d0*(ff(1)+ff(imax  )) &
101 &   + 59.0d0*(ff(2)+ff(imax-1)) &
102 &   + 44.0d0*(ff(3)+ff(imax-2)) )/ 48.d0  *hh
103 
104  else if (imax==5)then
105    ans = (     (ff(1)+ff(5)) &
106 &   + four*(ff(2)+ff(4)) &
107 &   + two * ff(3)         )/ three  *hh
108 
109  else if (imax==4)then
110    ans = (three*(ff(1)+ff(4)) &
111 &   + nine *(ff(2)+ff(3))  )/ eight  *hh
112 
113  else if (imax==3)then
114    ans = (     (ff(1)+ff(3)) &
115 &   + four* ff(2)         )/ three *hh
116 
117  else if (imax==2)then
118    ans = (ff(1)+ff(2))/ two  *hh
119 
120  else if (imax==1)then
121    ans = ff(1)*hh
122 
123  end if
124 
125 end subroutine ctrap

TABLE OF CONTENTS

ABINIT/ctrap [ Functions ]