PROC tc = (PROC(REAL)REAL sin, cos)VOID: (
PROC tz = (REAL tr, r)VOID: (
PROC sq = (REAL x)REAL: x*x;
REAL f, fp, fpp;
(REAL sl = cos(tr), cn = sin(tr) - tr*cos(tr);
PROC ff = (REAL r)REAL: sq(IF r < tr THEN sin(r) ELSE cn + sl*r FI),
    ffp = (REAL r)REAL: IF r < tr THEN sin(2*r)
                        ELSE 2*(cn*sl + sq(sl)*r) FI,
    ffpp= (REAL r)REAL: IF r < tr THEN 2*cos(2*r) ELSE 2*sl*sl FI;
f := ff(r); fp := ffp(r); fpp := ffpp(r);
PROC td = (PROC(REAL)REAL g, REAL pd, STRING s)VOID:
  IF ABS(g(r+0.00001)-g(r))/0.00001 - pd> 0.001
  THEN print(("Foul ", s, g(r+0.00001)-g(r), pd, newline)) FI;
td(ff, fp, "ff"); td(ffp, fpp, "ffp"));
print(("tr = ", tr, ", r = ", r, newline));

PROC sum = (PROC(INT)REAL f)REAL: (REAL s := 0;
   FOR j FROM 0 TO 3 DO s +:= f(j) OD; s);
PROC delgam = (INT i, j, k, l)REAL:
  IF i=1 AND j=1 AND k=2 AND l=2 THEN -0.5 * fpp
  ELIF i=1 AND j=2 AND (k=1 AND l=2 OR k=2 AND l=1)
     THEN 0.5*(fpp/f - sq(fp/f))
  ELSE 0 FI;
[0:3, 0:3, 0:3] REAL gam;
[0:3, 0:3, 0:3, 0:3] REAL rie;
[0:3, 0:3] REAL gl, gu;
FOR i FROM 0 TO 3 DO FOR j FROM 0 TO 3 DO gl[i,j]:=0; gu[i,j]:=0 OD OD;
  gl[0,0] := -1; gl[1,1] := 1; gl[2,2] := f; gl[3,3] := 1;
  gu[0,0] := -1; gu[1,1] := 1; gu[2,2] := 1.0/f; gu[3,3] := 1;
FOR i FROM 0 TO 3 DO FOR j FROM 0 TO 3 DO
FOR k FROM 0 TO 3 DO gam[i,j,k] :=
  IF i=1 AND j=2 AND k=2 THEN -0.5*fp
  ELIF i=2 AND (j=2 AND k=1 OR j=1 AND k=2) THEN 0.5*fp/f
  ELSE 0 FI OD OD OD;
FOR i FROM 0 TO 3 DO FOR j FROM 0 TO 3 DO
FOR k FROM 0 TO 3 DO FOR l FROM 0 TO 3 DO
   rie[i,j,k,l] := delgam(k, i, l, j) - delgam(l, i, k, j)
     + sum((INT n)REAL: gam[i,k,n]*gam[n,l,j] - gam[i,l,n]*gam[n,k,j])
   OD OD OD OD;
FOR i FROM 0 TO 3 DO FOR j FROM 0 TO 3 DO
FOR k FROM 0 TO 3 DO FOR l FROM 0 TO 3 DO
  REAL rr = sum((INT n)REAL: gl[i,n]*rie[n,j,k,l]);
  IF rr /= 0 THEN print(("R ", i, j, k, l, "= ", rr, newline))
  FI OD OD OD OD;
[0:3, 0:3] REAL rt; # Ricci tensor #
FOR i FROM 0 TO 3 DO FOR j FROM 0 TO 3 DO
  rt[i,j] := sum((INT k)REAL: rie[k, i, k, j]) OD OD;
FOR i FROM 0 TO 3 DO FOR j FROM 0 TO 3 DO
  IF rt[i, j] /= 0 THEN print(("Ricci[", i, ", ", j, "] = ",
     rt[i, j], newline)) FI OD OD;
REAL bigr = sum((INT i)REAL: sum((INT j)REAL: gu[i,j]*rt[j,i]));
print(("R = ", bigr, newline));
[0:3, 0:3] REAL bigg; # Einstein tensor #
FOR i FROM 0 TO 3 DO FOR j FROM 0 TO 3 DO
  bigg[i,j] := rt[i,j] - 0.5*bigr*gl[i,j] OD OD;
FOR i FROM 0 TO 3 DO FOR j FROM 0 TO 3 DO
  IF bigg[i,j] /= 0 THEN
    print(("G[", i, ",", j, "] = ", bigg[i,j], newline)) FI OD OD);
print(("Flat space", newline));
tz(0, 1.0); tz(0, 1.3); # The flat plane in polar coordinates. #
print(("flat part of curved space", newline));
tz(1, 1.4); # A conical spot outside the string #
print(("Curved part", newline));
tz(1, 0.5); tz(1, 0.6); tz(0.7, 0.31) # three places within the string .#
);
print(("Positive", newline)); tc(sin, cos);
print(("Negative", newline)); tc(sinh, cosh)