3D.MP : 3-DIMENSIONAL REPRESENTATIONS AND ANIMATIONS IN METAPOST

version 1.0, 8 April 1998 Denis Roegel (roegel@loria.fr) This package provides definitions enabling the manipulation and animation of 3-dimensional objects. Such objects can be included in a TEX file or used on web pages for instance. See the documentation enclosed in the distribution for more details. Thanks to John Hobby and Ulrik Vieth for helpful hints. LIMITATIONS: − an object can not include TEX labels; to overcome this limitation would need to wrap the metapost output with dvips so that the necessary fonts are loaded (easy to do, but not done here) PROJECTS FOR THE FUTURE: − take light sources into account and show shadows and darker faces − handle overlapping of objects (obj name can be used when going through all faces)

if known three d version:

expandafter endinput% avoids loading this package twice fi; message ; numeric three d version; three d version = 1.0; Vector operations

% definition of vector i by its coordinates def vect def (expr i, xi, yi, zi) = vect[i]x := xi; vect[i]y := yi; vect[i]z := zi; enddef;

% a point is stored as a vector let set point = vect def ;

% set obj point is like set point, but the first parameter % is a local point number def set obj point(expr i, xi, yi, zi) = set point(pnt(i), xi, yi, zi) enddef;

% definition of a vector by an other vector def vect def vect(expr i, j) =

vect[i]x := vect[j]x; vect[i]y := vect[j]y; vect[i]z := vect[j]z; enddef;

% vector sum: vect[k] ← vect[i]+vect[j] def vect sum(expr k, i, j) =

vect[k]x := vect[i]x + vect[j]x; vect[k]y := vect[i]y + vect[j]y; vect[k]z := vect[i]z + vect[j]z;

enddef;

% vector translation: vect[i] ← vect[i]+vect[v] def vect translate(expr i, v) = vect sum(i, i, v) enddef;

% vector difference: vect[k] ← vect[i]−vect[j] def vect diff (expr k, i, j) =

vect[k]x := vect[i]x − vect[j]x;

vect[k]y := vect[i]y − vect[j]y;
vect[k]z := vect[i]z − vect[j]z;

enddef;

% dot product of vect[i] and vect[j] vardef vect dprod(expr i, j) =

(vect[i]x ∗ vect[j]x + vect[i]y ∗ vect[j]y + vect[i]z ∗ vect[j]z) enddef;

% modulus of vect[i] vardef vect mod(expr i) = sqrt(vect dprod(i, i)) enddef;

% vector product: vect[k] ← vect[i] ∧ vect[j] def vect prod(expr k, i, j) =

vect[k]x := vect[i]y ∗ vect[j]z − vect[i]z ∗ vect[j]y;

vect[k]y := vect[i]z ∗ vect[j]x − vect[i]x ∗ vect[j]z;
vect[k]z := vect[i]x ∗ vect[j]y − vect[i]y ∗ vect[j]x;

enddef;

% scalar multiplication: vect[j] ← vect[i] ∗ v def vect mult(expr j, i, v) =

vect[j]x := v ∗ vect[i]x; vect[j]y := v ∗ vect[i]y; vect[j]z := v ∗ vect[i]z; enddef;

% middle of two points def mid point(expr k, i, j) = vect sum(k, i, j); vect mult(k, k, .5); enddef; Vector rotation

% Rotation of vect[v] around vect[axis] by an angle alpha

The vector v⃗ is first projected on the axis giving vectors a⃗ and h⃗:

If we set b⃗ = \frac{\overrightarrow{axis}}{∥\overrightarrow{axis}∥}, the rotated vector v⃗′ is equal to h⃗ + f⃗ where f⃗ = cos α · a⃗ + sin α · c⃗. and h⃗ = (v⃗ · b⃗)b⃗

% The rotation is independent of vect[axis]’s module. % v = old and new vector % axis = rotation axis % alpha = rotation angle % vardef vect rotate(expr v, axis, alpha) =

