PROGRAM INTEGRATION
  IMPLICIT NONE
  REAL(8), PARAMETER :: PI = 4.0D0 * ATAN(1.0D0)
  INTEGER :: I,NUM_POINTS
  REAL(8) :: XMIN, XMAX, DX, X, RANGE
  REAL(8) :: SUM1,SUM2,SUMM,INTEGRAL(2)
  REAL(8) :: MAX_RANGE,MIN_RANGE,DIFF_RANGE
  REAL(8) :: sqrtpinv

  NUM_POINTS = 10000
  MAX_RANGE = 0.1D0
  DIFF_RANGE = 0.10D0
  MIN_RANGE = DIFF_RANGE
  SQRTPINV = 1.77245385090552D0

  !OPEN(UNIT=1,FILE='r_int1_int2_diff.dat')
  !WRITE(1,'(4(A20,1X))')'RANGE','SIMPSON.S INT','RANDOM SAMPLING',' DIFFERENCE'
  WRITE(*,'(4(A20,1X))')'RANGE','SIMPSON.S INT','RANDOM SAMPLING',' DIFFERENCE'

  RANGE = 10.0D0
!  RANGE = MIN_RANGE
!  RANGE_LOOP: DO WHILE (RANGE <= MAX_RANGE)

     XMIN =   -RANGE ! SET IT TO ZERO FOR INTEGRATION FROM ZERO TO RANGE
!     XMIN = 0.D0 ! SET IT TO ZERO FOR INTEGRATION FROM ZERO TO RANGE
     XMAX =  RANGE
     DX = (XMAX - XMIN) / DBLE(NUM_POINTS)

     !!!!!!!!!!!!!!!!!!!!!!
     ! SIMPSON'S 1/3 RULE !
     !!!!!!!!!!!!!!!!!!!!!!
     X = XMIN                      ! SIMPSON'S 1/3 RULE
     SUM1 = 0.0D0;  SUM2 = 0.0D0
     DO I = 1, NUM_POINTS - 1
        X = X + DX
        SELECT CASE (MOD(I,2))
        CASE(1)
           SUM1 = SUM1 + GAUSSIAN(X)
        CASE(0)
           SUM2 = SUM2 + GAUSSIAN(X)
        END SELECT
     END DO
     SUMM = GAUSSIAN(XMIN) + GAUSSIAN(XMAX) + 4.0D0 * SUM1 + 2.0D0 * SUM2
     SUMM = SUMM * DX / 3.0D0
     INTEGRAL(1) = SUMM/SQRTPINV

     !!!!!!!!!!!!!!!!!!!!!!!!!!!
     ! UNIFORM RANDOM SAMPLING !
     !!!!!!!!!!!!!!!!!!!!!!!!!!!
     SUMM = 0.0                 !UNIFORM RANDOM SAMPLING
     DO I = 0, NUM_POINTS
        CALL RANDOM_NUMBER(X)
        X = XMIN + (XMAX - XMIN) * X
        SUMM = SUMM + GAUSSIAN(X)
     END DO
     SUMM = DX * SUMM
     INTEGRAL(2) = SUMM/SQRTPINV
!     WRITE(1,'(4(F20.15,1X))')RANGE,INTEGRAL,ABS(INTEGRAL(2)-INTEGRAL(1))
!     RANGE = RANGE + DIFF_RANGE

!  END DO RANGE_LOOP

!     WRITE(1,'(4(F20.15,1X))')RANGE,INTEGRAL,ABS(INTEGRAL(2)-INTEGRAL(1))
     WRITE(*,'(4(F20.15,1X))')RANGE,INTEGRAL,ABS(INTEGRAL(2)-INTEGRAL(1))

CONTAINS
  
  REAL(8) FUNCTION GAUSSIAN(VAR)
    REAL(8), INTENT(IN) :: VAR
    REAL(8) :: TMP
    TMP = X * X
    GAUSSIAN = EXP(-TMP)
  END FUNCTION GAUSSIAN

END PROGRAM INTEGRATION