← Back to branch summary

~maddevelopers/mg5amcnlo/new_clustering

« back to all changes in this revision

Viewing changes to Template/NLO/Source/derivative.f

  • Committer: Rikkert Frederix
  • Date: 2021-09-09 15:51:40 UTC
  • mfrom: (78.75.502 3.2.1)
  • Revision ID: frederix@physik.uzh.ch-20210909155140-rg6umfq68h6h47cf
merge with 3.2.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/Source/derivative.f
 
1
      function derivative(fun,xi,hi,idiri,xmini,xmaxi,erroro)
 
2
c returns the derivative of f
 
3
c
 
4
c fun    (INPUT)= function to differentiate
 
5
c xi     (INPUT)= point where to compute the derivative
 
6
c hi     (INPUT)= initial stepsize
 
7
c idiri  (INPUT)= type of derivative (see below)
 
8
c xmini  (INPUT)= lower bound in x
 
9
c xmaxi  (INPUT)= upper bound in x
 
10
c error  (OUTPUT)= estimated error
 
11
c
 
12
c idiri =-1 --> right derivative
 
13
c idiri = 0 --> both-sides derivative
 
14
c idiri = 1 --> left derivative
 
15
c idiri = 2 --> higher-order both-sides left derivative
 
16
      implicit none
 
17
      real*8 derivative,fun,xi,hi,xmini,xmaxi,erroro
 
18
      integer idiri
 
19
c
 
20
      integer idir,i,j,itn,jmax
 
21
      parameter(itn=2,jmax=2)
 
22
c itn = 2 guarantees a fast code. itns >> 2 gives more precise
 
23
c results in extreme configurations, but slows the code down
 
24
c linearly
 
25
      real*8 x,h,hh,error,con,big,dd,
 
26
     &xmin,xmax,tiny,res(0:itn,jmax),er(0:itn)
 
27
      parameter(con=3d0,big=1d40)
 
28
      parameter(tiny=1.d-8)
 
29
      external fun
 
30
c
 
31
      x=xi
 
32
      h=hi
 
33
      xmin=xmini
 
34
      xmax=xmaxi
 
35
      idir=idiri
 
36
      if(h.eq.0d0)then
 
37
         write(*,*)'Error #1 in function derivative'
 
38
         stop
 
39
      endif
 
40
      if(idir.lt.-1.or.idir.gt.2)then
 
41
         write(*,*)'Error #2 in function derivative',idir
 
42
         stop
 
43
      endif
 
44
      if(xmin.ge.xmax)then
 
45
         write(*,*)'Error #3 in function derivative',x,xmin,xmax
 
46
         stop
 
47
      endif
 
48
c Set the derivative equal to zero and error equal to one if outside range
 
49
      if(x.lt.xmin.or.x.gt.xmax)then
 
50
        derivative=0.d0
 
51
        erroro=1.d0
 
52
        return
 
53
      endif
 
54
c If close to borders of range, use left or right first-order derivative;
 
55
c may include higher orders if implemented from the left or right only
 
56
      if(x.le.(xmin+tiny))then
 
57
        x=xmin+tiny
 
58
        idir=-1
 
59
      elseif(x.ge.(xmax-tiny))then
 
60
        x=xmax-tiny
 
61
        idir=1
 
62
      endif
 
63
c
 
64
      h=min(h,min(xmax-x,x-xmin))
 
65
      error=big
 
66
      dd=0.d0
 
67
      do i=0,itn
 
68
         hh=h/1d1**(6d0*i/itn)
 
69
         do j=1,jmax
 
70
            if(idir.eq.0)then
 
71
               res(i,j)=(fun(x+hh)-fun(x-hh))/(2d0*hh)
 
72
            elseif(idir.eq.-1)then
 
73
               res(i,j)=(fun(x+hh)-fun(x))/hh
 
74
            elseif(idir.eq.1)then
 
75
               res(i,j)=(fun(x)-fun(x-hh))/hh
 
76
            elseif(idir.eq.2)then
 
77
c second-order formula: yet higher-order formulae exist
 
78
c but they require too many evaluations of f, which slows
 
79
c the code down quite significantly
 
80
               res(i,j)=(fun(x-2*hh)-8*fun(x-hh)
 
81
     &                  -fun(x+2*hh)+8*fun(x+hh))/(12d0*hh)
 
82
            endif
 
83
            if(j.lt.jmax)hh=hh/con
 
84
         enddo
 
85
         er(i)=abs(abs(res(i,1))-abs(res(i,2)))/
 
86
     &         max(abs(res(i,1)),abs(res(i,2)))
 
87
         if(er(i).lt.error.and.er(i).ne.0d0)then
 
88
            error=er(i)
 
89
            dd=(res(i,1)+res(i,2))/2d0
 
90
         endif
 
91
      enddo
 
92
      derivative=dd
 
93
      erroro=error*2*abs(dd)
 
94
c
 
95
      return
 
96
      end