save v a, v b, v c, v d, v e, v f ; v a := new vect; v b := new vect; v c := new vect; v d := new vect; v e := new vect; v f := new vect; v g := new vect; v h := new vect; vect mult(v b, axis, 1/vect mod(axis)); vect mult(v h, v b, vect dprod(v b, v)); % projection of v on axis vect diff (v a, v, v h); vect prod(v c, v b, v a); vect mult(v d, v a, cosd(alpha)); vect mult(v e, v c, sind(alpha)); vect sum(v f , v d, v e); vect sum(v, v f , v h); free vect(v h); free vect(v g); free vect(v f ); free vect(v e); free vect(v d); free vect(v c); free vect(v b); free vect(v a);

enddef; Operations on objects

% iname is the handler for an instance of an object of class name % iname must be a letter string vardef assign obj (expr iname, name) =

save tmpdef ; string tmpdef ; % we need to add double quotes (char 34) tmpdef = & iname & & ditto & name & ditto & ; scantokens tmpdef ; def obj (iname);

enddef;

% name is the the name of an object instance % It must be made only of letters (or underscores), but no digits. def def obj (expr name) =

scantokens begingroup

save tmpdef ; string tmpdef ;

tmpdef = & obj class (name) & & ditto & name & ditto & ; tmpdef

endgroup enddef;

% This macro puts an object back where it was right at the beginning. % iname is the name of an object instance. vardef reset obj (expr iname) =

save tmpdef ; string tmpdef ; define current point offset (iname); tmpdef = & obj class (iname) & ; scantokens tmpdef (iname);

enddef;

% Put an object at position given by pos (a vector) and % with orientations given by angles psi, theta, phi. % The object is scaled by scale. % iname is the name of an object instance. vardef put obj (expr iname, pos, scale, psi, theta, phi) =

reset obj (iname); scale obj (iname, scale); v x := new vect; v y := new vect; v z := new vect; vect def vect(v x, vect I ); vect def vect(v y, vect J ); vect def vect(v z, vect K ); rotate obj abs pv(iname, point null, v z, psi); vect rotate(v x, v z, psi); vect rotate(v y, v z, psi); rotate obj abs pv(iname, point null, v y, theta); vect rotate(v x, v y, theta); vect rotate(v z, v y, theta); rotate obj abs pv(iname, point null, v x, phi); vect rotate(v y, v x, phi); vect rotate(v z, v x, phi); free vect(v z); free vect(v y); free vect(v x); translate obj (iname, pos);

enddef; Rotation, translation and scaling of objects

% Rotation of an object instance name around an axis % going through a point p (local to the object) % and directed by vector vect[v]. The angle of rotation is a. vardef rotate obj pv(expr name, p, v, a) =

define current point offset (name); rotate obj abs pv(name, pnt(p), v, a);

enddef;

vardef rotate obj abs pv(expr name, p, v, a) = save v a; define current point offset (name); v a := new vect; for i := 1 upto obj points (name):

vect diff (v a, pnt(i), p); vect rotate(v a, v, a); vect sum(pnt(i), v a, p);

endfor; free vect(v a);

enddef;

% Rotation of an object instance name around an axis % going through a point p (local to the object) % and directed by vector \overrightarrow{pq}. The angle of rotation is a. vardef rotate obj pp(expr name, p, q, a) =

save v a, axis; define current point offset (name); v a := new vect; axis := new vect; vect diff (axis, pnt(q), pnt(p)); for i := 1 upto obj points (name):

vect diff (v a, pnt(i), pnt(p)); vect rotate(v a, axis, a); vect sum(pnt(i), v a, pnt(p));

endfor; free vect(axis); free vect(v a);

enddef;

% Translation of an object instance name by a vector vect[v]. vardef translate obj (expr name, v) =

define current point offset (name); for i := 1 upto obj points (name):

vect sum(pnt(i), pnt(i), v); endfor;

enddef;

% Scalar multiplication of an object instance name by a scalar v. vardef scale obj (expr name, v) =

define current point offset (name); for i := 1 upto obj points (name):

vect mult(pnt(i), pnt(i), v); endfor;

enddef; Functions to build new points in space

% Rotation in a plane: this is useful to define a regular polygon. % k is a new point obtained from point j by rotation around o % by a angle α equal to the angle from i to j.

