autint.f

网络带宽测试工具

      Function autint(y, aa, bb, eps)               Result( result )
      Use      numerics
      Implicit None
! --------------------------------------------------------------------
      Real(l_)            :: y, aa, bb, eps
      Real(l_)            :: a, b, gint, ha, hanul, upper,
     &                       tot, tota, totb, xx, eps1, eps2, rot,
     &                       delta, totlag, difn, dift, verh,
     &                       result
      Integer             :: i
      Real(l_)            :: lenght, ma, mb, nauwk
      Real(l_), Parameter :: d1 = 0.90617984593866_l_,
     &                       d2 = 0.23692688505618_l_,
     &                       d3 = 0.53846931010568_l_,
     &                       d4 = 0.47862867049936_l_,
     &                       d5 = 0.56888888888888_l_
      Real(l_)            :: lc(10)
      Logical             :: deel, wissel
      Real(l_)            :: f(12)
! --------------------------------------------------------------------
      a = aa
      b = bb
      wissel = .FALSE.
      If ( a /= b ) Go To 10
      gint = 0.0_l_
      nauwk = 0.0_l_
      Go To 190
! --------------------------------------------------------------------
! --- The weights of the lagrangian formula.

10    lc(1) = 0.4347779016587e-1_l_
      lc(10) = lc(1)
      lc(2) = 0.8166287908853e-1_l_
      lc(9) = lc(2)
      lc(3) = 0.2369317573982_l_
      lc(8) = lc(3)
      lc(4) = 0.2054446657835_l_
      lc(7) = lc(4)
      lc(5) = 0.4324829075631_l_
      lc(6) = lc(5)
      If ( b >= a ) Go To 20
      wissel = .TRUE.
      upper = a
      a = b
      b = upper
20    upper = b
      deel = .FALSE.
      lenght = b - a
      hanul = lenght/1.0e6_l_
      nauwk = 0.0_l_
      gint = 0.0_l_
      f(1) = y(a)
30    ha = (upper-a)/2.0_l_
      If ( ha - hanul ) 40, 60, 60
40    upper = a + 2.0_l_*hanul
      ha = hanul
      If ( upper - b ) 60, 60, 50
50    ha = (b-a)/2.0_l_
      upper = b
60    Continue
      ma = (upper+a)/2.0_l_
      tot = 0.0_l_
      f(2) = y(-d1*ha+ma)
      tot = d2*f(2) + tot
      f(4) = y(-d3*ha+ma)
      tot = d4*f(4) + tot
      tot = d5*y(ma) + tot
      f(7) = y(d3*ha+ma)
      tot = d4*f(7) + tot
      f(9) = y(d1*ha+ma)
      tot = d2*f(9) + tot
70    f(10) = y(upper)
      ha = ha/2.0_l_
      ma = (3.0_l_*a+upper)/4.0_l_
      mb = (3.0_l_*upper+a)/4.0_l_
      tota = 0.0_l_
      totb = 0.0_l_
      xx = -d1*ha
      f(11) = y(xx+ma)
      tota = d2*f(11) + tota
      totb = d2*y(xx+mb) + totb
      xx = -d3*ha
      f(3) = y(xx+ma)
      tota = d4*f(3) + tota
      f(6) = y(xx+mb)
      totb = d4*f(6) + totb
      tota = d5*y(ma) + tota
      totb = d5*y(mb) + totb
      f(5) = y(-xx+ma)
      tota = d4*f(5) + tota
      f(8) = y(-xx+mb)
      totb = d4*f(8) + totb
      xx = d1*ha
      f(12) = y(xx+ma)
      tota = d2*f(12) + tota
      totb = d2*y(xx+mb) + totb
      totb = tota + totb
      dift = totb - tot*2.0_l_
      eps1 = 2.0_l_*eps/lenght
      eps2 = 1.e-10_l_*abs(totb)/2.0_l_
      If ( eps1 < eps2 ) eps1 = eps2
      If ( eps1 /= 0.0_l_ ) Go To 80
      rot = 2.0_l_
      Go To 100
80    delta = 2046.0_l_*eps1
      rot = abs(dift)
      If ( rot <= delta/1023.0_l_ ) Go To 90
      rot = delta/rot
      Go To 100
90    rot = 1024.0_l_
100   If ( rot < 1.0_l_ .AND. ha > hanul/2.0_l_ ) Go To 170
      totlag = 0.0_l_
      Do i = 1, 10
         totlag = lc(i)*f(i) + totlag
      End Do
      difn = totlag*2.0_l_ - totb
      If ( Abs( difn*dift )<= 1.E-20_l_*totb**2.0_l_) Go To 130
      If ( difn*dift < 0.0_l_ ) Go To 120
      If ( difn /= 0.0_l_ ) verh = dift/difn
      If ( verh >= 2.0_l_ .AND. verh <= 2.3_l_ ) Go To 150
      If ( Abs( tot - totlag ) < eps1/2.0_l_ ) Go To 150
120   If ( ha > hanul/2.0_l_ ) Go To 140
130   nauwk = nauwk + Abs( totlag - tot )*ha*2.0_l_
      Go To 160
140   tot = tota
      upper = ( upper + a )/2.0_l_
      f(2) = f(11)
      f(4) = f(3)
      f(7) = f(5)
      f(9) = f(12)
      Go To 70
150   nauwk = nauwk + eps1*ha*2.0_l_
160   gint = totb*ha + gint
      If ( upper == b ) Go To 190
      xx = upper
      deel = .TRUE.
      upper = rot**(+0.1_l_)*(upper-a)*0.9_l_ + upper
      a = xx
      If ( upper > b ) upper = b
      f(1) = f(10)
      Go To 30
170   COntinue
      If ( .NOT. deel ) Go To 180
      upper = rot**(0.1_l_)*(upper-a)*0.9_l_ + a
      Go To 30
180   tot = tota
      upper = (upper+a)/2.0_l_
      f(2) = f(11)
      f(4) = f(3)
      f(7) = f(5)
      f(9) = f(12)
      Go To 70
190   If ( wissel ) gint = -gint
      result = gint
      eps = nauwk
! --------------------------------------------------------------------
      End Function autint