#### reacteqn.red4.5 KB Permalink History Raw

 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147 ``````module reacteqn; % REDUCE support for reaction equations. % Author: H. Melenk % January 1991 % Copyright (c) Konrad-Zuse-Zentrum Berlin, all rights reserved. % Introduce operators for chemical equations. algebraic operator rightarrow; newtok '((!- !>) rightarrow); infix rightarrow; precedence rightarrow,equal; algebraic operator doublearrow; newtok '((!< !>) doublearrow); infix doublearrow; precedence doublearrow,equal; algebraic operator rate; global '(species); share species; global '(rates); share rates; put('reac2ode,'psopfn,'r2oeval); symbolic procedure r2oeval u; begin scalar r,k,x,rhs,lhs,ratel,odel,oldorder,lhsl,rhsl; integer rc; if eqcar(species,'list) then odel:=for each x in cdr species collect reval x . 0; u := reval car u; if not eqcar(u,'list) then typerr(u,"list of reactions"); u := cdr u; loop: if null u then goto finis; r := reval car u; u := cdr u; if not pairp r or not memq(car r,'(rightarrow doublearrow)) then goto synerror; lhs := r2speclist cadr r; rhs := r2speclist caddr r; % include new species for each x in append(lhs,rhs) do odel:=r2oaddspecies(cdr x,odel); % generate contribution from forward reaction. k := if u and (x:=reval car u) and not(pairp x and memq(car x,'(rightarrow doublearrow))) then <> else list('rate,rc:=rc+1); ratel := k . ratel; r2oreaction(lhs,rhs,k,odel); % eventually generate backward reaction if car r='doublearrow then <> else list('rate,rc:=rc+1); ratel := k . ratel; r2oreaction(rhs,lhs,k,odel); >>; lhsl := lhs.lhsl; rhsl := rhs.rhsl; goto loop; finis: ratel := reversip ratel; rates := 'list. ratel; for each x in ratel do if numberp x or pairp x and get(car x,'dname) then ratel := delete(x,ratel); species := 'list. for each x in odel collect car x; r2omat(cdr species,reversip lhsl,reversip rhsl); for each r in ratel do if not idp r then ratel:=delete(r,ratel); if ratel then eval list('order,mkquote ratel); oldorder := setkorder append(ratel,cdr species); odel := 'list . for each x in odel collect list('equal,list('df,car x,'t),reval cdr x); setkorder oldorder; return odel; synerror: typerr(r,"reaction"); end; symbolic procedure r2omat(sp,lhsl,rhsl); % construct input and output matrices in REDUCE syntax. begin scalar m; integer nreac,nspec,j; nspec := length sp; nreac:= length lhsl; apply ('matrix,list list list('inputmat,nreac,nspec)); apply ('matrix,list list list('outputmat,nreac,nspec)); for i:=1:nreac do << for each x in nth(lhsl,i) do <>; for each x in nth(rhsl,i) do <>; >>; end; symbolic procedure r2findindex(a,l); r2findindex1(a,l,1); symbolic procedure r2findindex1(a,l,n); if null l then rederr "index not found" else if a=car l then n else r2findindex1(a,cdr l,n+1); symbolic procedure r2speclist u; % convert lhs/rhs to a list of pairs (multiplicity . spec). <>; symbolic procedure r2speclist1 x; if eqcar(x,'times) then r2speclist2(cadr x,caddr x,cdddr x) else 1 . x; symbolic procedure r2speclist2(x1,x2,rst); if not null rst or not fixp x1 and not fixp x2 then typerr(append(list('times,x1,x2),rst),"species") else if fixp x1 then x1.x2 else x2.x1; symbolic procedure r2oaddspecies(s,odel); % generate a new (empty) equation for a new species. if assoc(s,odel) then odel else <>; symbolic procedure r2oreaction(lhs,rhs,k,odel); % add the contribution of one reaction to the ode's. begin scalar coeff,e; coeff := k; for each x in lhs do coeff:=aeval list('times,coeff,list('expt,cdr x,car x)); for each x in lhs do <>; for each x in rhs do <>; return odel; end; endmodule; end; ``````