REAL U,V,MUG,SUMG,DIF,SUMD,SDG,BETA    
      INTEGER I,N ,YN
      DIMENSION FB(5000),M(5000),B(5000),D(5000),G(5000) 
      CALL SEED(-1) 
      PRINT*,'======================================================='
      PRINT*,'|               Monte Carlo Simulation                |'
      PRINT*,'|              By  Pracha Vorarittinapa               |'
      PRINT*,'======================================================='
    1 PRINT*,'Enter no. of samples to be generated = ?'
      READ*,N 
      SUMG=0 
      DO 5, I=1,N
      CALL RANDOM(V)
      U=V
      CALL STAND(U,Z)
      FB(I)=EXP(7.329+0.31223*Z)
      CALL RANDOM(V)
      U=V
      CALL STAND(U,Z)
      B(I)=5.6+0.224*Z
      CALL RANDOM(V)
      U=V
      CALL STAND(U,Z)
      D(I)=11.4+0.342*Z
      CALL RANDOM(V)
      U=V
      M(I)=94599-(LOG(0-LOG(U)))/0.00011  
      G(I)=FB(I)-6*M(I)/(B(I)*D(I)*D(I)) 
      SUMG=SUMG+G(I)
    5 CONTINUE  
      MUG=SUMG/N 
      SUMD=0
      DO 10, I=1,N
      DIF=(G(I)-MUG)*(G(I)-MUG)
      SUMD=SUMD+DIF
   10 CONTINUE
      SDG=SQRT(SUMD/(N-1))
      BETA=MUG/SDG 
      PRINT*,'--------------------------------------'  
      PRINT*,'    Statistics of Safety Margin(G)'
      PRINT*,'--------------------------------------'
      PRINT*,'Mean               =',MUG
      PRINT*,'Standard Deviation =',SDG
      PRINT*,'Reliability Index  =',BETA    
      PRINT*,'--------------------------------------' 
      PRINT*,'New run?(0=yes,1=no)'
      READ*,YN
      IF (YN.EQ.0)THEN
      GOTO 1  
      ELSE
      PRINT*,'Quit'
      END IF
      END  
      SUBROUTINE STAND(U,Z)  
      PARAMETER(C0=2.515517,C1=0.802853,C2=0.010328)
      PARAMETER(D1=1.432788,D2=0.189269,D3=0.001308)
      IF(U.GT.0.5)THEN
      T=SQRT(0-LOG((1-U)*(1-U)))
      Z=T-(C0+C1*T+C2*T*T)/(1+D1*T+D2*T*T+D3*T*T*T)
      ELSE 
      T=SQRT(0-LOG(U*U))
      Z=(C0+C1*T+C2*T*T)/(1+D1*T+D2*T*T+D3*T*T*T)-T 
      END IF
      END