Your MPW folder must be at the top level on the desk top, and must contain a folder called "fermat," which in turn must contain pow2028, all the errs files, the executable versions of Fermat, ferstartup, and the folder called ''object files", which contains all the Fermat object files. (You need the object files only if you are going to invoke Pascal or C programs, as described in "The Interface with C and Pascal.")
Open your UserStartup file and add the commands:
Alias ferm "{Boot}MPW:fermat:ferm"
Export ferm
for every version of Fermat (here called "ferm") that you have in the fermat folder. This makes it possible to invoke Fermat from inside any folder; i.e., no matter what the default directory is.
On 68K machines, Fermat assumes 68020 or better processor with 68881 or better coprocessor. On a PPC you must have a software fpu to simuate the 68881 calls.
Fermat needs a lot of stack space to run. Use Resedit to reset the MPW Shell resource HEXA 128 (or just plain HEXA) to 00064000 (00032000 in MPW 3.1 or earlier).
Fermat needs at least five meg of main memory, and you will need more if you do lots of graphics, especially movies.
Click on the MPW shell, then do command-I (information). Check the amount of space allocated to MPW. You want it to be as large as feasible.
Under system 6 (is anyone still using system 6?), make sure you have turned on the RAM cache with the control panel. 128K is enough.
With system 6, multifinder is a big memory hog. But if you have the memory, go ahead. Fermat works fine under systems 7, 8, and 9 if you've got at least 8 meg.
If you have a MacIIci (that was a great machine!), try turning off the color on your monitor and selecting only two gray scales with the control panel. This should provide a 30% speedup of all MPW tools, because there is no video card in the standard IIci. This is worth a try on any other old MacII as well.
Things will look best on the Worksheet if you choose the font Courier
12.
(includes both QFermat and FFermat)
Some of the symbols (glyphs) don't look quite right here, including
option-t, option-=, and option-shift-7.
| Symbol | Equivalent | Keystroke | Meaning | Use | Apply to arrays? |
| \ | - | \ | divide and truncate | x \y, [x]\y | yes |
| | | - | shift-\ | modulo | x | y, [x] |y | yes |
| |...| | - | shift-\ | absolute value | |(x+y)| | no |
| ? | - | option-= | not equal to | x ? y | no |
| > | - | option- > | greater than or equal | x > y | no |
| < | - | option- < | less than or equal | x < y | no |
| $ | - | shift-4 | greatest integer | $x | yes |
| ^ | - | shift-6 | exp(x) | ^x | no |
| ^ | - | shift-6 | raise to power | x^n | yes |
| _ | Ln | shift- - | natural log | _ x | no |
| : | - | shift- ; | suppress display; terminal only | 10^9999 : | no |
| ' | - | ' | start, end quote | ' ... ' | no |
| ¢ | - | option-shift-] | derivative of poly. or quoly | x¢ | no |
| ` | Sin | ` | sine | `x | no |
| " | Cos | shift-' | cosine | "x | no |
| % | Tan | shift-5 | tangent | %x | no |
| b | Arctan | option-s | arctangent | bx | no |
| p | - | option-p | pi = 3.14159... | p | no |
| p | - | option-p | pop | &p | no |
| P | - | option-shift-p | push | &P | no |
| P | - | option-shift-p | multiply out | P{i = 1,n } [...] | no |
| ÷ | Sqrt | option-v | square root | ÷x | no |
| ! | - | shift-1 | factorial | x! | no |
| ! | - | shift-1 | display | !(x,y,z), !!x, ![z | yes |
| ? | - | shift-/ | get from terminal | ?s, ?[x] | yes |
| © | - | option-g | complex conjugate | ©x | yes |
| ® | - | option-r | real part | ®x | yes |
| Ï | - | option-shift-f | imaginary part | Ï x | yes |
| ¡ | - | option-1 | ÷(-1) | ¡ | no |
| ¿ | - | option-shift-/ | random number | ¿, ¿¿ | yes |
| f | - | option-t | previous value | x := f, [f], f[n] | yes |
| S | - | option-w | add up | S{i = 1,n } [...] | no |
| Symbol | Equivalent | Keystroke | Meaning | Use | Apply to arrays? |
| · | - | option-8 | polynomial evaluation | x ·y , x·(u = y), x ·[y] | yes |
| · | - | option-8 | panic stop | &· | no |
| c | - | option-c | coefficient of poly. | c(x,n) | no |
| £ | Lterm | option-3 | leading term of poly. | £x | no |
| ¦ | Deg | option-d | degree of a polynomial or size of matrix | ¦x, ¦[x] | yes |
| & oe | - | option-q | monomial list (QFermat) | & oe | no |
| D | Det | option-j | determinant | D[x], others | yes |
| ¬ | - | option-shift-v | diagonal of matrix | ¬[x] | yes |
| TM | - | option-2 | transpose | TM[x] | yes |
| ¢ | Cols | option-4 | num. of columns | ¢ [x] | yes |
| f | - | option-f | concatenate arrays | [x] f [y] | yes |
| ± | - | option-shift-= | change Laurent | &± | no |
| @ | - | shift-2 | cancel | @[x], @F, @ < [x] > | yes |
| ¶ | - | option-7 | add (drop) poly. var. | &¶ | no |
| << | - | option-\ | start comment | << ... >> | no |
| >> | - | option-shift-\ | end comment | << ... >> | no |
| >> | - | option-shift-\ | exit function | & >> | no |
| > | - | shift-. | exit loop | & > | no |
| ] | - | ] | cycle (in loop) | &] | no |
| ... | - | option-; | ellipsis | x[2...,4...] | yes |
| ^ | - | shift-6 | boolean and | { i > 3 } ^ {j < 4 } | no |
| ò | - | option-shift- > | boolean or | { i > 3 } ò { x < 4 } | no |
| ÿ | - | option-l | boolean not | ÿ{ i > 3 } | no |
| ´ | - | option-5 | rational mode | &(p=´) | no |
| ,, | - | option-shift-w | line continuation | only from terminal | no |
| § | - | option-6 | system function | only from terminal | no |
| - | option-shift-7 | array of arrays | yes | ||
| Ø | - | option-shift-o | object parameter
(FFermat) |
F(Øx) = ... | yes |
| o | Codeg | option-k | codegree; in poly. | o x, o (x, 2) | no |
| Ç | - | option-shift-c | binomial coefficient | Ç(n, r) | no |
| ¥ | Var | option-y | poly var, by ordinal | ¥ 1 (= highest level poly var) | no |
Appendix 3. Some Significant Examples
These examples are not necessarily "best possible" solutions.
For typesetting reasons, all the · symbols
have been replaced by #.
Example 1.
;;** Miriam Duhan
;;** Numerical Analysis. Fordham University
;;** Final Project. May 14, 1990
;; Numerical analysis demonstration. Definite Integration.
;; How to use these functions:
;; (1) Use CompLeg to compute the Legendre polynomials and their derivatives.
;; Store in arrays leg and legp.
;; (2) Use Findrts and Findwts to fill the roots and weights arrays.
Write your own little program to do this.
;; (3) Then Quadra can be used to compute integrals. Make up your favorite
;; function F to be integrated.
;; (4) Compare for speed and accuracy to Romberg Integration.
;;************************************************************************
;;** QUADRATURE
;;** This uses formula 5.3.27 in Atkinson (An Introduction to
;;** Numerical Analysis, Wiley, 1989) pg 276 to do
;;** GAUSS-LEGENDRE QUADRATURE
;;**
;;** Input: eps the degree of accuracy required. This is not implemented
;;** since the best degree of accuracy I could get was n=2^5.
;;** n, lg of degree desired
;;** ll, lower limit of interval of integration
;;** up, upper limit of interval of integration
;;** Globals: F, the function is defined globally (by the user).
;;** arrays rt and wt already created and defined.
;;** Output: the value of the integral of the function on the interval
;;**
;;** The computation of 'newx' (and the multiplication of 'quad' by the
;;** factor 'z') adjusts the computation to the interval [-1,1].
Function Quadra(eps,n,ll,up,i,quad,maxi) =
if n > 4 then
!!('n = ', n:1, ' is out of bounds');
!!'Try again with n <= 4. Quadrature is done on 2^n points';
& >>
fi;
quad := 0;
z := up - ll;
row := n - 1;
maxi := 2^(n - 1) - 1;
for i = 0, maxi do
newx := (ll + up + rt[row,i]*z)/2;
quad := quad + wt[row,i]*F(newx)
od;
for i = 0, maxi do
newx := (ll + up - rt[row,i]*z)/2;
quad := quad + wt[row,i]*F(newx)
od;
quad := z*quad/2;
!!('for 2^n = ', 2^n:1, 'points');
!!' ';
!!(' ', up:1);
!!(' Integral F(x) dx = ', quad);
!!(ll:1).;
;;** This is the ROMBERG method of quadrature. I followed the algorithm
;;** in Atkinson on pg. 298 called 'Romberg'.
;;** Input: eps the degree of accuracy required.
;;** m, lg of degree desired
;;** ll, lower limit of interval of integration
;;** up, upper limit of interval of integration
;;** arrays a, b, t, already exist but are empty.
;;** Global: F the function is defined globally (by the user).
;;** Output: the value of the integral of the function on the interval.
Function Romberg(eps,maxk,ll,up,h,k,n,m,j,sum) =
n := 1;
t[0] := ((up - ll)/2)(F(ll) + F(up));
r[0] := t[0];
a[0] := t[0];
for k = 1, maxk - 1 do
n := 2n;
h := (up - ll)/n;
sum := 0;
for j = 1, n/2 do
sum := sum + F(ll + (2j - 1)h)
od;
t[k] := h*sum + t[k-1]/2;
for j = 0, k - 1 do
b[j] := a[j]
od;
a[0] := t[k];
m := 1;
for j = 1, k do
m := 4m;
a[j] := (m*a[(j-1)] - b[(j-1)])/(m - 1)
od;
r[k] := a[k];
if |(r[k] - r[(k-1)])| <= eps then
& >
fi
od;
if k > maxk - 1 then
!!'no convergence'
fi;
!!('k is', k:1);
!!(' ', up:1);
!!(' Integral F(x) dx = ', r[k-1]);
!!(ll:1).;
;;** This calculates the Legendre polynomials using the recursion formula:
;;** Pn+1(x) =(2n+1)/(n+1) x Pn - n/(n+1) Pn-1(x) Atkinson, pg 215 4.4.26
;;** It also calculates the derivatives of the polynomials, using the formula:
;;** P'n+1(x) = (2n+1)/(n+1) (Pn + x P'n) - n/(n+1) P'n-1(x)
;;**
;;** It saves the polynomials of degrees 2^n and 2^n + 1,
;;** for n = 1...m, m a parameter
;;** It saves the derivatives of the polys of degree 2^n
;;** (These three are needed to calculate roots and weights)
;;** p0,p1,p2 and pp0,pp1 and pp2 are three generations of polys and derivatives
;;** [leg] is an externally defined array which holds polys of degree 2^n
;;** [legp] is an externally defined array which holds derivatives
;;** [nextleg] is an externally defined array which holds polys of deg 2^n + 1
;;** Input: m the degree of the largest Legendre Polynomial required.
Function CompLeg(m,i,j,n,p0,p1,p2,pp0,pp1,pp2) =
p0 := 1;
p1 := x;
leg[0] := p1;
pp0 := 0;
pp1 := 1;
legp[0] := pp1;
n := 1;
for i = 1, m do
p2 := ((2i + 1)/(i + 1))*x*p1 - (i/(i + 1))*p0;
pp2 := ((2i + 1)/(i + 1))(p1 + x*pp1) - (i/(i + 1))*pp0;
<< save pprimen and pn >>
if i = 2^n - 1 then
legp[n] := pp2;
leg[n] := p2
fi;
<< save pn+1 >>
if i = 2^n then
nextleg[n] := p2;
n :+
fi;
p0 := p1;
p1 := p2;
pp0 := pp1;
pp1 := pp2
od;
!'end'.;
;;****************************************************************************
;;** ROOT FINDING
;;**
;;** Input eps degree of accuracy desired
;;** n, log of degree of poly (base 2)
;;** p2, polynomial of degree 2^n
;;** pp2, derivative of polynomial
;;**
;;** I ran into some problems with rootfinding.
;;** Using Newton's algorithm was fine. The problem was the initial guesses
;;** For polys of higher degrees (deg>8, not very high). It was not
;;** at all obvious how to divide the interval to make initial guesses.
;;** I started with the obvious one. Dividing the interval into 2^n-1 pieces.
;;**
;;** For n>3, this did not compute 2^n distinct roots.
Function Findrts(eps,n,p2,pp2,rt,j,nextrt,k) =
k := 2^(n - 1);
for j = 1, k do
nextrt := j/(k + 1);
rt := 0;
while |(rt - nextrt)| > eps do
rt := nextrt;
nextrt := rt - p2#rt/pp2#rt << # means poly evaluation >>
od;
rt[n-1,j-1] := nextrt
od;
!!('P(', 2^n:1, ')').;
;;** This worked to find distinct roots of the Legendre of degree 16,
;;** but it only found 30 distinct roots of the 32 degree Legendre.
Function Findrt2(eps,n,p2,pp2,rt,j,nextrt,k) =
k := 2^n;
rt := 0;
for j = 1, k do
nextrt := j/(k + 1);
while |(rt - nextrt)| > eps do
rt := nextrt;
nextrt := rt - p2#rt/pp2#rt << # means poly evaluation >>
od;
rt[n-1,j-1] := nextrt
od;
!!('P(', 2^n:1, ')').;
;;************************************************************************
;;** Weight finding.
;;** This function finds the weights associated with a given degree
;;** Legendre poly. It uses formula 5.3.28, pg 276 in Atkinson.
;;** Input: the Legendre Poly of degree n+1, the derivative of the
;;** Legendre poly of degree n.
;;** Output: the n-1st row of a table of weights corresponding to polys
;;** of degree n.
Function Findwts(pp1,p2,n,j,i) =
i := 2^n;
for j = 1, i/2 do
wt[n-1,j-1] := -2/((i + 1)*(pp1#rt[n-1,j-1])*(p2#rt[n-1,j-1]))
od.;
&(p = 1); ;; Move into fast float mode.
&(a = 0); ;; Set array initial index to 0.
&(¶ = x);
Array legp[10]; Array nextleg[10]; Array rt[6,16]; Array a[10]; Array b[10]; Array leg[10]; Array wt[6,16]; Array r[10]; Array t[10]; &x;
Example 2.
Included here are
First, the C program, in a file called ferc.c:
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
##pragma segment [segext]
Cfunc(m) { This name MUST be Cfunc. }
int *m;
{
int n;
char infname[12] = "instuff";
char outfname[12] = "outstuff";
FILE *infile;
FILE *outfile;
infile = fopen(infname, "r");
outfile = fopen(outfname, "a");
fscanf(infile, "%d", &n);
fprintf(outfile,":x = %d %d \n", n, m[0]);
m[0] = n+1;
fclose(outfile);
fclose(infile);
}
Next, the Pascal program, in a file called ferp.p:
unit ferp_ext;
interface
{$S[segext]}
{$R-}
uses paslibintf, sane;
procedure Pfunc(var m: longint); { This name MUST be Pfunc. }
implementation
procedure Pfunc;
var n: integer;
infname, outfname: text;
begin
open(infname, 'instuff');
open(outfname, 'outstuff');
rewrite(outfname);
readln(infname, n);
writeln(outfname, ':x = ', n, m); { a Fermat readable statement }
m := n + 1; { pass back n+1 to Fermat }
close(infname);
close(outfname)
end;
end.
Next, the make file, called ferext.make,
that compiles the above programs and links them with the Fermat object
files. This assumes that the default directory is
::mpw:fermat:,
the mpw folder is at the top level on the desk top, the mpw folder contains
all the right interface and library files, and that the object files for
the above Pascal and C programs are in the top level inside the fermat
folder:
# File: ferext.make
# Target: ferext
# Sources: ferc.c ferp.p
ferp.p.o Ú- ferext.make ferp.p
Pascal -MC68020 -MC68881 -D elems881=true -p -r ferp.p
ferc.c.o Ú- ferext.make ferc.c
C -MC68020 -MC68881 -D elems881=true -p -r ferc.c
SOURCES = ferc.c ferp.p
OBJECTS = ¶
::Libraries:Libraries:Interface.o ¶
::Libraries:Libraries:Runtime.o ¶
::Libraries:Libraries:ToolLibs.o ¶
::PLibraries:PasLib.o ¶
::PLibraries:SANELib881.o ¶
":object files:ferm114.p.o" ¶
":object files:fer414.p.o" ¶
":object files:fer314.p.o" ¶
":object files:fer214.p.o" ¶
ferp.p.o ¶
ferc.c.o
ferext Ú-Ú- ferext.make OBJECTS
Link -ss 200000 -p -w -t MPST -c 'MPS ' ¶
OBJECTS ¶
"{Libraries}" stubs.o ¶
"{CLibraries}"CRuntime.o ¶
"{Libraries}"Runtime.o ¶
"{CLibraries}"StdCLib.o ¶
"{CLibraries}"CSANELib.o ¶
"{CLibraries}"Math.o ¶
"{CLibraries}"CInterface.o ¶
"{CLibraries}"Complex.o ¶
"{Libraries}" Interface.o ¶
"{PLibraries}"PasLib.o ¶
"{PLibraries}"SANELib.o ¶
"{Libraries}"ToolLibs.o ¶
-o ferext
Example 3.
Included here are several functions that invoke graphics routines.
pi := p;
<< This is another typesetting-font cludge; not necessary really.
Just
use the symbol p
in the functions below. >>
Function Globe(r,inc,x0,y0,z0,x1,y1,z1,s,t,u,doit) = << Draw a three-dimensional stick figure globe centered at the origin. r is the radius. r, t, and s are the spherical coordinates of the successive vertices, t = polar angle. inc gives the angle (in radians) between successive points. Invoke this function with Globe(<radius>, <inc>). >> t := 0; doit := 1; << Boolean true >> while t <= pi do s := 0; x0 := r*Sin(t); << Cos(s) = 1 >> y0 := 0; << Sin(s) = 0 >> z0 := r*Cos(t); while s <= 2pi do s :+ inc; x1 := r*Sin(t)Cos(s); y1 := r*Sin(t)Sin(s); z1 := r*Cos(t); << Draw 'horizontal' latitude line. >> Line(x0, y0, z0, x1, y1, z1); if t = 0 then & > fi; if doit then u := t - inc; << previous t level >> x0 := r*Sin(u)*Cos(s); y0 := r*Sin(u)*Sin(s); z0 := r*Cos(u); doit := 0; << Draw 'vertical' longitude line back to previous t level. >> Line(x0, y0, z0, x1, y1, z1) else doit := 1 fi; x0 := x1; y0 := y1; z0 := z1 od; t :+ inc od.; Function Turn(n,j,i) = for i = 1, 2*pi*n do Copy(x, z[i]); ERotate(x, 1/n, j) od.;Use the above two functions in the following way:
Erase;
Setthree(300, 250);
Setscale(30, 30, 30)
Mark;
Globe(4, p/10);
Setobj(x);
Objarray(z[200]);
Turn(30, 3);
Animate(z);
Demonstration of complex function graphing:
Func F(x,y) = y * ^(2pix)/10.;
Func G(x) = Sin(x).; Function P(n,z,s) = << Plot the absolute value of the complex valued function Sin(x) on ever widening circles. n is the number of circles, z the increment, s the flag to select points or line segments. >> for i = 1, n do if i = 20 then z := z/2 fi; Cfuncplot(s, 0, 1, z, F(..., i), G) od.;
Planetary simulation. These functions and commands create objects to simulate four of the planets. The radii and periods of the orbits are not authentic.
pi := p;
<< This is another typesetting-font cludge; not necessary really.
Just
use the symbol p
in the functions below. >>
Function Setup = Erase; Setthree(300, 250); Setscale(35, 35, 35).; Function Run(n,i) = Penstyle(4, 4, 'O'); Plot(0, 0, 0); << the Sun >> for i = 1, n do Penstyle(4, 2); ERotate(mars, pi/100, 3); Penstyle(4, 4, 'J'); ERotate(jup, pi/200, 3); Penstyle(4, 2); ERotate(earth, pi/60, 3); Penstyle(4, 4, 'S'); ERotate(saturn, pi/350, 3) od; Penstyle(4, 0).; Function Grid(n,i) = Penstyle(1, 0); for i = -n, n do Line(i, n, 0, i, -n, 0); Line(n, i, 0, -n, i, 0) od; Penstyle(4, 0).; Setup; Grid(7); Mark; Plot(1,0,0); Setobj(earth); Mark; Plot(2,0,0); Setobj(mars); Mark; Plot(4,0,0); Setobj(jup); Mark; Plot(6,0,0); Setobj(saturn); Run(2000);
An example of a hidden line object:
[x] := [(-3, 3, -3, 3)];
[y] := [(-3, -3, 3, 3)];
Function F2(x,y,i,s,u) = << Compute the "distance" of (x, y) to each of the 4 points given by the arrays [x] and [y]. Add up the reciprocals. The 3.42 below was computed in advance to be right for squaring off the tops of the peaks in the picture. >> Compile([x], [y]); Real; for i from 1 to 4 do u := |(x - x[i])|; u := u + |(y - y[i])|; if u = 0 then s := 3.42; & > fi; s := s + 1/u od; if s > 3.42 then 25 else 25*s/3.42 fi.;The following commands invoking F2 produce the surface below:
Erase;
Setthree(300, 250);
Setscale(30, 30, 30)
Surface(x, 0, -5, 5, -5, 5, 1/3, F2);
ERotate(x, p/8,
3);
View(p/6,
x);
This group of functions traces a trefoil knot. Preface with Mark.
The resulting object (rotated) is below.
pi := p;
<< This is another typesetting-font cludge; not necessary really.
Just
use the symbol p
in the functions below. >>
Function Cross(x,y,z) =
<< Cross product of vectors. >>
Real;
z[1] := x[2]*y[3] - y[2]*x[3];
z[2] := x[3]*y[1] - y[3]*x[1];
z[3] := x[1]*y[2] - y[1]*x[2].;
Function Turn_into(x, y, rot) =
<< Two vectors are given, [x] and [y] (anchored at origin) with
[x] x [y] <> 0. Set [rot] = the rotation that carries [x] onto
[y]'s direction. [x] and [y] are made unit vectors. >>
Real;
Array sec[3,3];
Array z[3];
Array u[3];
Array v[3];
Cross([x], [y], [z]);
if Iszero[z] then
!!('bad data');
Return(0)
fi;
[x] := [x]/Mod([x]);
[y] := [y]/Mod([y]);
[z] := [z]/Mod([z]);
Cross([z], [x], [u]);
Cross([z], [y], [v]);
[sec[,1]] := [x];
[sec[,2]] := [z];
[sec[,3]] := [u];
[rot[,1]] := [y];
[rot[,2]] := [z];
[rot[,3]] := [v];
[rot] := [rot] * [sec]^-1;
@([v], [u], [z], [sec]).;
Function Path(t,x; u) =
<< Used by Thick below. Parameterizes a trefoil knot in three-space.
Note that x is an array. >>
Real;
t := 2*pi*t;
u := 2*(1 + 0.3*Cos(3*t));
x[1] := u*Cos(2*t);
x[2] := u*Sin(2*t);
x[3] := u*0.35*Sin(3*t).;
Function Mod(x) =
Real;
Sqrt((x[1]^2 + x[2]^2 + x[3]^2)).;
Function Color(x) =
Real;
4*x*(1 - x).;
Function Makebases =
Real;
[n0] := [x1] - [x0];
[n1] := [x2] - [x1];
Cross([n0], [n1], [b1]);
Cross([b1], [n0], [b2]);
[b1] := [b1]/Mod([b1]);
[b2] := [b2]/Mod([b2]).;
Function Thick(rad,n,m; i,j,k,alpha) =
<< A path in three-space is given by the function Path. A small circle
(actually a polygon) moves perpendicularly along it, "thickening" the path.
rad = radius of the small circle, n = number of copies of the circle
generated, m = number of sides to the polygon approximating the circle.
[v] is array of 12 items, i.e. 4 points.
Example: Invoke with Thick(2/5, 120, 16) and use Path above. >>
Real;
Array pts[3,m+1];
Array nxt[3,m+1];
alpha := 2*pi/m;
[pts] := 0;
Array x0[3];
Array x1[3];
Array x2[3];
Array b1[3];
Array b2[3];
Path(-2/n, [x0]);
Path(-1/n, [x1]);
Path(0, [x2]);
<< Create bases [b1] and [b2] orthogonal to [x1] - [x0]. >>
Makebases;
<< Create initial circle of points, centered at [x1], on plane orthogonal
to [x1]-[x0]. >>
for i = 0, m do
[pts[,i+1]] := [x1] + rad*Cos(alpha*i)*[b1] + rad*Sin(alpha*i)*[b2]
od;
<< "Translate" circle n times, following the path. Can reuse [b1] and
[b2]. >>
Array rot[3,3];
for i = 1, n do
<< At start of loop, [pts] are a circle in the plane orthog to [x1]-[x0]
centered at [x1], and [x2] has been computed. >>
[b2] := [x2] - [x1];
[b1] := [x1] - [x0];
<< Bring the circle down to the origin, then rotate it to become the
next circle and move it up to [x2]. >>
for j = 0, m do
[nxt[,j+1]] := [pts[,j+1]] - [x1]
od;
Turn_into([b1], [b2], [rot]);
for j = 0, m do
[nxt[,j+1]] := [rot]*[nxt[,j+1]] + [x2]
od;
SColor(60000*Color(i/n), 20000, 60000*(1 - Color(i/n)));
for j = 1, m do
[v[1...3]] := [pts[,j]];
[v[4...6]] := [pts[,j+1]];
[v[7...9]] := [nxt[,j+1]];
[v[10...]] := [nxt[,j]];
Polygon([v], -1)
od;
<< Prepare for next pass. >>
[x0] := [x1];
[x1] := [x2];
Path(i/n, [x2]);
[pts] := [nxt]
od.;
This object was created with the functions above.
Finally, here is a function that creates a torus. Preface with Mark. Below are two toruses created this way, then rotated, translated, shaded, colored, and merged into one object with Union.
Function Torus(r1,r2,n,m,i,j,h,alpha,delta,rad) = << Create torus centered at origin parallel to x-y plane. r1 = outer radius, r2 = inner radius, n and m = number of pieces the two perpendicular circles are cut into. For example, invoke with Torus(4,1.3,32,16) >> Real; Array pts[3,m+1]; Array nxt[3,m+1]; Array rotate[3,3]; delta := 2*pi/n; alpha := 2*pi/m; h := (r1 + r2)/2; rad := (r1 - r2)/2; [rotate] := [1]; rotate[1,1] := Cos(delta); rotate[1,2] := -Sin(delta); rotate[2,2] := Cos(delta); rotate[2,1] := Sin(delta); [pts] := 0; << Create initial circle of points. >> for i = 0, m do pts[1,i+1] := - h - rad*Cos(alpha*i); pts[3,i+1] := rad*Sin(alpha*i) od; << Rotate circle through angle delta n times. >> for i = 1, n do [nxt] := [rotate]*[pts]; for j = 1, m do [v[1...3]] := [nxt[,j]]; [v[4...6]] := [nxt[,j+1]]; [v[7...9]] := [pts[,j+1]]; [v[10...]] := [pts[,j]]; << Color each polygon red. >> SColor(60000,0,0); Polygon[v] od; [pts] := [nxt] od.;
Two toruses created with the functions above. On a color monitor, these are shades of purple and green.