      PROGRAM GAUSSIT
C=======================================================================
C  INTEGRATION DEMONSTRATION MAKING USE OF GAUSS LEGENDRE METHOD    
C
C     EIJI FUKUMORI NOVEMBER 28, 1995 
C=======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER (MXINTEPT=6)
C=======================================================================
      DIMENSION SAI(MXINTEPT), W(MXINTEPT)
C=======================================================================
      CALL INPUT ( INTEPT, NSEG, XST, XEN )
C=======================================================================
      CALL GRULE ( MXINTEPT, INTEPT , SAI , W )
C=======================================================================
      CALL INTE ( MXINTEPT,INTEPT,SAI,W,NSEG,XST,XEN,AREA )
      WRITE (*,*) 'AREA = ', AREA
C=======================================================================
      OPEN ( 1, FILE='GAUSS.SIN', STATUS = 'UNKNOWN' )
	WRITE (1,*) 'INTEPT =', INTEPT
	WRITE (1,*) 'NSEG =' , NSEG
	WRITE (1,*) 'XST =', XST
	WRITE (1,*) 'XEN =',XEN
	WRITE (1,*)
	WRITE (1,*) 'GAUSS-LEGENDRE INTEGRATION COEFFICINTS(I,SAI(I),W(I))'
	WRITE (1,*) '   I            SAI(I)              W(I)'
	WRITE (1,*) '----------------------------------------'
	DO I = 1 , INTEPT
      WRITE (1,'(I5,2F18.14)') I,SAI(I),W(I)
	END DO
	WRITE (1,*)
      WRITE (1,*) 'AREA = ', AREA
      CLOSE (1)
      STOP
      END
C
C
      SUBROUTINE INTE ( MXINTEPT,INTEPT,SAI,W, NSEG,XST,XEN,AREA )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION SAI(MXINTEPT), W(MXINTEPT)
      DX = ( XEN - XST ) / NSEG
      H = DX / 2.
      AREA = 0.
      DO I = 1 , NSEG
      X1 = ( I - 1 ) * DX
      X2 = I * DX
      XM = ( X1 + X2 ) / 2. + XST
	CALL GAUSS ( MXINTEPT,INTEPT,SAI,W, H,XM, SUBAREA )
      AREA = AREA + SUBAREA
      END DO
      RETURN
      END
C
C
      SUBROUTINE GAUSS ( MXINTEPT,INTEPT,SAI,W, H,XM, SUBAREA )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION SAI(MXINTEPT), W(MXINTEPT)
      SUBAREA = 0.
      DO I = 1 , INTEPT
      X = H * SAI(I) + XM
      SUBAREA = SUBAREA + W(I)*F(X)
      END DO
	SUBAREA = SUBAREA * H
      RETURN
      END
C
C
      SUBROUTINE INPUT ( INTEPT, NSEG, XST, XEN )
      IMPLICIT REAL*8 ( A-H , O-Z )
      WRITE (*,210)
  210 FORMAT( 'NUMBER OF GAUSS POINTS = ' $ )
      READ (*,*) INTEPT
      IF ( INTEPT .LE. 0 ) INTEPT = 1
      IF ( INTEPT .GT. 6 ) INTEPT = 6
      WRITE (*,220)
  220 FORMAT( 'NUMBER OF SEGMENTS = '  $ )
      READ (*,*) NSEG
      IF ( NSEG .LE. 0 ) NSEG = 1
      WRITE (*,230)
  230 FORMAT( 'LOWER LIMIT OF INTEGRATION = '  $ )
      READ (*,*) XST
      WRITE (*,240)
  240 FORMAT( 'UPPER LIMIT OF INTEGRATION = '  $ )
      READ (*,*) XEN
      IF ( XST .GE. XEN ) THEN
      WRITE (*,*) 'XST >= XEN'
      STOP 'DATA DO NOT MAKE SENSE.'
	END IF
      RETURN
      END
C
C
      SUBROUTINE GRULE ( MXINTEPT, INTEPT, SAI , W )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION SAI(MXINTEPT) , W(MXINTEPT)
      IF ( INTEPT .LT. 2 ) STOP'INTEPT<2'
      IF ( INTEPT .GT. 6 ) STOP'INTEPT>6'
      GO TO ( 99, 20, 30, 40, 50, 60 ) , INTEPT
   99 STOP
   20 SAI(1) = DSQRT(3.D0)/3.D0
      W(1) = 1.D0
      GO TO 88
   30 SAI(1) = DSQRT(15.D0)/5.D0
      SAI(2) = 0.D0
      W(1) = 5.D0/ 9.D0
      W(2) = 8.D0/ 9.D0
      GO TO 88
   40 SAI(1) = 0.33998104358485D0
      SAI(2) = 0.86113631159405D0
        W(1) = 0.65214515486254D0
        W(2) = 0.34785484513745D0
      GO TO 88
   50 SAI(1) = 0.90617984593866D0
      SAI(2) = 0.53846931010568D0
      SAI(3) = 0.D0
        W(1) = 0.23692688505619D0
        W(2) = 0.47862867049937D0
        W(3) = 5.12D0 / 9.D0
      GO TO 88
   60 SAI(1) = 0.23861918608320D0
      SAI(2) = 0.66120938646626D0
      SAI(3) = 0.93246951420315D0
        W(1) = 0.46791393457269D0
        W(2) = 0.36076157304814D0
        W(3) = 0.17132449237917D0
   88 NN = INTEPT / 2
      DO I = 1 , NN
      J = INTEPT - I + 1
      SAI(J) = - SAI(I)
      W(J) = W(I)
      END DO
      RETURN
      END
C
C
      FUNCTION F (X)
      IMPLICIT REAL*8 ( A-H , O-Z )
      F = DSIN(X)
      RETURN
      END