vardef rotate in plane (expr k, o, i, j) = save cosalpha, sinalpha, alpha, v a, v b, v c; v a := new vect; v b := new vect; v c := new vect; vect diff (v a, i, o); vect diff (v b, j, o); vect prod(v c, v a, v b); cosalpha = vect dprod(v a, v b)/vect mod(v a)/vect mod(v b); sinalpha = sqrt(1 − cosalpha ∗∗ 2); alpha = angle((cosalpha, sinalpha)); vect rotate(v b, v c, alpha); vect sum(k, o, v b); free vect(v c); free vect(v b); free vect(v a);

enddef;

vardef rotate in plane(expr k, o, i, j) =

rotate in plane (pnt(k), o, pnt(i), pnt(j)) enddef;

% Build a point on a adjacent face. The middle m of points i and j is such that ˆ (\overrightarrow{om}, \overrightarrow{mc}) = α This is useful to define regular polyhedra

vardef new face point (expr c, o, i, j, alpha) = save v a, v b, v c, v d, v e; v a := new vect; v b := new vect; v c := new vect; v d := new vect; v e := new vect; vect diff (v a, i, o); vect diff (v b, j, o); vect sum(v c, v a, v b); vect mult(v d, v c, .5); vect diff (v e, i, j); vect sum(c, v d, o); vect rotate(v d, v e, alpha); vect sum(c, v d, c); free vect(v e); free vect(v d); free vect(v c); free vect(v b); free vect(v a);

enddef;

vardef new face point(expr c, o, i, j, alpha) =

new face point (pnt(c), pnt(o), pnt(i), pnt(j), alpha) enddef;

vardef new abs face point(expr c, o, i, j, alpha) =

new face point (c, o, pnt(i), pnt(j), alpha) enddef; Computation of the projection of a point on the “screen”

% p is the projection of m % m = point in space (3 coordinates) % p = point of the intersection plane

vardef project point(expr p, m) = save tmpalpha, v a, v b; v a := new vect; v b := new vect; vect diff (v b, m, Obs); % vector Obs-m

% vect[v a] is vect[v b] expressed in (ObsI ,ObsJ ,ObsK )

% coordinates. vect[v a]x := vect[IObsI ]x ∗ vect[v b]x + vect[IObsJ ]x ∗ vect[v b]y + vect[IObsK ]x ∗ vect[v b]z; vect[v a]y := vect[IObsI ]y ∗ vect[v b]x + vect[IObsJ ]y ∗ vect[v b]y + vect[IObsK ]y ∗ vect[v b]z; vect[v a]z := vect[IObsI ]z ∗ vect[v b]x + vect[IObsJ ]z ∗ vect[v b]y + vect[IObsK ]z ∗ vect[v b]z; if vect[v a]x < Obs dist: % then, point m is too close

message & decimal m & ; x[p] := too big ; y[p] = too big ;

else:

if (angle(vect[v a]x, vect[v a]z) > h field/2)

or (angle(vect[v a]x, vect[v a]y) > v field/2):

message & decimal m & ; x[p] := too big ; y[p] = too big ;

else:

tmpalpha := Obs dist/vect[v a]x;

y[p] := drawing scale ∗ tmpalpha ∗ vect[v a]y; x[p] := drawing scale ∗ tmpalpha ∗ vect[v a]z;

fi; fi; free vect(v b); free vect(v a);

enddef;

% Object projection % This is a mere iteration on project point def project obj (expr name) =

define current point offset (name); for i := 1 upto obj points (name): project point(ipnt (i), pnt(i)); endfor;

enddef; Draw one face, hiding it if it is hidden

% The order of the vertices determines what is the visible side % of the face. The order must be clockwise when the face is seen. % drawhidden is a boolean; if true only hidden faces are drawn; if false, % only visible faces are drawn. Therefore, draw face is called twice % by draw faces. vardef draw face(text vertices)(expr col, drawhidden) =

save p, num, normal vect, v a, v b, v c, overflow; path p; boolean overflow; overflow = false; forsuffixes $ = vertices:

if z[ipnt ($)] = (too big , too big ): overflow := true; fi; exitif overflow;

endfor; if overflow: message ;

else:

p = forsuffixes $ = vertices: z[ipnt ($)] - - endfor cycle; num_{0} := 0; num_{1} := 0; num_{2} := 0;

