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"}
procedure set_MENU_choices; {this procedure (=subroutine) asks the user to assign values}
begin
writeln('What should the saturation concentration be? (50)');
writeln('(The numbers in parentheses are the values used in the original program.');
readln(SATURATION);
writeln('What should the supersaturation concentration be? (100)');
readln(SUPERSATURATION);
writeln('What diffusion constant should A have? (0.5)');
readln( DiffConst_A);
writeln('What diffusion constant should B have? (0.5)');
readln(DiffConst_B);
writeln('The initial concentration of A should be what? (200)');
readln(ACONC);
writeln('The initial concentration of A should be what? (20)');
readln(BCONC);
writeln('A new graph should be drawn every how many recalculation cycles? (5)');
skip:=5; readln(skip);
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_MENU_choices;
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;
readln;
end. {There, now! You understood it all, didn't you?}
|