|
DataMuseum.dkPresents historical artifacts from the history of: RC4000/8000/9000 |
This is an automatic "excavation" of a thematic subset of
See our Wiki for more about RC4000/8000/9000 Excavated with: AutoArchaeologist - Free & Open Source Software. |
top - metrics - download
Length: 12288 (0x3000)
Types: TextFile
Names: »chstdemotx«
└─⟦621cfb9a2⟧ Bits:30002817 RC8000 Dump tape fra HCØ. Detaljer om "HC8000" projekt.
└─⟦0364f57e3⟧
└─⟦this⟧ »chstdemotx«
chem=set 200
scope day chem
chem=algol list.no
program for demonstration of chemical kinetics
Calahans method for stiff differential equations
version.5
fp: nmol.20 npl.30 blistx.false blistf.false bplot.true
test1.false test2.false testx.false testh.false testfj.false
plotter.<:tek4006c:>
26 01 83 11 00 00
begin
integer ntext,nmol,npl,ns,s;
boolean blistx,blistf,bplot,test1,test2,testx,testh,testfj;
real array PLOTTER(1:2);
readifp(<:ntext:>,ntext,100);
readifp(<:nmol:>,nmol,20);
readifp(<:npl:>,npl,30);
readbfp(<:blistx:>,blistx,false);
readbfp(<:blistf:>,blistf,false);
readbfp(<:bplot:>,bplot,true);
readbfp(<:testx:>,testx,false);
readbfp(<:testh:>,testh,false);
readbfp(<:test1:>,test1,false);
readbfp(<:test2:>,test2,false);
readbfp(<:testfj:>,testfj,false);
readsfp(<:plotter:>,PLOTTER,<:tek4006c:>);
s:= 1;
setplotname(string PLOTTER(increase(s)),0);
begin
integer mi,mn,ri,rn,vki,n,r;
integer array V,K(1:ntext),MOL,REA(1:nmol);
procedure writevkmr;
begin
integer i;
write(out,<:<10> V , K:>);
for i:= 1 step 1 until vki do
write(out,<:<10>:>,<< ddddddd>,V(i),<< d>,K(i));
write(out,<:<10>MOL:>);
for i:= 1 step 1 until mn do
write(out,<< d>,MOL(i));
write(out,<:<10>REA:>);
for i:= 1 step 1 until rn do
write(out,<< d>,REA(i));
end;
procedure syntaxerror;
begin write(out,<:<10>syntax error:>);
if test1 then write(out,
<:<10>vki = :>,vki,<: mn = :>,mn,<: rn = :>,rn);
if test1 then writevkmr;
alarm(<:<10>error:>,0);
end;
integer procedure searchmol(a);
integer a;
begin integer i,mi;
boolean found,comp;
found:= false; searchmol:= 0;
mi:= 0;
for mi:= mi+1 while MOL(mi)>0 and -,found do
begin
comp:= true;
for i:= 0,i+1 while K(a+i)=6 and K(MOL(mi)+i)=6 do
if V(a+i)<>V(MOL(mi)+i) then comp:= false;
if comp then
begin found:= true; searchmol:= mi end;
end
end;
begin
integer i;
integer array alfabet(0:255);
for i:= 0 step 1 until 9, 11 step 1 until 24,
26 step 1 until 42, 44, 46, 47, 58, 59, 60, 61,
63, 64, 94, 95, 96, 126 step 1 until 255 do
alfabet(i):= i;
for i:= 48 step 1 until 57 do alfabet(i):= 2 shift 12+i;
for i:= 65 step 1 until 93, 97 step 1 until 125 do
alfabet(i):= 6 shift 12 + i;
for i:= 10, 43, 45, 62 do alfabet(i):= 7 shift 12 + i;
alfabet(25):= 8 shift 12 + 25;
intable(alfabet);
write(out,<:reactions<10>:>); setposition(out,0,0);
i:= readall(in,V,K,1);
if i<0 then alarm(<:textbuf size:>,i);
unstackcuri;
intable(0);
end;
mn:= rn:= 1; vki:= 1;
cleararray(MOL); cleararray(REA);
REA(1):= 1;
if -,(V(1)=10 or K(1)=2 or K(1)=6 or V(1)=45) then
syntaxerror;
for vki:= vki while -,K(vki)=8 do
if K(vki)=2 and K(vki+1)=6 then
vki:= vki+1 else
if K(vki)=6 then
begin
if searchmol(vki)=0 then
begin MOL(mn):= vki; mn:= mn+1 end;
for vki:= vki+1 while K(vki)=6 do;
if V(vki)=62 and K(vki)<>7 and K(vki)<>8 then
syntaxerror;
end else
if V(vki)=10 and V(vki+1)=10 then vki:= vki+1 else
if V(vki)=10 and
(K(vki+1)=2 or K(vki+1)=6 or K(vki+1)=8 or V(vki+1)=45) then
begin vki:= vki+1; REA(rn):= vki end else
if V(vki)=43 and (K(vki+1)=2 or K(vki+1)=6) then
vki:= vki+1 else
if V(vki)=45 and V(vki+1)=62 and
(K(vki+2)=2 or K(vki+2)=6
or V(vki+2)=10 or K(vki+2)=8) then
begin rn:= rn+1; vki:= vki+2; REA(rn):= vki end else
syntaxerror;
n:= mn:= mn-1; r:= rn:= rn-1;
begin
integer i,rho,ny,ny1,ny2,ni,nm,jpl,mfac;
real xm,g,tt,hh,cpu,time;
boolean bfirst,reactants;
integer array N,R(1:n,1:r),LIST(1:n);
real array G,RK(1:r),XX,C(1:n),H(1:r,1:n),J(1:n,1:n),
PLT(0:npl),PLX(1:n,0:npl);
real array field f;
procedure writern;
begin write(out,<:<10><10>R:>);
for mi:= 1 step 1 until mn do
begin write(out,<:<10>:>);
for ri:= 1 step 1 until rn do
write(out,<< d>,R(mi,ri)) end;
write(out,<:<10><10>N:>);
for mi:= 1 step 1 until mn do
begin write(out,<:<10>:>);
for ri:= 1 step 1 until rn do
write(out,<< d>,N(mi,ri)) end;
end;
procedure outpl;
begin integer i,j;
if bplot then
begin
for i:= 1 step 1 until n do
begin plotmove(PLT(0),PLX(i,0)); pendown;
for j:= 1 step 1 until jpl-1 do
plotmove(PLT(j),PLX(i,j));
penup;
end;
plotend;
PLT(0):= PLT(npl);
for i:= 1 step 1 until n do
PLX(i,0):= PLX(i,npl);
end bplot;
jpl:= 1;
end;
procedure xin(X,tb,dt,te,h,epsr,epsa);
real tb,dt,te,h,epsr,epsa;
real array X;
begin
procedure plotnew;
begin
xm:= readr(<:maximal concentration in plot:>);
plotform(0,15,10);
scalexcoor(0,te,2,15-0.5);
scaleycoor(0,xm,1,10-2.5);
plotframe(0.0,0.0);
plotend;
end;
procedure parameters;
begin
dt:= readr(<:time interval:>);
te:= readr(<:maximal time:>);
epsr:= readr(<:relative accuracy:>);
epsa:= readr(<:absolute accuracy:>);
end;
procedure listing;
begin write(out,<:<10>listing af species<10><10>:>);
for mi:= 1 step 1 until mn do
begin f:= 2*MOL(mi)-2; mfac:= 1;
write(out,string V.f(increase(mfac)),<: = :>);
setposition(out,0,0); read(in,LIST(mi)); end;
end;
procedure concentrations;
begin write(out,<:<10>start concentrations<10><10>:>);
for mi:= 1 step 1 until mn do
begin f:= 2*MOL(mi)-2; mfac:= 1;
write(out,string V.f(increase(mfac)),<: = :>);
setposition(out,0,0); read(in,C(mi));
if C(mi)<0 then goto lb; end;
end;
procedure rateconstants;
begin write(out,<:<10>rateconstants<10><10>:>);
for ri:= 1 step 1 until rn do
begin
write(out,<<dd>,ri,<:: :>);
setposition(out,0,0); read(in,RK(ri));
if RK(ri)<0 then goto lb; end;
end;
if bfirst then
begin rateconstants; concentrations;
bfirst:= false;
if blistx then listing; parameters;
for ny:= 1 step 1 until n do X(ny):= C(ny);
PLT(0):= tb:= 0; jpl:= 1;
for ny:= 1 step 1 until n do PLX(ny,0):= X(ny);
h:= dt/2;
if bplot then plotnew; end
else begin
if readb(<:continue:>) then
begin tb:= tt; h:= hh;
for ny:= 1 step 1 until n do X(ny):= XX(ny)
end
else begin if readb(<:new start concentrations:>) then
concentrations;
for ny:= 1 step 1 until n do X(ny):= C(ny);
PLT(0):= tb:= 0; jpl:= 1;
for ny:= 1 step 1 until n do PLX(ny,0):= X(ny);
end continue;
if readb(<:new rateconstants:>) then rateconstants;
if blistx then
begin if readb(<:new listing of species:>) then listing end;
parameters;
if bplot then begin if readb(<:new plot:>) then plotnew end;
h:= dt/2;
end bfirst;
cpu:= systime(1,0,time);
end xin;
procedure xout(X,t,h,s);
integer s;
real t,h;
real array X;
begin
procedure listf;
begin
integer i,j;
boolean bf;
real array F(1:n),J(1:1);
fjac(false,X,t,F,J);
bf:= true; j:= 0;
for i:= 1 step 1 until n do
if (LIST(i)=1 and s=0) or (LIST(i)=2 and abs(s)=i) then
begin if bf then write(out,"nl",1,<:F =:>,"sp",8);
if j=5 or j=10 then write(out,"nl",1,"sp",11);
j:= j+1; bf:= false; write(out,<< -d.dddd'-dd>,F(i));
end;
if j>0 then setposition(out,0,0);
end;
procedure writefjac(bjac);
boolean bjac;
begin real array F(1:n),J(1:n,1:n);
fjac(bjac,X,t,F,J);
write(out,<:<10>F<10>:>);
for ny:= 1 step 1 until n do
begin write(out,<< -d.dddd'-dd>,F(ny));
if ny=5 or ny=10 then write(out,"nl",1);
end;
if bjac then
begin write(out,<:<10>J<10>:>);
for ny1:= 1 step 1 until n do
begin
for ny2:= 1 step 1 until n do
begin write(out,<< -d.dddd'-dd>,J(ny1,ny2));
if ny2=5 or ny2=10 then write(out,"nl",1);
end;
write(out,"nl",2);
end
end
end writefjac;
procedure listx;
begin integer i,j;
boolean btime;
btime:= true; j:= 0;
for i:= 1 step 1 until n do
if (LIST(i)=1 and s=0) or (LIST(i)=2 and abs(s)=i) then
begin if btime then write(out,"nl",1,<<d.ddd'-dd>,t,"sp",2);
if j=5 or j=10 then write(out,"nl",1,"sp",11);
j:= j+1; btime:= false; write(out,<< -d.dddd'-dd>,X(i));
end;
if j>0 then setposition(out,0,0);
end listx;
procedure writex;
begin integer i;
write(out,<:<10>:>,<<-d>,s,"sp",1,<<d.ddd'-dd>,t,"sp",2);
for i:= 1 step 1 until n do
begin if i=5 or i=9 then write(out,"nl",1,"sp",14);
write(out,<< d.dddd'-dd>,X(i));
end; setposition(out,0,0);
end;
PLT(jpl):= tt:= t; hh:= h;
for ny:= 1 step 1 until n do PLX(ny,jpl):= XX(ny):= X(ny);
jpl:= jpl+1;
if jpl>npl then outpl;
if blistx then listx;
if blistf then listf;
if testx then writex;
if testh then
begin write(out,<:<10>h = :>,<<d.d'-dd>,hh); setposition(out,0,0) end;
if testfj and s=0 then writefjac(true);
if opkey then goto le;
end xout;
procedure fjac(bjac,X,t,F,J);
boolean bjac;
real t;
real array X,F,J;
begin
for rho:= 1 step 1 until r do
begin G(rho):= RK(rho);
for ny:= 1 step 1 until n do
if R(ny,rho)=0 then else
if R(ny,rho)=1 then G(rho):= G(rho)*X(ny) else
if R(ny,rho)=2 then G(rho):= G(rho)*X(ny)*X(ny) else
G(rho):= G(rho)*X(ny)**R(ny,rho);
if bjac then
for ny2:= 1 step 1 until n do
if R(ny2,rho)=0 then H(rho,ny2):= 0
else if X(ny2)<>0 then
H(rho,ny2):= G(rho)*R(ny2,rho)/X(ny2)
else if R(ny2,rho)>1 then H(rho,ny2):= 0
else begin g:= RK(rho);
for ny:= 1 step 1 until ny2-1,
ny2+1 step 1 until n do
if R(ny,rho)=0 then
else if R(ny,rho)=1 then g:= g*X(ny)
else if R(ny,rho)=2 then g:= g*X(ny)*X(ny)
else g:= g*X(ny)**R(ny,rho);
H(rho,ny2):= g;
end
end;
for ny:= 1 step 1 until n do
begin F(ny):= 0;
for rho:= 1 step 1 until r do
F(ny):= F(ny)+N(ny,rho)*G(rho) end;
if bjac then
for ny1:= 1 step 1 until n do
for ny2:= 1 step 1 until n do
begin J(ny1,ny2):= 0;
for rho:= 1 step 1 until r do
J(ny1,ny2):= J(ny1,ny2)+N(ny1,rho)*H(rho,ny2);
end;
end;
cleararray(N); cleararray(R);
for ri:= 1 step 1 until rn do
begin vki:= REA(ri); mfac:= 1; reactants:= true;
for vki:= vki while -,(K(vki)=8 or V(vki)=10
or(V(vki)=45 and -,reactants)) do
if K(vki)=2 then
begin mfac:= V(vki); vki:= vki+1 end else
if K(vki)=6 then
begin mi:= searchmol(vki);
if reactants then R(mi,ri):= R(mi,ri)+mfac else
N(mi,ri):= N(mi,ri)+mfac;
mfac:= 1; for vki:= vki+1 while K(vki)=6 do;
end else
if V(vki)=43 then vki:= vki+1 else
if V(vki)=45 then
begin reactants:= false; vki:= vki+2 end else
begin write(out,<:<10>programerror:>);
writern; alarm(<:<10>error:>,0); end;
end;
if test1 then begin writevkmr; writern end;
for mi:= 1 step 1 until mn do
for ri:= 1 step 1 until rn do
N(mi,ri):= N(mi,ri)-R(mi,ri);
bfirst:= true;
if test2 then writern;
lb:chstiffpp(n,fjac,xin,xout);
le:cpu:= systime(1,time,time)-cpu;
write(out,<<dddd.dd>,<:<10><10>cpu = :>,cpu,<: time = :>,time);
if bplot then begin outpl; plotend end;
if -,readb(<:end:>) then goto lb;
lend:if bplot then plotclose;
end;
end;
end
▶EOF◀