% get the three last suffixes: forsuffixes $ = vertices: num_{0} := num_{1}; num_{1} := num_{2}; num_{2} := $; endfor; normal vect := new vect; v a := new vect; v b := new vect; v c := new vect; vect diff (v a, pnt(num_{1}), pnt(num_{0})); vect diff (v b, pnt(num_{2}), pnt(num_{1})); vect prod(normal vect, v a, v b); vect diff (v c, pnt(num_{1}), Obs);

if filled faces:

if vect dprod(normal vect, v c) < 0:

fill p withcolor col; if draw contours: drawcontour(p, contourwidth); fi; else: % draw p dashed evenly; if this is done, you must ensure

% that hidden faces are (re)drawn at the end

fi;

else:

if vect dprod(normal vect, v c) < 0: %visible if not drawhidden: draw p; fi;

else: % hidden

if drawhidden:

draw p withcolor background; % avoid strange overlapping dashes draw p dashed evenly;

fi; fi;

fi; free vect(v c); free vect(v b); free vect(v a); free vect(normal vect);

fi; enddef;

% p is the path to draw (a face contour) and thickness is the pen width def drawcontour(expr p, thickness) =

pickup pencircle scaled thickness; draw p; pickup pencircle scaled .4pt;

enddef; % Variables for face handling. First, we have an array for lists of vertices % corresponding to faces. string face points []; % analogous to vect arrays

% Then, we have an array of colors. A color needs to be a string % representing an hexadecimal RGB coding of a color. string face color [];

% name is the name of an object instance vardef draw faces(expr name) =

save tmpdef ; string tmpdef ; define current face offset (name);

% first the hidden faces (dashes must be drawn first): for i := 1 upto obj faces (name):

tmpdef := & face points [face(i)] & & ditto & face color [face(i)] & ditto & ; scantokens tmpdef ;

endfor;

% then, the visible faces: for i := 1 upto obj faces (name):

tmpdef := & face points [face(i)] & & ditto & face color [face(i)] & ditto & ; scantokens tmpdef ;

endfor; enddef;

% Draw point n of object instance name vardef draw point(expr name, n) =

define current point offset (name); project point(ipnt (n), pnt(n)); if z[ipnt (n)] ≠ (too big , too big ):

pickup pencircle scaled 5pt;

drawdot(z[ipnt (n)]);

pickup pencircle scaled .4pt; fi;

enddef;

vardef draw axes(expr r, g, b) = project point(1, vect null); project point(2, vect I ); project point(3, vect J ); project point(4, vect K ); if (z_{1} ≠ (too big , too big )):

if (z_{2} ≠ (too big , too big )):

drawarrow z_{1} - - z_{2} dashed evenly withcolor r; fi; if (z_{3} ≠ (too big , too big )):

drawarrow z_{1} - - z_{3} dashed evenly withcolor g; fi; if (z_{4} ≠ (too big , too big )):

drawarrow z_{1} - - z_{4} dashed evenly withcolor b; fi;

fi; enddef; % Definition of a macro obj name returning an object name % when given an absolute % face number. This definition is built incrementally through a string, % everytime a new object is defined. % obj name is defined by redefine obj name .

% Initial definition string index to name ; index to name = ;

% name is the name of an object instance % n is the absolute index of its last face def redefine obj name (expr name, n) =

index to name := index to name & & decimal n & & ditto

& name & ditto;

scantokens begingroup index to name & endgroup; enddef;

% i is an absolute face number % vertices is a string representing a list of vertices % rgbcolor is a string representing a color in rgb hexadecimal def set face(expr i, vertices, rgbcolor) =

face points [i] := vertices; face color [i] := rgbcolor; enddef;

% i is a local face number %vertices is a string representing a list of vertices %rgbcolor is a string representing a color in rgb hexadecimal def set obj face(expr i, vertices, rgbcolor) = set face(face(i), vertices, rgbcolor) enddef; Compute the vectors corresponding to the observer’s viewpoint

% (vectors ObsI ,ObsJ and ObsK in the vect I ,vect J , % vect K reference; and vectors IObsI ,IObsJ and IObsK % which are vect I ,vect J ,vect K % in the ObsI ,ObsJ ,ObsK reference)

(here, ψ > 0, θ < 0 and ϕ > 0; moreover, |θ| ≤ 90^{∘})

def compute reference(expr psi, theta, phi) = % ObsI defines the direction of observation; % ObsJ and ObsK the orientation % (but one of these two vectors is enough, % since ObsK = ObsI ∧ ObsJ ) % The vectors are found by rotations of vect I ,vect J ,vect K .

