PROGRAM APREQI
C
C     APPROACH TO EQUILIBRIUM
C     -----------------------------------------------
C     AUTHOR: ROBERT EHRLICH (GEORGE MASON U    1973)
C     Revised: Peter Signell (Mich. St. U.      3/86)
C     Revised: Gary Herzenstiel (Mich. St. U.   3/86)
C     Revised: E.J.D.Kales (Mich. St. U.       11/87)
C     Revised: Peter Signell (Mich. St. U.      3/02)
C
C     This program calculates and plots the fraction of
C     coins in a set that are tails after specified number
C     of random tosses are made. It is assumed that all
C     coins are initially heads.
C
C     SUBROUTINES:
C          GETINFO - get values for R,L,C,DT,N,V0,W
C          CALC    - flip the coins and put in graph line
C          DISPLA  - display the graph 
C          QUERY   - asks if user wants to try again
C
C                            A P R 2 E Q I
C                                  :
C            ----------------------------------------------
C            :                :             :             :
C         GETINFO           CALC          DISPLAY      QUERY
C      (COINS,TURNS) (COINS,TURNS,    (NS,LINE,ARRAY) (AGAIN)
C                    NS,LINE,ARRAY)                  
C   
C     ***** DECLARATIONS *****
      REAL COINS,TURNS,ARRAY(501)
      INTEGER NS, IX
      CHARACTER LINE(51)*81,AGAIN
      OPEN(UNIT=3,FILE='(T:81)m357p1f.out')
      IX=13
C                    
C     ***** M A I N   B O D Y *****
C
      WRITE(3,'(/)')
 10   CALL GETINF(                 COINS,TURNS)
      CALL CALC   (COINS,TURNS,IX,   NS,LINE,ARRAY,NLINES)
      CALL DISPLA (NS,LINE,ARRAY,NLINES            )  
      CALL QUERY  (                       AGAIN)
      IF(AGAIN.EQ.'Y') THEN
          WRITE(3,'(A1)') 12
          GOTO 10
      ENDIF
      IF(AGAIN.EQ.'y') THEN
          WRITE(3,'(A1)') 12
          GOTO 10
      ENDIF
      CLOSE(UNIT=3)
      STOP
      END
C
C     ***** G E T I N F O *****
C
      SUBROUTINE GETINF(    COINS,TURNS)
C 
C     get values for COINS and TURNS
C
C     ***** DECLARATIONS *****
      REAL COINS,TURNS
C     
C     COINS = number of coins
C     TURNS = number of random turns
C
      PRINT *,'Please supply COINS & TURNS values:'
      WRITE(*,20) 'COINS='
      READ *,COINS
      WRITE(*,20) 'TURNS='
      READ *,TURNS
  20  FORMAT(/A6)
      PRINT *
      WRITE(*,24) COINS,TURNS
      WRITE(3,24) COINS,TURNS
  24  FORMAT('m357p1.for: THE APPROACH TO EQUILIBRIUM'/
     +'     Coins=',F7.0/'     Turns=',F7.0/)
      IF (COINS.GE.1001) NCOINS = 1001 - 1
      RETURN                                         
      END
C
C     ***** C A L C *****
C
      SUBROUTINE CALC(COINS,TURNS,IX,       NS,LINE,ARRAY,NLINES)
C 
C     turns the coins and makes the lines of graph output
C
C     ***** DECLARATIONS *****
      REAL COINS,TURNS,RAN,F,ARRAY(51)
      INTEGER NC,NT,NS,TAILS,K,J,NL,L,IX,NLINES
      CHARACTER LINE(51)*81,U*1001
C
C     LINE  = a line of the output graph
C     U     = H for heads or T for tails
C
      TAILS=0
      NC=COINS
      NT=TURNS
      NS=(NT+49)/50
      NLINES=NT/NS
C
C     initialize the first output line (turn 0)
      NL=1
      LINE(NL)(1:64)=' '
      LINE(NL)(1:1)='*'
      LINE(NL)(51:51)='|'
      LINE(NL)(64:64)='.'
      ARRAY(NL)=0.000
C
C     Initially assume all coins are Heads
C
      DO 30,K=1,NC
           U(K:K)='H'
 30   CONTINUE
C
C     loop over the number of turns
      DO 40,K=1,NT
C
C          generate a random number R (within the dashed
C          lines)
C
C          ---------------------------------------------
           IX=899*IX
           IX=IX-32768*INT(IX/32768)
           RAN=IX
           RAN=RAN/32767.0
C          ---------------------------------------------
C
C          J is a random integer from 1 to NC
C
           J=RAN*(COINS-1.0)+1.0
           IF(J.LT.1)J=1
           IF(J.GT.1001)J=1001
C
C          turn the Jth coin over and revise the number of
C          tails
C
           IF(U(J:J).EQ.'H')THEN
                U(J:J)='T'
                TAILS=TAILS+1
           ELSE
                U(J:J)='H'
                TAILS=TAILS-1
           ENDIF
           F=TAILS/COINS
C
C          only display results if K is a multiple of NS
C
           IF(K.EQ.NS*(K/NS))THEN
                NL=NL+1
C               initialize the NL-th output line to blanks
                LINE(NL)(1:64)=' '
                LINE(NL)(1:1)='.'
                LINE(NL)(51:51)='|'
                LINE(NL)(64:64)='.'
                L=INT(100.0*F+1)
                LINE(NL)(L:L)='*'
                ARRAY(NL)=F
           ENDIF
C
 40   CONTINUE
      RETURN
      END
C
C     ***** D I S P L A Y *****
C
      SUBROUTINE DISPLA(NS,LINE,ARRAY,NLINES    )
C 
C     display the graphical output from the program
C
C     ***** DECLARATIONS *****
      REAL ARRAY(51)
      INTEGER NS,X
      CHARACTER LINE(51)*81
C
C     write the table header
      WRITE(3,52)
 52   FORMAT('       Fraction',49x,'50%')
      WRITE(3,55) '  Turn  Tails  '
 55   FORMAT(A15,30H..............................,
     +34H..................................)

C     write the table lines
      WRITE(3,56) ((X-1)*NS,ARRAY(X),LINE(X)(1:64),X=1,NLINES+1)
 56   FORMAT(I6,F7.3,2X,A64)

C     write the table footer
      WRITE(3,55) ' '
      RETURN
      END
C
C     ***** Q U E R Y *****
C
      SUBROUTINE QUERY(    AGAIN)
C 
C     ***** DECLARATIONS *****
      CHARACTER AGAIN
C
      PRINT *
      WRITE(*,62)
 62   FORMAT(
     +'Would you like to do this again using new values? (Y/N)? '/)
      READ( *,'(A1)')AGAIN
      RETURN
      END
C     ***** P R O G R A M    E N D *****