COMMENT THESE ROUTINES WERE SCANNED IN 1998 FROM LASER 9700 OUTPUT;
COMMENT
COPYRIGHT (C) 1981
COLUMBIA UNIVERSITY BIOENGINEERING INSTITUTE
&
VASOS-PETER JOHN PANAGIOTOPOULOS II
TECHNICAL STANFORD-SAIL-ALGOL SUBROUTINE PACKAGE ONE SBRPC1
This project was funded in part (since April, 1979) by
NSF Rat Skull Grant-Professor Melvin L. Moss-Fall 1979 only
NIH HD14371 & HL16851 Professor Richard Skalak
CU-CCA NCA-ID fund and other funds for educational computing
Samani International Enterprises Siengg Computational Research Fund
Vasos-Feter John Panagiotopoulos II
I wish to acknowledge the unlimited and most kind assistance of
Professor Gautam Dasgupta.
References utilised:
Numerical Analysis, Scheid, Schaum/McGrawHill
Applied Numerical Methods, Carnahan, Luther&Wilkes,Wiley
Introductory Computer Methods and Numerical Analysis,Pennington
Advanced Calculus for Applications,Hildebrand
Handbook of Mathematical Functions,Abramowitz&Stegun
The following Columbia courses:
ChE3100X79,ChE4510Y81,ChE4520YS1
ENTRY Factorial,Lagrange!lnterpolation,Average,Norm,Quadratic,matsum,
Summation!Of!Function,secant!root,Newton!Raphson,Regula!Falsi,
Power!Method,RUNGEK,realrd,read!arai,File!Read!Arai,Matrix!Print,
FIBONACCI,gauss!jordan,Simpson,Recirculation,Max!Array,Min!Array,
InstaPlot,cc!trapper,sink,Structural!Section,Monte!Carlo,
Binary!Monte!Carlo,Gradient!Search,Gamma,Legendre!P,
Weighted!Gauss!Legendre!Integration,
C!Conjugate,C!Add,C!Subtract,C!Abs,C!Multiply,C!Divide,C!ExP,CvCmplx
;
BEGIN "SBRPCl"
Require "sai:abbrev.sai" source!file;comment sail abbreviations package.;
Require "omni:disprm.sai" source!file;comment sail omnigraph package.;
Require "sai:routin.hdr" source!file;comment sail routine package.;
Require "sai:ttyio.hdr" source!file;comment sail i/o package.;
Define Tan(x)= {sin(x)/cos(x)},
Degrees (radians)={radians*180/4/atan(1)},
Round!Off(x)= {2*(x%2)},
Even!Odd(x)= {2*abs(x/2-x%2)},
Kroenecker(i,j)= {if i=j then 1 else 0},
Heaviside(x)= {if x>O then 1 else 0},
Dirac(x)= {if x=O then 1 else 0},
Permutation(i.J,k)={if i=j or j=k or i=k then 0 else
if (i=1 and j=2 and k=3) or(i=2 and j=3 and k=1)
or (i=3 and j=l and k=2)
then 1 else -1},
ARCOS(x)={4*ATAN(1)+ATAN((1+x^2)^.5/X)},
ARSIN(x)={-ARCOS((1-x^2)^.5)},
INFINITY={.1701412@39},
COMPLEX= {RECORD!POINTER(CMPLX)},
NEW!COMPLEX= {NEW!RECORD(CMLX)};
RECORD!CLASS CMPLX (REAL Re.lm);
Internal Real Procedure Factorial (Integer x);
BEGIN "Factorial Recursive (self-calling) procedure"
return(if x=O then I else x*factorial (x-1));
END "Factorial Recursive (self-calling) procedure";
Internal Real procedure Lagrange!Interpolation (Real Array x,y;
Integer n;
Real xx);
begin"Lagrange!Interpolation"
real array 1 [1:n]
integer i,j;
real yy;
for i:=1 step 1 until n do
begin "Lagrange looping i=/=j"
for j:=1 step 1 until n do if (j neq i) then
l[i] =l[i]*(xx-x(j])/(x[i]-x[j])
yy:=y[i]*l[i]+yy;
end "Lagrange looping i=/=j";
return (yy);
end"Lagrange!Interpolation";
Internal Real Procedure Integral!of!Polynomial!Term(real a,b,n,u1,u2);
begin "Integral Of Polynomial Term"
comment Integral of (A+BU)+N from U1 to U2;
return(if n neq -1 then((a+b*u2)f(n+1))/(n+1)/b-((a+b*u1)+(n+1))/(n+1)/b
else 1/b*(log(a+b*u2)-log(a+b*u1));
end "lntegral Of Polynomial Term";
Internal Real Procedure Average(real array a:integer n);
begin "avg"
real sum; integer m;m:=n;
while (n:=n-1) do sum:=sum+a[n];
return(sum/m);
end "avg";
Internal Real Procedure Norm (Real Array Vector;Integer Rank);
BEGIN "Norm of vector "
REAL Scalar!Norm;
while (rank:=rank-1) do Scalar!Norm:=Scalar!Norm+Vector[rank]+2;
return(Scalar!Norm^.5);
END "Norm of vector ";
Internal Procedure Quadratic(real a,b,c;
reference real first!root, second!root;
reference string root!type);
begin "Quadratic"
real g,d;
g:=b^2-4*a*c;
if (a=O) then
begin "real1"
first!root:=second!root:=-c/b;
root!type:="single real root";
end "real1";
else
begin "two-root"
if (g>0) then
begin "real2"
d:=sqrt(g);
first!root:=(-b-d)/(2.*a);
second!root:=( b-d)/(2.*a);
root!type:="double real roots";
end "real2"
else
begin "cmp2"
first!root:=-b/2./a;
second!root:=sqrt(-g)/2/a;
root!type:="imaginary first +/-second*i";
end "cmp2";
end "two-root";
return;
end "Quadratic";
Internal Real procedure matsum (real array x; integer dim);
begin"matsum" real ssum;
for dim:=dim step -l until 1 do ssum:=ssum+x [dim] ;return(ssum)
end"matsum";
Internal Real procedure
Summation!Of!Function(INTEGER Start!Index,End!Index,Interval;
REAL X;
REAL PROCEDURE Funct);
begin"Summation!Of!Function"
integer icount;real g;
for icount:=start!index step Interval until End!lndex do
g:=g+funct(x,icount);
return (g)
end"Summation!Of!Function";
Internal Real Procedure secant!root(REAL X0,X1,Accuracy;REAL PROCEDURE F);
begin"sec" real x2;
while (abs(f(xO))geq accuracy ) do
x0:=(f(x1:=x0)*(x2:=x1)*f(x2)*x1)
/(f(x1)-f(x2));
return (xO)
end"sec";
Internal Real procedure Newton!Raphson(REAL X0,X1,Accuracy;REAL PROCEDURE F);
Begin "Newton Raphson Root Determination Reverse diff'
Real x2;
While (abs(f(x2)-f(x1)) geq accuracy) do
x0:=((x2:=x1)-(x1:=x0))*f(x1)/(f(x1)-f(x2))+x1;
Return (xO)
End"Newton Raphson Root Determination Reverse diff";
Internal Real procedure Regula!Falsi(REAL XO,Xl,Accuracy;REAL PROCEDURE F):
Begin"Regula Falsi According to Firouztale, Spencer & Wright"
Real x2;
While (abs(f(x2)-f(x1)) geq accuracy) do
x0:=((x2:=x1)-(x1:=x0))*f(x1)/(f(x1)-f(x2))+x1;
Return (x0);
End"Regula Falsi According to Firouztale, Spencer & Wright";
lnternai Procedure Power!Method( Real array a;
Reference Real Array Eigen!Vector;
Integer Rank;
Reference Real Eigen!Value;
Real Accuracy);
BEGIN "POWERM: Power Method of Mises for Eigen Systems"
Comment Panagiotopoulos Strikes Again
Real Array Temp!Eigen!Pseudo[1:Rank];
Integer I,J;
Real Old!Eigen!Value;
While ( Eigen!Value-Old!Eigen!Value>Accuracy) do
BEGIN "Power Method Inner Block"
Old!Eigen!Value:=Eigen!Value;
Eigen!Value:=0;
For I:=1 step 1 until Rank do
BEGIN "Mises"
For J:=1 step 1 until Rank do
Temp!Eigen!Pseudo[i]:=Temp!Eigen!Pseudo[i]
+a[j,3]*Eigen!Vector[3];
Eigen!value:=Eigen!Value+Temp!Eigen!Pseudo[i]^2;
END "Mises";
Eigen!Value:=Eigen!Value^.5;
For I:=1 step 1 until Rank do
Eigen!Vector[i]:=Temp!Eigen!Pseudo[i]/Eigen!Value:
END "Power Method Inner Block";
If (Eigen!Vector[1) neq 0) then
While(Rank:=Rank-1) do
Eigen!Vector[Rank]:=Eigen!Vector[Rank]/Eigen!Vector[1];
Return;
END "POWERM: Power Method of Mises for Eigen Systems";
Internal Procedure RUNGEK(REFERENCE REAL ARRAY X;
INTEGER N!Equations;
REAL PROCEDURE DRV;
REAL Delta!X, End!x);
BEGIN "RUNGEK: Runge-Kutta-Gill fourth-order Integration Program"
Integer i,j,k;
Real XX;
Comment : 1)Derived from Taylor expansion of ODE
2)DRV must return derivative of Jth function
3)X has independant variable&'solutions
4)IC's must be entered as x[O, _];
for xx:=x[O,C] step DeltaIX until EndIX do
BEGIN "Integration Stepping"
Real Array C[1:4,1:N!Equations];
i:=i+1;x[i,0]:=xx;
for j:=1 step 1 until N!Equations do
BEGIN "Gill constants"
c[1,j]:=DRV(x[i,0],x[i,j],j);
c[2,j]:=DRV(x[i,0]+.5*Delta!X,x[i,j]+.5*Delta!X*C[1,j],j);
c[3,j]:=DRV(x[i,0]+.5*Delta!X,x[i,j]+.2071068*Delta!X*C[1,j]
+.2928932*Delta!X*C[2,j] ,j);
c[4,j]:=DRV(x[i,0]+Delta!X,x[i,j]-.7071068*Delta!X*C[1,j]
+1.7071068*Delta!X*C[3,j] ,j);
c[i+1,j]:=x[i,j]
+Delta!X/6*
( C[1,j]
+.2928932*C[2,j]
+3.4142136*c[3,j]
+C[4,j]
);
END "Gill constants":
for j:=1 step 1 until N!Equations do
If(ABS( (C[3,j]-C[2,j])
/(C[2,j]-C[1,j])
)
>.01
)
Then Delta!X:=Delta!x/10;
END "Integration Stepping";
END "RUNGEK: Runge-Kutta-Gill fourth-order Integration Program";
Internal Real procedure realrd(string s);
begin "realrd"
! input real constants and assign to real variables;
real f; string a;
print(13&10&"Please input ",s,"> ");
a:=inchwl;f:=realscan(a,32);
return (f)
end "realrd";
Internal Procedure read!arai (reference real array a; integer n;string name);
begin"araird" string ary; integer i;
!similarly input arrays;
print (13&10&"Type in array,"&name&", elements sep by space,",
" when finished;",n,"elements "&13&10&" > ");
ary:=inchwl;
for i:=1 step 1 until n do a[i]:=realscan(ary,32);
end"araird";
Internal Procedure File!Read!Arai (REFERENCE REAL Array A;
INTEGER N,Channel,BreakTable);
begin"array-file" STRING Ary; INTEGER I;Ary:=NULL;
for i:=0 step 1 until n do
begin "ary-in"
if ary=null then ary:=input(Channel,Breaktable);
a[i]:=realscan(ary,32)
end "ary-in";
end"array-file";
Internal Simple Procedure Matrix!Print(REAL ARRAY A;
INTEGER Down,Accross;
STRING Name);
BEGIN "Print a two-dimensional array" Integer I,J;
Print (10&13&9&"Matrix:"&Name&13&10);
For I:=1 step 1 until Down do
BEGIN "Matrix looping"
for j:=1 step 1 until accross do
Print (a[i,j] ,32&null)
Print (13&10&Null);
END "Matrix looping";
END "Print a two-dimensional array";
Internal Real Procedure FIBONACCI(REAL Aa,Bb,Accuracy;
REAL PROCEDURE F);
begin"fibonacci"
! MAXIMISATION -give neg fct for min;
real new!fibonacci,fibonacci!old,temporary,tau,l,r,a,b;
fibonacci!old:=0;new!fibonacci:=1;l:=a:=aa;r:=b:=bb;
while abs((l-r)/r)>accuracy do
begin"fibonacci section"
new!fibonacci:=fibonacci!old+new!fibonacci
+O*(fibonacci!old:=new!fibonacci);
tau:=fibonacci!old/new!fibonacci;
l:=b-tau*abs(b-a);r:=a+tau*abs (b-a);
if f(l)>f(r) then b:=r else a:=l;
if abb then done;
end"fibonacci section";
return((a+b)/2);
end"fibonacci";
Internal Procedure gauss!jordan (reference real array a,bb;integer n);
begin "gauss" real array c[1:n,1:n+1],b,d,e [1:n+1]
integer i,ipl,j,m,k,l,ko,km;real r,rl,s;label unitise,skip;
comment The Panagiotopoulos-Snyder matrix zonker;
for i:=1 step 1 until n do for j:=i step 1 until n do c[i,j]:=a[i,j];
for i:=1 step 1 until n do c[i,n+1]:=bb[i];
for i:=1 step 1 until n do
begin "gj1"
r:=c[i,i]; if (i geq n) then go unitise; m:=i;ip1:=i+1;
for j:=ip1 step 1 until n do
begin "gj2"
r1:=abs(c[m,i]); s:=abs(c[j,i]);if (r1Maximum then Maximum:=Our!Vector[Dimension];
Return (Maximum)
END "Maximum of Vector";
Internal Real Procedure Min!Array(Real Array Our!Vector;
Integer Dimension);
BEGIN "Minimum of Vector" Real Minimum: Minimum:=INFlNITY;
Dimension:=Dimension+1;
While (Dimension:=Dimension-1) do
if Our!Vector[Dimension]x!upper[i] or x[i]