vect def vect(ObsI , vect I ); vect def vect(ObsJ , vect J ); vect def vect(ObsK , vect K ); vect rotate(ObsI , ObsK , psi); vect rotate(ObsJ , ObsK , psi); % gives (u,v,z) vect rotate(ObsI , ObsJ , theta); vect rotate(ObsK , ObsJ , theta); % gives (Obs_{x},v,w) vect rotate(ObsJ , ObsI , phi); vect rotate(ObsK , ObsI , phi); % gives (Obs_{x},Obs_{y},Obs_{z})

% The passage matrix P from vect I ,vect J ,vect K % to ObsI ,ObsJ ,ObsK is the matrix % composed of the vectors ObsI ,ObsJ and ObsK expressed % in the base vect I ,vect J ,vect K . % We have X = P X′ where X are the coordinates of a point % in vect I ,vect J ,vect K % and X′ the coordinates of the same point in ObsI ,ObsJ ,ObsK . % In order to get P^{−1}, it suffices to build vectors using % the previous rotations in the inverse order.

vect def vect(IObsI , vect I ); vect def vect(IObsJ , vect J ); vect def vect(IObsK , vect K ); vect rotate(IObsK , IObsI , −phi); vect rotate(IObsJ , IObsI , −phi); vect rotate(IObsK , IObsJ , −theta); vect rotate(IObsI , IObsJ , −theta); vect rotate(IObsJ , IObsK , −psi); vect rotate(IObsI , IObsK , −psi);

enddef; Point of view

% This macro computes the three angles necessary for compute reference % name = name of an instance of an object % target = target point (local to object name) % phi = angle vardef point of view obj (expr name, target, phi) =

define current point offset (name); % enables pnt point of view abs(pnt(target), phi);

enddef;

% Compute absolute perspective. target is an absolute point number vardef point of view abs(expr target, phi) =

save v a, psi, theta; v a := new vect; vect diff (v a, target, Obs); vect mult(v a, v a, 1/vect mod(v a)); psi = angle((vect[v a]x, vect[v a]y)); theta = − angle((vect[v a]x ++ vect[v a]y, vect[v a]z)); compute reference(psi, theta, phi); free vect(v a);

enddef;

% Distance between the observer and point n of object name % Result is put in dist vardef obs distance(text dist)(expr name, n) =

save v a; v a := new vect; define current point offset (name); % enables pnt dist := vect mod(v a, pnt(n), Obs); free vect(v a);

enddef; Vector and point allocation

% Allocation is done through a stack of vectors numeric last vect ; last vect = 0;

% vector allocation def new vect = incr(last vect )

% message & decimal last vect & ;

enddef;

let new point = new vect;

def new points(text p)(expr n) = save p; numeric p[]; for i := 1 upto n: p[i] := new point; endfor;

enddef;

% Free a vector % A vector can only be freed safely when it was the last vector created. def free vect(expr i) =

if i = last vect : last vect := last vect − 1;

else: errmessage( & decimal i & ); fi;

enddef;

let free point = free vect;

def free points(text p)(expr n) =

for i := 10 step −1 until 1: free point(p[i]); endfor; enddef;

Debugging

def show vect(expr t, i) = message & t & & & decimal vect[i]x & & decimal vect[i]y &

& decimal vect[i]z & ;

enddef;

let show point = show vect;

def show pair(expr t, zz) =

message t & & decimal xpart(zz) & & decimal ypart(zz) & ; enddef; Access to object features

% a must be a string representing an class name, such as . % b is the tail of a macro name.

def obj (expr a, b, i) = scantokens begingroup save n; string n; n = a & b & i; n endgroup

enddef;

def obj points (expr name) =

obj (obj class (name), , name) enddef;

def obj faces (expr name) =

obj (obj class (name), , name) enddef;

vardef obj point offset (expr name) =

obj (obj class (name), , name) enddef;

vardef obj face offset (expr name) =

obj (obj class (name), , name) enddef;

def obj class (expr name) = obj (name, , ) enddef; def define point offset (expr name, o) =

