showtime;

comment some examples of the FOR statement;

comment summing the squares of the even positive integers
        through 50;

for i:=2 step 2 until 50 sum i**2;

comment to set  w  to the factorial of 10;

w := for i:=1:10 product i;

comment alternatively, we could set the elements a(i) of the
        array  a  to the factorial of i by the statements;

array a(10);

a(0):=1$

for i:=1:10 do a(i):=i*a(i-1);

comment the above version of the FOR statement does not return
        an algebraic value, but we can now use these array
        elements as factorials in expressions, e. g.;

1+a(5);

comment we could have printed the values of each a(i)
        as they were computed by writing the FOR statement as;

for i:=1:10 do write a(i):= i*a(i-1);

comment another way to use factorials would be to introduce an
operator FAC by an integer procedure as follows;

integer procedure fac (n);
   begin integer m;
        m:=1;
    l1: if n=0 then return m;
        m:=m*n;
        n:=n-1;
        go to l1
   end;

comment we can now use  fac  as an operator in expressions, e. g.;

z**2+fac(4)-2*fac 2*y;

comment note in the above example that the parentheses around
the arguments of FAC may be omitted since it is a unary operator;

comment the following examples illustrate the solution of some
        complete problems;

comment the f and g series (ref  Sconzo, P., Leschack, A. R. and
         Tobey, R. G., Astronomical Journal, Vol 70 (May 1965);

deps:= -sig*(mu+2*eps)$
dmu:= -3*mu*sig$
dsig:= eps-2*sig**2$
f1:= 1$
g1:= 0$
 
for i:= 1:8 do 
 <<f2:= -mu*g1 + deps*df(f1,eps) + dmu*df(f1,mu) + dsig*df(f1,sig);
   write "F(",i,") := ",f2;
   g2:= f1 + deps*df(g1,eps) + dmu*df(g1,mu) + dsig*df(g1,sig);
   write "G(",i,") := ",g2;
   f1:=f2;
   g1:=g2>>;

comment a problem in Fourier analysis;

factor cos,sin;

on list;

(a1*cos(wt) + a3*cos(3*wt) + b1*sin(wt) + b3*sin(3*wt))**3
        where {cos(~x)*cos(~y) => (cos(x+y)+cos(x-y))/2,
               cos(~x)*sin(~y) => (sin(x+y)-sin(x-y))/2,
               sin(~x)*sin(~y) => (cos(x-y)-cos(x+y))/2,
               cos(~x)**2 => (1+cos(2*x))/2,
               sin(~x)**2 => (1-cos(2*x))/2};

remfac cos,sin;

off list;

comment end of Fourier analysis example;

comment the following program, written in  collaboration  with  David
Barton  and  John  Fitch,  solves a problem in general relativity. it
will compute the Einstein tensor from any given metric;

on nero;

comment here we introduce the covariant and contravariant metrics;

operator p1,q1,x;

array gg(3,3),h(3,3);

gg(0,0):=e**(q1(x(1)))$
gg(1,1):=-e**(p1(x(1)))$
gg(2,2):=-x(1)**2$
gg(3,3):=-x(1)**2*sin(x(2))**2$

for i:=0:3 do h(i,i):=1/gg(i,i);

comment generate Christoffel symbols and store in arrays cs1 and cs2;

array cs1(3,3,3),cs2(3,3,3);

for i:=0:3 do for j:=i:3 do
   <<for k:=0:3 do
        cs1(j,i,k) := cs1(i,j,k):=(df(gg(i,k),x(j))+df(gg(j,k),x(i))
                                     -df(gg(i,j),x(k)))/2;
        for k:=0:3 do cs2(j,i,k):= cs2(i,j,k) := for p := 0:3
                                 sum h(k,p)*cs1(i,j,p)>>;

comment now compute the Riemann tensor and store in r(i,j,k,l);

array r(3,3,3,3);

for i:=0:3 do for j:=i+1:3 do for k:=i:3 do
   for l:=k+1:if k=i then j else 3 do
      <<r(j,i,l,k) := r(i,j,k,l) := for q := 0:3
                sum gg(i,q)*(df(cs2(k,j,q),x(l))-df(cs2(j,l,q),x(k))
                + for p:=0:3 sum (cs2(p,l,q)*cs2(k,j,p)
                        -cs2(p,k,q)*cs2(l,j,p)));
        r(i,j,l,k) := -r(i,j,k,l);
        r(j,i,k,l) := -r(i,j,k,l);
        if i neq k or j>l
          then <<r(k,l,i,j) := r(l,k,j,i) := r(i,j,k,l);
                 r(l,k,i,j) := -r(i,j,k,l);
                 r(k,l,j,i) := -r(i,j,k,l)>>>>;

comment now compute and print the Ricci tensor;

array ricci(3,3);

for i:=0:3 do for j:=0:3 do  
    write ricci(j,i) := ricci(i,j) := for p := 0:3 sum for q := 0:3
                                        sum h(p,q)*r(q,i,p,j);

comment now compute and print the Ricci scalar;

rs := for i:= 0:3 sum for j:= 0:3 sum h(i,j)*ricci(i,j);

comment finally compute and print the Einstein tensor;

array einstein(3,3);

for i:=0:3 do for j:=0:3 do
         write einstein(i,j):=ricci(i,j)-rs*gg(i,j)/2;

comment end of Einstein tensor program;

clear gg,h,cs1,cs2,r,ricci,einstein;

comment an example using the matrix facility;

matrix xx,yy,zz;

let xx= mat((a11,a12),(a21,a22)),
   yy= mat((y1),(y2));

2*det xx - 3*w;

zz:= xx**(-1)*yy;

1/xx**2;

comment end of matrix examples;

comment a physics example;

on div; comment this gives us output in same form as Bjorken and Drell;

mass ki= 0, kf= 0, p1= m, pf= m;

vector ei,ef;

mshell ki,kf,p1,pf;

let p1.ei= 0, p1.ef= 0, p1.pf= m**2+ki.kf, p1.ki= m*k,p1.kf=
    m*kp, pf.ei= -kf.ei, pf.ef= ki.ef, pf.ki= m*kp, pf.kf=    
    m*k, ki.ei= 0, ki.kf= m*(k-kp), kf.ef= 0, ei.ei= -1, ef.ef=
    -1; 

operator gp;

for all p let gp(p)= g(l,p)+m;

comment this is just to save us a lot of writing;

gp(pf)*(g(l,ef,ei,ki)/(2*ki.p1) + g(l,ei,ef,kf)/(2*kf.p1))
  * gp(p1)*(g(l,ki,ei,ef)/(2*ki.p1) + g(l,kf,ef,ei)/(2*kf.p1))$

write "The Compton cross-section is ",ws;

comment end of first physics example; 

off div;

comment another physics example;

index ix,iy,iz;

mass p1=mm,p2=mm,p3= mm,p4= mm,k1=0;

mshell p1,p2,p3,p4,k1;

vector qi,q2;

factor mm,p1.p3;

operator ga,gb;

for all p let ga(p)=g(la,p)+mm, gb(p)= g(lb,p)+mm; 

ga(-p2)*g(la,ix)*ga(-p4)*g(la,iy)* (gb(p3)*g(lb,ix)*gb(qi)
    *g(lb,iz)*gb(p1)*g(lb,iy)*gb(q2)*g(lb,iz)   +   gb(p3)
    *g(lb,iz)*gb(q2)*g(lb,ix)*gb(p1)*g(lb,iz)*gb(qi)*g(lb,iy))$

let qi=p1-k1, q2=p3+k1;

comment it is usually faster to make such substitutions after all the
        trace algebra is done;

write "CXN =",ws;

comment end of second physics example; 

showtime;

end;