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]