begingroup save n, tmpdef ; string n, tmpdef ; n = obj class (name) & & name; expandafter numeric scantokens n; scantokens n := last point offset ; last point offset := last point offset + o; tmpdef = & obj class (name) & & name&

& decimal o & ; scantokens tmpdef ;

endgroup enddef;

def define face offset (expr name, o) =

begingroup save n, tmpdef ; string n, tmpdef ; n = obj class (name) & & name; expandafter numeric scantokens n; scantokens n := last face offset ; last face offset := last face offset + o; tmpdef = & obj class (name) & & name&

& decimal o & ; scantokens tmpdef ;

endgroup enddef;

def define current point offset (expr name) = save current point offset ; numeric current point offset ; current point offset := obj point offset (name);

enddef;

def define current face offset (expr name) = save current face offset ; numeric current face offset ; current face offset := obj face offset (name);

enddef; Drawing an object

% name is an object instance def draw obj (expr name) = project obj (name); draw faces(name); enddef;

Normalization of an object

% This macro translates an object so that a list of vertices is centered % on the origin, and the last vertex is put on a sphere whose radius is 1. % name is the name of the object and vertices is a list % of points whose barycenter will define the center of the object. % (vertices need not be the list of all vertices) vardef normalize obj (expr name)(text vertices) =

save v a, nvertices, last;

numeric v a, last; nvertices = 0; v a = new vect; vect def (v a, 0, 0, 0) forsuffixes $ = vertices:

vect sum(v a, v a, pnt($)); nvertices := nvertices + 1; last := $;

endfor;

vect mult(v a, v a, −1/nvertices);

translate obj (name, v a); % object centered on the origin scale obj (name, 1/vect mod(pnt(last))); free vect(v a);

enddef; General definitions

% Vector arrays numeric vect[]x, vect[]y, vect[]z;

% Observer numeric Obs; Obs = new point; % default value: set point(Obs, 0, 0, 20);

% Observer’s vectors ObsI = new vect; ObsJ = new vect; ObsK = new vect; IObsI = new vect; IObsJ = new vect; IObsK = new vect;

% distance observer/plane (must be > 0) numeric Obs dist; % represents Obs dist × drawing scale % default value: Obs dist = 2; % means Obs dist × drawing scale

% Screen Size % The screen size is defined through two angles: the horizontal field % and the vertical field numeric h field, v field; h field = 100; % degrees v field = 70; % degrees

% Reference vectors 0, i⃗, j⃗ and k⃗ % and their definition numeric vect null, vect I , vect J , vect K ; vect null = new vect; vect I = new vect; vect J = new vect; vect K = new vect; vect def (vect null, 0, 0, 0); vect def (vect I , 1, 0, 0); vect def (vect J , 0, 1, 0); vect def (vect K , 0, 0, 1); numeric point null; point null = vect null;

% Observer’s orientation, defined by three angles numeric Obs psi, Obs theta, Obs phi; % default value: Obs psi = 0; Obs theta = 90; Obs phi = 0;

% Points for the figures numeric points [];

% name is the name of an object instance % npoints is its number of defining points def new obj points(expr name, npoints) =

define point offset (name, npoints); define current point offset (name); for i := 1 upto obj points (name): pnt(i) := new point; endfor;

enddef;

% name is the name of an object instance % nfaces is its number of defining faces def new obj faces(expr name, nfaces) =

define face offset (name, nfaces); define current face offset (name); redefine obj name (name, current face offset + nfaces);

enddef; % Absolute point number corresponding to object point number i % This macro must only be used within the function defining an object % (such as def cube) or the function drawing an object (such as % draw cube). def ipnt (expr i) = i + current point offset enddef; def pnt(expr i) = points [ipnt (i)] enddef;

def face(expr i) = (i + current face offset ) enddef;

% Absolute point number corresponding to local point n % in object instance name vardef pnt obj (expr name, n) =

hide(define current point offset (name); )pnt(n) enddef;

% Absolute face number corresponding to local face n % in object instance name vardef face obj (expr name, n) =

hide(define current face offset (name); )face(n) enddef;

% Scale numeric drawing scale; drawing scale = 2cm;

% Color % This function is useful when a color is expressed in hexadecimal. def hexcolor(expr s) =

(hex(substring(0, 2) of s)/255, hex(substring(2, 4) of s)/255,

hex(substring(4, 6) of s)/255)

