'******************************************************************************
'Subject: 1-D continuous-valued cellular automaton:
' the Fermi-Pasta-Ulam quadratic nonlinear wave model.
'Author : Sjoerd.J.Schaper
'Date : 02-26-2006
'Code : all QBasic's, FreeBasic extendable
'Keys : [Esc] quit program
' [Enter] change palette
' [+] and [-] adjust wave amplitude
' [Tab] graph or space-time diagram
' [Backspc] reverse the simulation
' [s] scroll/ lock display
' [l] linear/ nonlinear rule
' [Up] and [Down arrow] tune non-linearity
' [Space] freeze
' press any other key to reset
'******************************************************************************
'This program is copyright (c) 2006 by the author. It is made available as is,
'and no warranty - about the program, its performance, or its conformity to any
'specification - is given or implied. It may be used, modified, and distributed
'freely, as long as the original author is credited as part of the final work.
'******************************************************************************
REM $STATIC
DEFINT A-Z
DECLARE SUB populate (t&)
'display the automaton at time t
DECLARE SUB evaluate (fl)
'apply the quadratic wave rule
DECLARE FUNCTION mix ()
'mix 8-bit colours
DECLARE FUNCTION Fsum! ()
'seed with Fourier sine waves sum,
'return non-linearity
DECLARE FUNCTION lcm& (x, y)
CONST MODE = 13, dW = 320, dH = 200
' set screenmode, width, height
'CONST MODE = 18, dW = 640, dH = 480 ' FreeBasic only
CONST MAX = 198
' max. grid size 198 for Qbasic's
CONST WAVE = .6666667
CONST NOLI = MAX * .002
CONST GAIN = 1.41, DNOL = 1.035
' parameter defaults
r = 1 + dW * .25 * .618 ' viewport size
s = r + dW * .75 - 1 ' FreeBasic workaround:
r = s - r + 1: s& = dH - 2 ' declaration of local array img()
DIM SHARED img(r * s& \ 2) ' isn't allowed inside Sub populate.
DIM SHARED c(MAX, 1) AS SINGLE, sw, plot, roll
' cell states, generation/ display mode switches
DIM SHARED Wav AS SINGLE, noL(2) AS SINGLE
' wave and non-linearity parameters
RANDOMIZE TIMER
DIM f AS SINGLE, t AS LONG
fl = 0: f = 1
dt = 1: t = 0
DO
populate t
t = t + dt
evaluate fl
IF t = 0 AND dt = -1 THEN
LOCATE 7, 1: PRINT "key..."
SLEEP
LOCATE 7, 1: PRINT " "
WHILE INKEY$ <> "": WEND
sw = 1 - sw: dt = 1: t = 1
END IF
g$ = INKEY$
IF g$ <> "" THEN
sc = ASC(RIGHT$(g$, 1))
IF LEN(g$) > 1 THEN sc = sc + 128
SELECT CASE sc
CASE 8
sw = 1 - sw
dt = -dt: t = t + dt ' reverse time
CASE 9
plot = NOT plot ' display toggle
CASE 13
IF NOT plot THEN COLOR mix
CASE 27
EXIT DO ' quit
CASE 32
LOCATE 4, 1: PRINT LTRIM$(STR$(t));
SLEEP ' freeze
CASE 43, 45
a! = GAIN
IF sc = 45 THEN a! = 1 / a!
FOR k = 0 TO MAX - 1 ' gain
c(k, 0) = c(k, 0) * a!
c(k, 1) = c(k, 1) * a!
NEXT k
CASE 108
fl = 1 - fl: g$ = "lin."
IF fl THEN g$ = "noLi" ' rule toggle
LOCATE 7, 1: PRINT "mode";
LOCATE 8, 1: PRINT g$;
CASE 115
roll = NOT roll ' scroll toggle
CASE 200, 208
a! = DNOL
IF sc = 208 THEN a! = 1 / a!
f = f * a!
IF f > 10 THEN f = 10
noL(1) = noL(2) * f ' tuning
CASE ELSE
noL(2) = Fsum ' reset
noL(1) = noL(2) * f
dt = 1: t = 0
END SELECT
WHILE INKEY$ <> "": WEND
END IF
LOOP
SYSTEM
'CGA colour components:
'blue
DATA 60,00,60,00,60,00,25,10
DATA 60,25,60,25,60,25,60,60
'green
DATA 00,60,60,00,00,60,25,10
DATA 25,60,60,25,25,60,60,00
'red
DATA 00,00,00,60,60,60,25,10
DATA 25,25,25,60,60,60,60,00
SUB evaluate (fl)
DIM dl AS SINGLE, dr AS SINGLE, s AS SINGLE
DIM u AS SINGLE, ul AS SINGLE, ur AS SINGLE
ul = c(MAX - 1, sw) ' wrap boundary
c(MAX, sw) = c(0, sw)
FOR k = 0 TO MAX - 1
u = c(k, sw)
ur = c(k + 1, sw)
dl = u - ul: dr = ur - u
s = ul + ur + (dr * dr - dl * dl) * noL(fl)
ul = u: u = u + u
s = (s - u) * Wav + u - c(k, 1 - sw)
IF ABS(s) > 2 THEN s = SGN(s) * 2 ' clip
c(k, 1 - sw) = s
NEXT k
sw = 1 - sw ' switch generations
END SUB
FUNCTION Fsum!
h = 8: DIM d AS SINGLE, i(h) ' 3 octave harmonics
DIM s AS SINGLE, w AS SINGLE
FOR r = 1 TO h
i(r) = r
NEXT r
FOR r = 1 TO h ' shuffle
t = 1 + INT(RND * h)
SWAP i(r), i(t)
NEXT r
d = 2 * ATN(1)
w = 4 * d / MAX ' velocity
d = d * INT(RND * 2) ' phase shift
r1 = RND < .67: r2 = RND < .5
FOR k = 0 TO MAX - 1
s = SIN(d + k * w * i(1)) / i(1)
IF r1 THEN s = s + SIN(d + k * w * i(2)) / i(2)
IF r2 THEN s = s + SIN(d + k * w * i(3)) / i(3)
c(k, 0) = s * 1.09: c(k, 1) = c(k, 0)
NEXT k
r = i(1)
IF r1 THEN r = lcm(r, i(2))
IF r2 THEN r = lcm(r, i(3))
Fsum = NOLI / SQR(r) ' non-linearity
END FUNCTION
FUNCTION lcm& (x, y)
a& = x: b = y
DO
r = x MOD y
x = y: y = r
LOOP UNTIL r = 0
lcm = a& * (b \ x) ' least common multiple
END FUNCTION
FUNCTION mix
DIM i(15), col(255) AS LONG
bands = 4
a! = 253 / (15 * bands)
FOR s = 0 TO 15 ' indices
i(s) = 2 + s * a!
NEXT s
col(0) = &H40406
col(1) = &H242426
FOR r = 0 TO 2
i0 = i(0)
READ c0: c0 = c0 * RND
FOR s = 1 TO 15
READ c1: c1 = c1 * RND
a! = (c1 - c0) / (i(s) - i0)
FOR t = i0 TO i(s) ' interpolate
u = c0 + (t - i0) * a!
col(t) = col(t) * 256 + u
NEXT t
i0 = i(s) + 1: c0 = c1
NEXT s
NEXT r
FOR t = i0 TO 255 ' duplicate
col(t) = col(2 + t - i0)
NEXT t
PALETTE USING col(0)
RESTORE
mix = 1
END FUNCTION
SUB populate (gen&) STATIC
IF (gen& = 0) OR (pl XOR plot) THEN
IF NOT fl THEN
SCREEN MODE, 8 ' 8 bpp
COLOR mix
r = 1 + dW * .25 * .618 ' viewport
s = r + dW * .75 - 1: fl = -1
VIEW (r, 1)-(s, dH - 2), , 1
a! = 253 / 4: d! = .491
LOCATE 3, 1: PRINT "genr.";
m = MAX / SQR(WAVE) ' initialize
Wav = MAX / m: Wav = Wav * Wav
noL(2) = Fsum
noL(1) = noL(2)
noL(0) = 0: sw = 0
plot = 0: roll = 0
END IF
pl = plot: li = 0
LOCATE 4, 1: PRINT " ";
IF plot THEN
WINDOW SCREEN (0, -2)-(MAX - 1, 2)
LINE (0, -2)-(MAX - 1, 2), 0, BF
ELSE
WINDOW SCREEN (-1, -1)-(MAX, m)
IF r THEN LINE (-1, -1)-(MAX, m), 0, BF
END IF
END IF
IF plot THEN
WAIT &H3DA, 8 ' slow down...
FOR k = 1 TO MAX - 1 ' plot curve
t = k - 1
LINE (t, c(t, 1 - sw))-(k, c(k, 1 - sw)), 0
LINE (t, c(t, sw))-(k, c(k, sw)), 1
NEXT k
ELSE
IF roll THEN
IF li < m - 1 THEN
li = li + 1
ELSE
WAIT &H3DA, 8 ' scroll up
GET (-d!, 0)-(MAX - d!, m - d!), img
PUT (-d!, -1), img, PSET
END IF
ELSE
li = (li + 1) MOD m ' cycle screen
END IF
FOR k = 0 TO MAX - 1 ' update CA
col = 2 + (c(k, sw) + 2) * a!
LINE (k - d!, li - d!)-(k + d!, li + d!), col, BF
NEXT k
END IF
dt! = TIMER - t0!
IF dt! > 1 THEN
LOCATE 4, 1: PRINT LTRIM$(STR$(gen&));
t0! = TIMER
END IF
END SUB