C=====Integrated Function===============================================
      REAL FUNCTION fun(x)
      REAL x
      !2*x
      !3*x**2
      !4*x**3
      !5*x**4
      !6*x**5
      !7*x**6
      !8*x**7
      !exp**(sin(3.1415926*(15*x)/180.0)+15*cos(3.1415926*(15*x)/180.0)
      fun=4*x**3
      END

C====Read parametres====================================================
      SUBROUTINE InputParams
      REAL a, b, dx
      COMMON/params/a, b, dx, n_steps
      OPEN(1, FILE = 'params.txt')
      READ(1, *)a, b, n_steps
      dx = (b - a)/n_steps
      CLOSE(1)
      END

C====Create grid========================================================
      SUBROUTINE CreateGrid(grid, n_steps)
      DIMENSION grid(n_steps + 1)
      CHARACTER f*16, s*16
      COMMON/params/a, b, dx
      n_steps = n_steps + 1
      point = a
      grid(1) = point
      DO i = 2, n_steps
        point = point + dx
        grid(i) = point
      END DO
      END

C=====Newtone-Cotes=====================================================
      REAL FUNCTION CotesNewtone(grid, n_steps)
      COMMON/params/a, b
      DIMENSION grid(n_steps)
      DIMENSION coef(6)
      DATA coef / 19.0, 75.0, 50.0, 50.0, 75.0, 19.0 /
      sum = 0.0
      DO i = 1, n_steps - 1
        a1 = grid(i)
        b1 = grid(i + 1)
              h = (b1 - a1)/5.0
              sum_int = 0.0
              DO j = 1, 6
                sum_int = sum_int + coef(j)*fun(a1 + (j - 1)*h)
              END DO
              sum = sum + sum_int*(b1 - a1)*(1.0/288.0)
      END DO
      CotesNewtone = sum
      END

C=====Gauss=============================================================
      REAL FUNCTION Gauss(grid, n_steps)
      DIMENSION grid(n_steps), q(3),x(3)
      REAL q, x, h, sum, sum_int, a, b
      COMMON/params/a, b
      DATA x / -0.7745967, 0.0, 0.7745967 /
      DATA q / 0.5555556, 0.8888889, 0.5555556 /
      h = (b-a)/(n_steps - 1)
      DO i = 1,3
        sum_int = 0.0
        DO j = 2, n_steps
            x1 = (grid(j-1)+grid(j)+x(i)*h)/2
            sum_int=sum_int+h*fun(x1)
        END DO
        sum = sum + q(i) * sum_int
      END DO
      Gauss = sum/2.0
      END

C=====Method of trapezoids==============================================
      REAL FUNCTION Trap(grid, n_steps)
      DIMENSION grid(n_steps)
      COMMON/params/a, b, dx
      sum = 0.0
      DO i = 2, n_steps - 1
        sum = sum + fun(grid(i))
      END DO
      Trap = (sum + (fun(grid(1)) + fun(grid(n_steps)))/2.0)*dx
      END

C=====Main function=====================================================
      PROGRAM Main
      REAL array
      DIMENSION array(1000000)
      REAL a, b, dx, ans
      COMMON/params/a, b, dx, n_steps
      CALL InputParams
      CALL CreateGrid(array(1), n_steps)
      ans = Trap(array(1), n_steps)
      PRINT *, ans
      ans = CotesNewtone(array(1), n_steps)
      PRINT *, ans
      ans = Gauss(array(1), n_steps)
      PRINT *, ans
      PAUSE
      END