enddef;

% Filling and contours boolean filled faces, draw contours; filled faces = true; draw contours = true; numeric contourwidth; % thickness of contours contourwidth = 1pt;

% Overflow control % An overflow can occur when an object is too close from the observer % or if an object is out of sight. We use a special value to mark % coordinates which would lead to an overflow. numeric too big ; too big = 4000;

% Object offset (the points defining an object are arranged % in a single array, and the objects are easier to manipulate % if the point numbers are divided into a number and an offset). numeric last point offset , last face offset ; last point offset = 0; last face offset = 0; Computation of field parameters of an animation

numeric xmin , ymin , xmax , ymax ;

def compute bbox =

if known xmin :

xmin := min(xmin , xpart(llcorner(currentpicture))); ymin := min(ymin , ypart(llcorner(currentpicture))); xmax := max(xmax , xpart(urcorner(currentpicture))); ymax := max(ymax , ypart(urcorner(currentpicture)));

else:

xmin = xpart(llcorner(currentpicture)); ymin = ypart(llcorner(currentpicture)); xmax = xpart(urcorner(currentpicture)); ymax = ypart(urcorner(currentpicture));

fi; enddef;

boolean show animation parameters; show animation parameters = false; numeric paper height; paper height = 29.7; % paper height in cm

% show bounding box of an animation, in PostScript points % and parameters for animation script vardef show animation bbox =

save trx, try, h, w, delta, pnmx, pnmy, pnmw, pnmh, res; res = 36; % 36 dots per inch in bitmap w = xmax − xmin ; h = ymax − ymin ;

if show animation parameters:

message & decimal round(xmin ) & & decimal round(ymin ) & & decimal round(w) & & decimal round(h) & ;

fi; if xmin ≤ 20: trx = 50 − xmin ; else: trx = 0; fi; if ymin ≤ 20: try = 50 − ymin ; else: try = 0; fi;

if show animation parameters:
message

& decimal round(trx) & & decimal round(try); fi; xmin := xmin + trx; ymin := ymin + try; delta = 10; pnmx = round(xmin ∗ (res/72) − delta); pnmy = round((paper height/2.54 ∗ 72 − ymin − h) ∗ (res/72) − delta); pnmw = round(w ∗ (res/72) + 2 ∗ delta); pnmh = round(h ∗ (res/72) + 2 ∗ delta);

if show animation parameters:
message

& decimal pnmx & & decimal pnmy & & decimal pnmw & & decimal pnmh;

fi; write script(round(trx), round(try),

pnmx, pnmy, pnmw, pnmh, res, jobname, );

enddef; Creation of a shell script to automate the animation

% This is UNIX targetted and may need to be customized.

vardef write script(expr trx, try, xmin, ymin, w, h, res, output, file) = save s; string s; def write to file(text arg) = write arg to file; enddef; write to file( ); write to file( );

write to file(  & output &  );
write to file(  & output &  & output &  );

if false: fi% indentation hack for meta-mode.el write to file( ); write to file( ); s := & ditto; s := s & decimal trx & & decimal try & & ditto & ; write to file(s);

% ghostscript PostScript into ppm s := ; s := s & & decimal res & ; write to file(s); write to file( );

% possible alternative:

% s :=  & decimal(w) &  & decimal(h);
% s := s &  & decimal(xmin) &  & decimal(ymin);

% s := s & ;

s :=  & decimal(xmin) &  & decimal(ymin);
s := s &  & decimal(w) &  & decimal(h) &  ;

s := s & ; write to file(s); write to file( ); write to file( ); write to file( & output & ); s := ; s := s & output & & output & ; write to file(s); write to file( & output & ); write to file(EOF ); % end of file

enddef; Standard animation definitions

% These definitions produce one image of some kind.

extra endfig := ;

% In the standard animations, the observer follows a circle, shown below:

% Standard image 1: this is an example and may be adapted. % name is an object instance def one image(expr name, i, a, rd, ang) =

beginfig(i); set point(Obs, −rd ∗ cosd(a ∗ ang), −rd ∗ sind(a ∗ ang), 1); Obs phi := 90; Obs dist := 2; point of view obj (name, 1, Obs phi); % fix point 1 of object name draw obj (name); rotate obj pv(name, 1, vect I , ang); draw point(name, 1); % show the rotation point draw axes(red, green, blue);

