------
------
PROGRAM LOGIS1
C THE LOGISTIC MAP
C THIS PROGRAM GENERATES POINTS FOR THE LOGISTIC MAP
C X(N+1) = AX(1-X)
C WHERE A IS THE BIFURCATION PARAMETER. AT THE PROMPT, SELECT 1 FOR
C THE FIRST ITERATE, 2 FOR THE SECOND ITERATE, 3 FOR THE THIRD
C ITERATE, AND 4 FOR THE FOURTH ITERATE.
C THE VALUES OF X ARE TO BE FOUND IN A FILE TITLED 'RESULT.'
C
DIMENSION X(101),Y(101)
OPEN(2,FILE='RESULT',STATUS='NEW')
WRITE(*,10)
10 FORMAT(1X,'SELECT THE VALUE OF THE CONTROL PARAMETER. ',\)
READ(*,*) A
WRITE(*,20)
20 FORMAT(1X,'SELECT THE ITERATIVE MAP. ',\)
READ(*,*) N
X(1) = 0.0
GO TO (25,35,45,55) N
C FIRST ITERATE MAP
25 DO 30 I=1,100
X(I+1) = X(I)+0.01
30 Y(I) = A*X(I)*(1.0-X(I))
GO TO 65
C SECOND ITERATE MAP
35 DO 40 I=1,100
X(I+1) = X(I)+0.01
Y(I) = A*X(I)*(1.0-X(I))
40 Y(I) = A*Y(I)*(1.0-Y(I))
GO TO 65
C THIRD ITERATE MAP
45 DO 50 I=1,100
X(I+1) = X(I)+0.01
Y(I) = A*X(I)*(1.0-X(I))
Y(I) = A*Y(I)*(1.0-Y(I))
50 Y(I) = A*Y(I)*(1.0-Y(I))
GO TO 65
C FOURTH ITERATE MAP
55 DO 60 I=1,100
X(I+1) = X(I)+0.01
Y(I) = A*X(I)*(1.0-X(I))
Y(I) = A*Y(I)*(1.0-Y(I))
Y(I) = A*Y(I)*(1.0-Y(I))
60 Y(I) = A*Y(I)*(1.0-Y(I))
65 WRITE(2,70) (Y(I), I=1,101)
70 FORMAT(E14.5)
CLOSE(2)
WRITE(*,80)
80 FORMAT(1X,'PRESS RETURN TO TERMINATE.')
PAUSE
END
------
PROGRAM LOGIS2
C THE LOGISTIC TIME SERIES
C THIS PROGRAM GENERATES POINTS FOR THE LOGISTIC TIME SERIES X(N)
C VERSUS N FROM THE LOGISTIC EQUATION
C X(N+1) = AX(1-X)
C WHERE A IS THE BIFURCATION PARAMETER. THE NUMBER OF POINTS TO BE
C CALCUATED CANNOT BE MORE THAN 999.
C THE VALUES OF X ARE TO BE FOUND IN A FILE TITLED 'RESULT.'
C
DIMENSION X(1000)
OPEN(2,FILE='RESULT',STATUS='NEW')
WRITE(*,10)
10 FORMAT(1X,'SELECT THE NUMBER OF POINTS TO BE CALCULATED. ',\)
READ(*,*) N
WRITE(*,20)
20 FORMAT(1X,'SELECT THE VALUE OF THE INITIAL POINT. ',\)
READ(*,*) X(1)
WRITE(*,30)
30 FORMAT(1X,'SELECT THE VALUE OF THE CONTROL PARAMETER. ',\)
READ(*,*) A
DO 50 I=1,N+1
50 X(I+1) = A*X(I)*(1.0-X(I))
WRITE(2,60) (X(I), I=1,N+1)
60 FORMAT(E14.5)
CLOSE(2)
WRITE(*,70)
70 FORMAT(1X,'PRESS RETURN TO TERMINATE.')
PAUSE
END
------
PROGRAM LOTKA
C SOLUTION OF THE LOTKA-VOLTERRA EQUATIONS BY THE RUNGE-KUTTA METHOD
C THE EQUATIONS ARE
C DX/DT = X(A-BY)
C DY/DT = Y(BX-C)
C WHERE A,B,C ARE PARAMETERS. THE SOLUTIONS ARE ADVANCED BY A
C TIME INCREMENT DT. THE PRINTOUTS OF X AND Y ARE CONTROLLED BY
C THE NUMBER OF TIME INCREMENTS NUM BETWEEN WRITES TO THE OUTPUT
C FILE. THIS IS SET BY THE VALUE OF M, WHERE M = 1/NUM*DT.
C TOO LARGE A TIME INCREMENT WILL CAUSE ERRATIC BEHAVIOR OF THE
C PROGRAM. AN INITIAL VALUE OF 0.01 IS RECOMMENDED.
C WHEN ENTERING VALUES OF PARAMETERS UPON SEEING THE SCREEN PROMPT,
C YOU MUST LEAVE A SPACE BETWEEN NUMBERS.
C THE OUTPUT FILES ARE LABELED 'XRESULT' AND 'YRESULT.'
C
F1(X,Y,A,B) = X*(A-B*Y)
F2(X,Y,B,C) = Y*(B*X-C)
OPEN(2,FILE='XRESULT', STATUS='NEW')
OPEN(3,FILE='YRESULT', STATUS='NEW')
WRITE(*,10)
10 FORMAT(1X,'SELECT THE NUMBER OF POINTS TO BE CALCULATED, THE VALUE
XS OF A,B,C AND THE',/1X,'INITIAL VALUES OF X AND Y. ',\)
READ(*,*) N,A,B,C,X,Y
WRITE(*,25)
25 FORMAT(1X,'SELECT THE TIME INCREMENT AND THE VALUE OF M. ',\)
READ(*,*) DT, M
WRITE(2,40) X
WRITE(3,40) Y
C NUM IS A COUNTER FOR TIME INTERVALS BETWEEN OUTPUT WRITES
NUM = 0
DO 50 I=1,N
AK1 = F1(X,Y,A,B)*DT
AM1 = F2(X,Y,B,C)*DT
X1 = X+AK1/2.0
Y1 = Y+AM1/2.0
AK2 = F1(X1,Y1,A,B)*DT
AM2 = F2(X1,Y1,B,C)*DT
X2 = X+AK2/2.0
Y2 = Y+AM2/2.0
AK3 = F1(X2,Y2,A,B)*DT
AM3 = F2(X2,Y2,B,C)*DT
X3 = X+AK3
Y3 = Y+AM3
AK4 = F1(X3,Y3,A,B)*DT
AM4 = F2(X3,Y3,B,C)*DT
X1 = X+(AK1+2.0*AK2+2.0*AK3+AK4)/6.0
Y1 = Y+(AM1+2.0*AM2+2.0*AM3+AM4)/6.0
X = X1
Y = Y1
IF(NUM.EQ.0) GO TO 35
IF(NUM*M*DT.EQ.1) GO TO 30
GO TO 45
30 NUM = 0
35 WRITE(2,40) X
WRITE(3,40) Y
40 FORMAT(E14.5)
45 NUM = NUM + 1
50 CONTINUE
CLOSE(2)
CLOSE(3)
WRITE(*,60)
60 FORMAT(1X,'PRESS RETURN TO TERMINATE')
PAUSE
END
------
PROGRAM MCINT
C MONTE CARLO INTEGRATION
C THIS PROGRAM INTEGRATES 1/(1+X*X) USING THE RANDOM NUMBER
C GENERATOR RAN2(IDUM) FROM PRESS. THE PROGRAM ASKS
C FOR THE INPUT, WHICH IS THE NUMBER OF SAMPLE POINTS
C (RAMDOM NUMBERS) IN THE EVALUATION
C OF THE INTEGRAL.
C THE OUTPUT IS TO THE SCREEN.
C NOTE HOW MANY POINTS MUST BE SAMPLED TO APPROACH THE TRUE VALUE OF
C 2.65163.
C
IMPLICIT REAL*8(A-H,O-Z)
CHARACTER*3 ANS
3 WRITE (*,5)
5 FORMAT(1X,'ENTER THE NUMBER OF SAMPLE POINTS DESIRED. ',\)
READ (*,*) N
WRITE (*,7)
7 FORMAT(1X,'THE PROGRAM IS RUNNING.')
IDUM = -1
AINT = 0.0
DO 10 I=1,N
X = -4.0+8.0*RAN2(IDUM)
C THE FOLLOWING LINE IS THE FUNCTION, WHICH MAY BE CHANGED.
F = 1.0/(1.0+X*X)
AINT = AINT+F
10 CONTINUE
AINT = AINT*8.0/N
20 FORMAT(6E12.5)
WRITE (*,12) AINT
12 FORMAT(1X,'THE VALUE OF THE INTEGRAL IS',2X,E10.5)
WRITE (*,25)
25 FORMAT(1X,'DO YOU WISH TO DO ANOTHER CALCULATION? (yes OR no) ',\)
READ (*,30) ANS
30 FORMAT(A3)
IF (ANS.EQ.'yes') GO TO 3
WRITE (*,35)
35 FORMAT(1X,'PRESS RETURN TO TERMINATE')
PAUSE
END
------
PROGRAM POLYMER
C THIS FORTRAN PROGRAM ILLUSTRATES POLYMER FOLDING BY GENERATING A
C SELF-AVOIDING RANDOM WALK. THE POLYMER MODEL IS THE BEAD-SPRING
C MODEL. THE PROGRAM PROMPTS FOR THE NUMBER OF BEADS AND FOR THE
C SEED VALUE FOR THE RANDOM NUMBER GENERATOR.
C THE NUMBER OF BEADS CANNOT EXCEED 500.
C IN SOME CASES THE POLYMER CHAIN MUST BE TERMINATED BECAUSE IT
C CANNOT GROW FARTHER. SUCH CASES ARE INDICATED ON THE SCREEN AND
C THE POLYMER IS DRAWN AS FAR AS IT CAN GO.
C THIS IS NOT A STATISTICAL SELF-AVOIDING WALK, SINCE THE MOVES ARE
C NOT ALL EQUALLY PROBABLE. IT IS INTENDED ONLY TO SHOW POLYMER
C FOLDING.
C THIS PROGRAM REQUIRES A RANDOM NUMBER GENERATOR SUCH AS RAN2 FOUND
C IN PRESS.
C THIS PROGRAM WRITES THE OUTPUT TO FILES 'RESULTX' AND 'RESULTY.'
C THE OUTPUT FILES CAN BE EXPORTED TO A PLOTTING ROUTINE, OR THE
C RESULTS CAN BE PLOTTED ON GRAPH PAPER.
C
C
DIMENSION IX(500),IY(500)
OPEN(2,FILE='RESULTX',STATUS='NEW')
OPEN(3,FILE='RESULTY',STATUS='NEW')
C NUMBER OF BEADS IS N
WRITE(*,5)
5 FORMAT(1X,'PLEASE TYPE IN THE NUMBER OF BEADS. '\)
READ(*,*) N
WRITE(*,10)
10 FORMAT(1X,'PLEASE TYPE IN THE VALUE OF THE RANDOM NUMBER SEED. '\
X)
READ(*,*) IDUM
C SET UP INITIAL SITES
L = 1
IX(1) = 100
IY(1) = 100
IX(2) = 100+L
IY(2) = 100
C BEGIN THE RANDOM WALK
DO 150 I=3,N
100 A = RAN2(IDUM)
102 NCOUNT = 0
IF(A.LT.0.25) GO TO 115
IF(A.LT.0.5) GO TO 125
IF(A.LT.0.75) GO TO 135
105 NCOUNT = NCOUNT+1
IF(NCOUNT.GT.4) GO TO 500
DO 110 K=1,I-1
110 IF(IX(K).EQ.IX(I-1).AND.IY(K).EQ.IY(I-1)-L) GO TO 115
IX(I) = IX(I-1)
IY(I) = IY(I-1)-L
GO TO 150
115 NCOUNT = NCOUNT+1
IF(NCOUNT.GT.4) GO TO 500
DO 120 K=1,I-1
120 IF(IX(K).EQ.IX(I-1)-L.AND.IY(K).EQ.IY(I-1)) GO TO 125
IX(I) = IX(I-1)-L
IY(I) = IY(I-1)
GO TO 150
125 NCOUNT = NCOUNT+1
IF(NCOUNT.GT.4) GO TO 500
DO 130 K=1,I-1
130 IF(IX(K).EQ.IX(I-1).AND.IY(K).EQ.IY(I-1)+L) GO TO 135
IX(I) = IX(I-1)
IY(I) = IY(I-1)+L
GO TO 150
135 NCOUNT = NCOUNT+1
IF(NCOUNT.GT.4) GO TO 500
DO 140 K=1,I-1
140 IF(IX(K).EQ.IX(I-1)+L.AND.IY(K).EQ.IY(I-1)) GO TO 105
IX(I) = IX(I-1)+L
IY(I) = IY(I-1)
GO TO 150
C NCOUNT > 4; CHAIN IS STUCK. REEVALUATE PREVIOUS POSITION.
500 MCOUNT = 0
IF(A.LT.0.25) GO TO 510
IF(A.LT.0.5) GO TO 520
IF(A.LT.0.75) GO TO 530
540 MCOUNT = MCOUNT+1
IF(MCOUNT.GT.4) GO TO 590
DO 505 K=1,I-1
505 IF(IX(K).EQ.IX(I-2).AND.IY(K).EQ.IY(I-2)-L) GO TO 510
GO TO 550
510 MCOUNT = MCOUNT+1
IF(MCOUNT.GT.4) GO TO 590
DO 515 K=1,I-1
515 IF(IX(K).EQ.IX(I-2)-L.AND.IY(K).EQ.IY(I-2)) GO TO 520
GO TO 560
520 MCOUNT = MCOUNT+1
IF(MCOUNT.GT.4) GO TO 590
DO 525 K=1,I-1
525 IF(IX(K).EQ.IX(I-2).AND.IY(K).EQ.IY(I-2)+L) GO TO 530
GO TO 570
530 MCOUNT = MCOUNT+1
IF(MCOUNT.GT.4) GO TO 590
DO 535 K=1,I-1
535 IF(IX(K).EQ.IX(I-2)+L.AND.IY(K).EQ.IY(I-2)) GO TO 540
GO TO 580
550 IX(I-1) = IX(I-2)
IY(I-1) = IY(I-2)-L
GO TO 102
560 IX(I-1) = IX(I-2)-L
IY(I-1) = IY(I-2)
GO TO 102
570 IX(I-1) = IX(I-2)
IY(I-1) = IY(I-2)+L
GO TO 102
580 IX(I-1) = IX(I-2)+L
IY(I-1) = IY(I-2)
GO TO 102
590 WRITE(*,595)
595 FORMAT('CHAIN GROWTH TERMINATED')
GO TO 205
150 CONTINUE
C WRITE FINAL CONFIGURATION TO FILES "RESULTX" AND "RESULTY"
205 WRITE(2,645) (IX(I),I=1,N)
WRITE(3,645) (IY(I),I=1,N)
645 FORMAT(I5)
CLOSE(2)
CLOSE(3)
WRITE(*,650)
650 FORMAT(1X,'PRESS RETURN TO TERMINATE')
PAUSE
END
------
PROGRAM PWRSPR
C POWER SPECTRUM
C THE INPUT FOR THE POWER SPECTRUM CONSISTS OF A TIME SERIES,
C EITHER CALCULATED OR MEASURED IN AN EXPERIMENT. THESE DATA
C ARE IN A FILE NAMED 'DATA.' THE NUMBER OF DATA POINTS
C OF THE TIME SERIES MUST BE PROVIDED ACCORDING TO THE PROMPT.
C THE OUTPUT IS THE POWER SPECTRUM IN A FILE NAMED 'RESULT.'
C THERE IS NO WINDOW PROVIDED FOR THE FOURIER TRANSFORM, SO THERE
C MAY BE LEAKAGE IN THE TRANSFORM AND THE POWER SPECTRUM.
C ALSO, CONSTANT COEFFICIENTS SUCH AS THE TIME INCREMENT AND THE
C TOTAL TIME HAVE BEEN OMITTED.
C THE OUTPUT DATA MAY BE EXPORTED TO YOUR FAVORITE GRAPHICS PROGRAM.
C
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION X(2700),Y(2700),YIMAG(2700)
DATA PI/3.141593/
OPEN(2,FILE='DATA',STATUS='OLD',FORM='FORMATTED')
OPEN(3,FILE='RESULT',STATUS='NEW',FORM='FORMATTED')
WRITE(*,10)
10 FORMAT(1X,'SELECT THE NUMBER OF POINTS FOR THE POWER SPECTRUM. TH
XIS SHOULD AGREE WITH',/1X,'THE OUTPUT (LESS ONE) FROM THE TIME SER
XIES. ',\)
READ(*,*) N
READ (2,20) (X(I), I=1,N)
20 FORMAT (E14.5)
C FOURIER TRANSFORM
DO 50 J=1,N
Y(J) = 0.0
YIMAG(J) = 0.0
DO 40 I=1,N
Y(J) = X(I)*COS(2.0*PI*J*(I-1)/(N+1))+Y(J)
40 YIMAG(J) = X(I)*SIN(2.0*PI*J*(I-1)/(N+1))+YIMAG(J)
50 CONTINUE
DO 60 J=1,N
Y(J) = Y(J)*Y(J)
YIMAG(J) = YIMAG(J)*YIMAG(J)
60 Y(J) = Y(J)+YIMAG(J)
WRITE(3,20) (Y(J), J=1,N)
CLOSE(2)
CLOSE(3)
WRITE(*,70)
70 FORMAT(1X,'PRESS RETURN TO TERMINATE')
PAUSE
END
------
PROGRAM REVCA
C REVERSIBLE CELLULAR AUTOMATA
C THIS PROGRAM CALCULATES ONE-DIMENSIONAL REVERSIBLE
C CELLULAR AUTOMATA USING THE NEIGHBOR-THREE RULES OF
C [WOLFRAM, 1983]. THE INPUT CONSISTS OF THE NUMBER OF
C SITES (CELLS) N IN THE CA PLUS TWO TO TAKE INTO ACCOUNT THE ENDS
C FOR CYCLIC BOUNDARY CONDITIONS, THE NUMBER OF TIME VALUES
C (INCLUDING TIME ZERO) NT, AND THE ORIGINAL CONFIGURATION IX(NT,N)
C COMPRISED OF 1 OR 0. THE VALUES OF N, NT, AND THE INITIAL
C CONFIGURATION MUST BE ENTERED IN A FILE LABELED 'INPUT.'
C SEE THE PROGRAM FOR THE FORMAT OF THE INPUT DATA.
C THE CA RULES IN THE PROGRAM MUST BE CHANGED FOR EACH NEW CA AND
C THE PROGRAM MUST BE RECOMPILED FOR EACH NEW CA.
C THE VALUE OF N MUST NOT EXCEED 100 AND THE VALUE OF NT MUST NOT
C EXCEED 500.
C THE CA APPEAR DIRECTLY ON THE SCREEN.
C
IMPLICIT REAL*8(A-H,O-V), INTEGER(I-N)
DIMENSION IX(500,100),IY(500,100)
OPEN(2, FILE='INPUT',STATUS='OLD')
C READ IN THE NUMBER OF CELLS (+2) AND THE NUMBER OF TIME STEPS
READ(2,*) N,NT
C READ IN THE INITIAL CONFIGURATION.
READ(2,*) (IX(1,J), J=1,N)
DO 20 J=2,N-1
C*****THE CA RULES MUST BE SPECIFIED BELOW. THE NEIGHBORHOODS ARE
C*****LABELED. IF DESIRED, THE CA NUMBER CAN BE RECORDED IN THE
C*****NEXT STATEMENT.
C CELLULAR AUTOMATA RULES FOR CA 126
IF(IX(1,J).EQ.1) GO TO 15
C NEIGHBORHOOD 101
IF(IX(1,J-1).EQ.1.AND.IX(1,J+1).EQ.1) IX(1+1,J) = 1
C NEIGHBORHOOD 100
IF(IX(1,J-1).EQ.1.AND.IX(1,J+1).EQ.0) IX(1+1,J) = 1
C NEIGHBORHOOD 001
IF(IX(1,J-1).EQ.0.AND.IX(1,J+1).EQ.1) IX(1+1,J) = 1
C NEIGHBORHOOD 000
IF(IX(1,J-1).EQ.0.AND.IX(1,J+1).EQ.0) IX(1+1,J) = 0
GO TO 20
C NEIGHBORHOOD 111
15 IF(IX(1,J-1).EQ.1.AND.IX(1,J+1).EQ.1) IX(1+1,J) = 0
C NEIGHBORHOOD 110
IF(IX(1,J-1).EQ.1.AND.IX(1,J+1).EQ.0) IX(1+1,J) = 1
C NEIGHBORHOOD 011
IF(IX(1,J-1).EQ.0.AND.IX(1,J+1).EQ.1) IX(1+1,J) = 1
C NEIGHBORHOOD 010
IF(IX(1,J-1).EQ.0.AND.IX(1,J+1).EQ.0) IX(1+1,J) = 1
20 CONTINUE
IX(1+1,1) = IX(1+1,N-1)
IX(1+1,N) = IX(1+1,2)
DO 120 I=2,NT-1
DO 118 J=2,N-1
C*****THE CA RULES MUST BE SPECIFIED BELOW. THE NEIGHBORHOODS ARE
C*****LABELED. IF DESIRED, THE CA NUMBER CAN BE RECORDED IN THE
C*****NEXT STATEMENT.
C CELLULAR AUTOMATA RULES FOR CA 126
IF(IX(I,J).EQ.1) GO TO 108
C NEIGHBORHOOD 101
IF(IX(I,J-1).EQ.1.AND.IX(I,J+1).EQ.1) IY(I+1,J) = 1
C NEIGHBORHOOD 100
IF(IX(I,J-1).EQ.1.AND.IX(I,J+1).EQ.0) IY(I+1,J) = 1
C NEIGHBORHOOD 001
IF(IX(I,J-1).EQ.0.AND.IX(I,J+1).EQ.1) IY(I+1,J) = 1
C NEIGHBORHOOD 000
IF(IX(I,J-1).EQ.0.AND.IX(I,J+1).EQ.0) IY(I+1,J) = 0
C ADD LINE NT-2 MOD 2 (DO NOT CHANGE THE NEXT FOUR STATEMENTS.)
IF(IX(I-1,J).EQ.0.AND.IY(I+1,J).EQ.0) IX(I+1,J) = 0
IF(IX(I-1,J).EQ.0.AND.IY(I+1,J).EQ.1) IX(I+1,J) = 1
IF(IX(I-1,J).EQ.1.AND.IY(I+1,J).EQ.0) IX(I+1,J) = 1
IF(IX(I-1,J).EQ.1.AND.IY(I+1,J).EQ.1) IX(I+1,J) = 0
GO TO 118
C NEIGHBORHOOD 111
108 IF(IX(I,J-1).EQ.1.AND.IX(I,J+1).EQ.1) IY(I+1,J) = 0
C NEIGHBORHOOD 110
IF(IX(I,J-1).EQ.1.AND.IX(I,J+1).EQ.0) IY(I+1,J) = 1
C NEIGHBORHOOD 011
IF(IX(I,J-1).EQ.0.AND.IX(I,J+1).EQ.1) IY(I+1,J) = 1
C NEIGHBORHOOD 010
IF(IX(I,J-1).EQ.0.AND.IX(I,J+1).EQ.0) IY(I+1,J) = 1
C ADD LINE NT-2 MOD 2 (DO NOT CHANGE THE NEXT FOUR STATEMENTS.)
IF(IX(I-1,J).EQ.0.AND.IY(I+1,J).EQ.0) IX(I+1,J) = 0
IF(IX(I-1,J).EQ.0.AND.IY(I+1,J).EQ.1) IX(I+1,J) = 1
IF(IX(I-1,J).EQ.1.AND.IY(I+1,J).EQ.0) IX(I+1,J) = 1
IF(IX(I-1,J).EQ.1.AND.IY(I+1,J).EQ.1) IX(I+1,J) = 0
118 CONTINUE
C ADD CYCLIC BOUNDARY CONDITIONS
IX(I+1,1) = IX(I+1,N-1)
IX(I+1,N) = IX(I+1,2)
120 CONTINUE
DO 125 I=1,NT
DO 127 J=1,N
IF(IX(I,J).EQ.0) GO TO 128
WRITE(*,130)
130 FORMAT(''\)
GOTO 127
128 WRITE(*,129)
129 FORMAT(' '\)
127 CONTINUE
WRITE(*,126)
126 FORMAT(' ')
125 CONTINUE
CLOSE(2)
C SCREEN REMARK
WRITE(*,160)
160 FORMAT(1X,'PRESS RETURN TO TERMINATE')
PAUSE
END
------
PROGRAM SQRCA
C SQUARE CELLULAR AUTOMATON
C THIS PROGRAM FOLLOWS THE EVOLUTION OF FIVE-NEIGHBOR, SUM MODULO 2
C RULE FOR CA 614. CYCLIC BOUNDARY CONDITIONS ARE USED. THE (1,1),
C (1,N), (N,1), AND (N,N) ELEMENTS ARE NEVER USED AND CAN BE READ IN
C AS 0.
C THE INPUT IS IN A FILE LABELED 'INPUT' AND CONTAINS THE SIZE OF
C THE CA, THE NUMBER OF TIME STEPS, AND THE INITIAL CONFIGURATION.
C SEE THE PROGRAM FOR THE FORMAT OF THE INPUT DATA.
C THE OUTPUT IS DIRECTLY TO THE SCREEN.
C
IMPLICIT REAL*8(A-H,O-V), INTEGER(I-N)
DIMENSION IX(100,100),IY(100,100)
OPEN(2, FILE='INPUT',STATUS='OLD')
C READ IN THE SIZE OF THE CA (+2) AND THE NUMBER OF TIME STEPS (NT)
READ(2,*) N,NT
C READ IN THE INITIAL CONFIGURATION
C THE ARRAY IS READ IN WITH THE COLUMN INDEX INCREASING MOST RAPIDLY
READ(2,*) ((IX(I,J), J=1,N), I=1,N)
DO 100 IT=1,NT
DO 30 I=2,N-1
DO 20 J=2,N-1
IY(I,J) = IX(I,J-1)+IX(I,J)+IX(I,J+1)+IX(I-1,J)+IX(I+1,J)
20 IY(I,J) = MOD(IY(I,J),2)
30 CONTINUE
C SET UP NEW CONFIGURATION
DO 60 I=2,N-1
DO 50 J=2,N-1
50 IX(I,J) = IY(I,J)
60 CONTINUE
C IMPOSE CYCLIC BOUNDARY CONDITIONS
DO 80 I=2,N-1
IX(I,1) = IX(I,N-1)
80 IX(I,N) = IX(I,2)
DO 90 J=2,N-1
IX(1,J) = IX(N-1,J)
90 IX(N,J) = IX(2,J)
100 CONTINUE
C PRINT OUT IN A SQUARE ARRAY
DO 140 I=2,N-1
DO 110 J=2,N-1
IF(IX(I,J).EQ.0) GO TO 120
WRITE(*,115)
115 FORMAT(''\)
GOTO 110
120 WRITE(*,125)
125 FORMAT(' '\)
110 CONTINUE
WRITE(*,130)
130 FORMAT(' ')
140 CONTINUE
CLOSE(2)
C SCREEN REMARK
WRITE(*,160)
160 FORMAT(1X,'PRESS RETURN TO TERMINATE')
PAUSE
END