#### Computer simulation of formation of Liesegang bands (written in Pascal)

```

program Liesegang_With_Menu;                {by Albert Harris, 1994}

uses memtypes, quickdraw;                    {in Pascal programs, you start by }
const              MN = 200; SPACE=2;        {naming all the variables you will use}
type ary = array [0..MN] of real;            {and whether they will be integers, etc.}
var
SATURATION, SUPERSATURATION, DiffConst_A, DiffConst_B, ACONC, BCONC : real;
Skip : integer;
a, b, na, nb : ary;
n, cycle : integer;
c :array [0..MN] of boolean;             {arrays are sort of like rows of mailboxes}
t :array [0..MN] of integer;             {in which you can store a value in each "box"}

begin
writeln('What should the saturation concentration be? (50)');
writeln('(The numbers in parentheses are the values used in the original program.');

writeln('What should the supersaturation concentration be? (100)');

writeln('What diffusion constant should A have? (0.5)');

writeln('What diffusion constant should B have? (0.5)');

writeln('The initial concentration of A should be what? (200)');

writeln('The initial concentration of A should be what? (20)');

writeln('A new graph should be drawn every how many recalculation cycles? (5)');
end;

procedure set_start;  {this procedure assigns initial values to appropriate "mailboxes"}
begin
for n:=1 to MN do begin;
a[n]:= 0.0;
b[n]:=BCONC;
c[n]:=false;
cycle:= -1;
end;
a[1]:=ACONC  ;b[1]:=0;
end; {of the set_start procedure}

procedure diffuse; {this procedure recalculates values as if they were diffusing}
begin;
for n:=2 to MN-1 do begin;
na[n]:=a[n]+ DiffConst_A * (a[n-1]+a[n+1]-2*a[n])/2;
nb[n]:=b[n]+ DiffConst_B * (b[n-1]+b[n+1]-2*b[n])/2;
end;
end; {of the diffuse procedure}

procedure plot_sisyphisAxB;   {this procedure draws one of the graphs of concentrations}
begin;
moveto(250+ SPACE,250-round(a[1]*b[1]/5));
for n:=2 to MN do begin;
lineto(250+ n*SPACE,250-round(a[n]*b[n]/5));
end;
end; {you don't really have to write which procedure it's the end of}
procedure react;  {this procedure causes crystalization where concentrations are high enough}
begin; cycle:= cycle+1;
for n:=2 to MN-1 do begin;
if a[n]*b[n]>SUPERSATURATION then c[n]:=true ;
if c[n]=true then
begin;
if a[n]*b[n]>SATURATION then
begin
t[n]:=cycle;
{plot_sisyphisAxB;}
repeat
a[n]:=a[n]-1;
b[n]:=b[n]-1;
until a[n]*b[n]< SATURATION {or b[n]<0.0};
end;
end;
end;
end; {of the react procedure, although the program would run just as well without this}

procedure new_old; {this procedure resets values, before another round of recalculation}
begin;
for n:=2 to MN-1 do begin;
a[n]:= na[n];
b[n]:= nb[n];
end;
end;  { can you figure out which procedure this is the end of? }

procedure plot_A;    {another graph plotting procedure}
begin;
moveto(10+ SPACE,150- round(a[1]));
for n:=2 to MN do begin;
lineto(10+ n*SPACE,150- round( a[n]));
end;
end; {of procedure plot_A}

procedure plot_B;   {yet another graph plotting procedure}
begin;
moveto(10+ SPACE,150-round( b[1]));
for n:=2 to MN do begin;
lineto(10+ n*SPACE,150-round( b[n]));
end;
end;  {of procedure plot_B}

procedure plot_C;  {this procedure marks location where crystal precipitation occurred}
begin;
for n:=2 to MN do begin;
if c[n]=true then
begin;
moveto(10+ n*SPACE,155);
lineto(10+ n*SPACE,170);
end;
end;
end;  {of procedure plot_C}

procedure plot_AxB;  {yet another graph plotting procedure}
begin;
moveto(10+SPACE, 250-round(SUPERSATURATION/5));
lineto(10+SPACE*MN,250-round(SUPERSATURATION/5));

moveto(10+SPACE, 250-round(SATURATION/5));
lineto(10+SPACE*MN,250-round(SATURATION/5));

moveto(10+ SPACE,250-round(a[1]*b[1]/5));
for n:=2 to MN do begin;
lineto(10+ n*SPACE,250-round(a[n]*b[n]/5));
end;
end; {of procedure plot_AxB}

begin  {main program};  {now the program uses all the subroutines defined above}
set_start;
while not keypressed do
begin;  {the following is a loop that keeps repeating forever, until you type}
react;
diffuse;
new_old;
if cycle mod skip=0 then
begin
clearscreen;
writeln('                                                   cycle number ',cycle);
writeln('                                        Diffusion constant of A   ',DiffConst_A:2:2);
writeln('                                        Diffusion constant of B   ',DiffConst_B:2:2);
writeln('                                       A concentration at margin  ',ACONC:2:2);
writeln('                                       concentration of B in gel  ',BCONC:2:2);
writeln('                                                     SATURATION =  ',SATURATION:2:2);
writeln('                                                SUPERSATURATION =  ',SUPERSATURATION:2:2);
plot_A;
plot_B;
plot_C;
plot_AxB;
end;end;
readln; writeln('The following are the positions and times at which rings formed.');
for n:= 1 to MN do begin

if c[n]=true then
writeln('                                           n= ',n,'  time = ',t[n]);
end;