endfig; enddef;

% Standard image 2: this is an example and may be adapted. % name a and name b are object instances. def one image two objects(expr name a, name b, i, a, rd, ang) =

beginfig(i); set point(Obs, −rd ∗ cosd(a ∗ ang), −rd ∗ sind(a ∗ ang), 1); Obs phi := 90; Obs dist := 2; point of view obj (name a, 1, Obs phi); % fix point 1 of object name a draw obj (name a); draw obj (name b); rotate obj pv(name a, 1, vect I , ang); rotate obj pv(name b, 13, vect J , −ang);

%rotate obj pp(name b, 13, 7, −ang); draw point(name a, 1); % show the rotation point draw axes(red, green, blue);

endfig; enddef; % Standard image 3: this is an example and may be adapted. % name a, name b and name c are object instances. def one image three objects(expr name a, name b, name c, i, a, rd, ang) =

beginfig(i); set point(Obs, −rd ∗ cosd(a ∗ ang), −rd ∗ sind(a ∗ ang), 1); Obs phi := 90; Obs dist := 2; h field := 100; v field := 150; point of view obj (name a, 1, Obs phi); % fix point 1 of object name a draw obj (name a); draw obj (name b); draw obj (name c); v a := new vect; vect def (v a, .03 ∗ cosd(−a ∗ ang + 90), .03 ∗ sind(−a ∗ ang + 90), 0); translate obj (name c, v a); free vect(v a); rotate obj pv(name a, 1, vect I , ang); rotate obj pv(name b, 13, vect J , −ang);

%rotate obj pp(name b, 13, 7, −ang); draw point(name a, 1); % show the rotation point draw axes(red, green, blue);

endfig; enddef;

% Standard image 4: this is an example and may be adapted. % name a and name b are object instances. def one image two identical objects(expr name a, name b, i, a, rd, ang) =

beginfig(i); set point(Obs, −rd ∗ cosd(a ∗ ang), −rd ∗ sind(a ∗ ang), 2); Obs phi := 90; Obs dist := 2; point of view obj (name a, 1, Obs phi); % fix point 1 of object name a draw obj (name a); draw obj (name b); rotate obj pv(name a, 1, vect I , ang); rotate obj pv(name b, 13, vect J , −ang);

%rotate obj pp(name a, 13, 7, −ang); draw point(name a, 1); % show the rotation point draw axes(red, green, blue);

endfig; enddef; % An animation is a series of images, and these series are produced here.

% Standard animation 1 % name is a class name def animate object(expr name, imin, imax, index) =

numeric ang; ang = 360/(imax − imin + 1); assign obj ( , name); for i := imin upto imax: one image( , i + index, i, 5, ang); endfor; show animation bbox;

enddef;

% Standard animation 2 % name a and name b are class names def animate two objects(expr name a, name b, imin, imax, index) =

numeric ang; ang = 360/(imax − imin + 1);

assign obj (  , name a); assign obj (  , name b);
translate obj (  , vect K ); translate obj (  , vect K );

for i := imin upto imax:

one image two objects( , , i + index, i, 5, ang); endfor; show animation bbox;

enddef;

% Standard animation 3 % name a, name b and name c are class names vardef animate three objects(expr name a, name b, name c, imin, imax, index) =

numeric ang; ang = 360/(imax − imin + 1); assign obj ( , name a); assign obj ( , name b); assign obj ( , name c); scale obj ( , .7); numeric v a; v a := new vect; vect def vect(v a, vect K ); vect mult(v a, v a, 4); put obj ( , v a, 1, 0, 0, 0); free vect(v a); scale obj ( , .5); translate obj ( , vect K ); translate obj ( , vect K ); for i := imin upto imax:

one image three objects( , , , i + index, i, 7, ang); endfor; show animation bbox;

enddef;

% Standard animation 4 % name is a class name def animate two identical objects(expr name, imin, imax, index) =

numeric ang; ang = 360/(imax − imin + 1);

assign obj (  , name); assign obj (  , name);
translate obj (  , vect K ); translate obj (  , vect K );

for i := imin upto imax:

one image two identical objects( , , i + index, i, 10, ang); endfor; show animation bbox;

enddef;

endinput