PROGRAM NumInt
**********************************************************************
*     Berechnung eines eindimensionalen Integrals mittels
*     Trapezverfahrens und Monte Carlo Technik 
*********************************************************************
 
      REAL R(97)     
      EXTERNAL RAN

      F(X) = 6.*X**2-10.*X+12.
      STF(X) = 2.*X**3-5*X**2+12.*X

      PRINT *, 
     &  'Numerische Berechnung der Integrals ueber die Funktion'
      PRINT *, 
     &  'F(X) = 6.*X**2-10.*X+12. von A=2 bis B=12'



      PRINT *, 
     &  'Eingabe Anzahl der Stützpunkte und ISEED fuer den',
     &  ' Zufallszahlengenerator:'
      READ *, M,ISEED

      A=2
      B=12
      RINT=0.     


C     RAN Zufallszahl zwischen 0 und 1
      DO 100, I=1,M
          X=A+RAN(ISEED,IX1,IX2,IX3,R)*(B-A)
          RINT=RINT+F(X)
 100  CONTINUE

      RINT=(B-A)*RINT/REAL(M)

C     Trapezregel

      TINT = 0. 
      TINT = TINT + F(A) + F(B)
     
      IF (M.GT.0) THEN
        DO 200, I=1,M-1
          X=A+ REAL(I)*(B-A)/REAL(M)
          TINT=TINT+2.*F(X)
 200    CONTINUE
        TINT = 0.5*TINT*(B-A)/REAL(M)


        EXAKT=STF(B)-STF(A)   

        PRINT *, 'Wert des Integrals nach der Trapezregel und ',
     &    'dem Monte Carlo Verfahren:'
        PRINT *,    TINT,RINT

      ELSE  
        PRINT *, 'Anzahl der Stuetzpunkte muss groesser als Null sein'         
      ENDIF      

      END       


*
*
      REAL FUNCTION RAN(IDUM,IX1,IX2,IX3,R)
**********************************************************************
*     Returns a random deviate between 0.0 and 1.0. Set IDUM to
*     any negative value to initialize or reinitialize the sequence.
*     Numerical Recipes p.196
**********************************************************************
      REAL R(97)
      INTEGER  IX1,IX2,IX3,IC1,IC2,IC3,IA1,IA2,IA3,IFF
      PARAMETER (M1=259200,IA1=7141,IC1=54733,RM1=1./M1)
      PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=1./M2)
      PARAMETER (M3=243000,IA3=4561,IC3=51349)

      DATA IFF /0/

      IF (IDUM.LT.0.OR.IFF.EQ.0) THEN
         IFF=1
         IX1=MOD(IC1-IDUM,M1)
         IX1=MOD(IA1*IX1+IC1,M1)
         IX2=MOD(IX1,M2)
         IX1=MOD(IA1*IX1+IC1,M1)
         IX3=MOD(IX1,M3)
         DO 11, J=1,97
           IX1=MOD(IA1*IX1+IC1,M1)
           IX2=MOD(IA2*IX2+IC2,M2)
           R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
 11      CONTINUE
         IDUM=1

       ENDIF

       IX1=MOD(IA1*IX1+IC1,M1)
       IX2=MOD(IA2*IX2+IC2,M2)
       IX3=MOD(IA3*IX3+IC3,M3)
       
       J=1+(97*IX3)/M3
       
*       IF (J.GT.97.OR.J.LT.1) PAUSE

       RAN=R(J)  


       R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1

       RETURN
       END
*
*