summaryrefslogtreecommitdiff
path: root/Master/texmf-dist/metapost/mp3d
diff options
context:
space:
mode:
authorKarl Berry <karl@freefriends.org>2006-01-11 22:58:36 +0000
committerKarl Berry <karl@freefriends.org>2006-01-11 22:58:36 +0000
commitac3c55a3216b5988f0e48ba9414ddb059f19a699 (patch)
treea752ab12de05a9ac4511903abc09675172018fd6 /Master/texmf-dist/metapost/mp3d
parentd087712418726a64822e40ce1c0627a514d17975 (diff)
trunk/Master/texmf-dist/metapost
git-svn-id: svn://tug.org/texlive/trunk@104 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Master/texmf-dist/metapost/mp3d')
-rw-r--r--Master/texmf-dist/metapost/mp3d/3d.mp1259
-rw-r--r--Master/texmf-dist/metapost/mp3d/3danim.mp269
-rw-r--r--Master/texmf-dist/metapost/mp3d/3dgeom.mp931
-rw-r--r--Master/texmf-dist/metapost/mp3d/3dpoly.mp294
-rw-r--r--Master/texmf-dist/metapost/mp3d/3dutil.mp37
-rw-r--r--Master/texmf-dist/metapost/mp3d/animpoly.mp35
-rw-r--r--Master/texmf-dist/metapost/mp3d/create_animation.sh15
7 files changed, 2840 insertions, 0 deletions
diff --git a/Master/texmf-dist/metapost/mp3d/3d.mp b/Master/texmf-dist/metapost/mp3d/3d.mp
new file mode 100644
index 00000000000..7335a384d2c
--- /dev/null
+++ b/Master/texmf-dist/metapost/mp3d/3d.mp
@@ -0,0 +1,1259 @@
+%%\input epsf
+%%\def\newpage{\vfill\eject}
+%%\advance\vsize1in
+%%\let\ora\overrightarrow
+%%\def\title#1{\hrule\vskip1mm#1\par\vskip1mm\hrule\vskip5mm}
+%%\def\figure#1{\par\centerline{\epsfbox{#1}}}
+%%\title{{\bf 3D.MP: 3-DIMENSIONAL REPRESENTATIONS IN METAPOST}}
+
+%% version 1.34, 17 August 2003
+%% {\bf Denis Roegel} ({\tt 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.
+
+%% PROJECTS FOR THE FUTURE:
+
+%% $-$ take light sources into account and show shadows and darker faces
+
+%% $-$ handle overlapping of objects ({\it obj\_name\/} can be used when
+%% going through all faces)
+
+if known three_d_version:
+ expandafter endinput % avoids loading this package twice
+fi;
+
+message "*** 3d, v1.34 (c) D. Roegel, 17 August 2003 ***";
+numeric three_d_version;
+three_d_version=1.34;
+
+% This package needs |3dgeom| in a few places. |3dgeom| also loads |3d|
+% but that's not a problem.
+%
+input 3dgeom;
+
+%%\newpage
+%%\title{Vector operations}
+
+% components of vector |i|
+def xval(expr i)=vec[i]x enddef;
+def yval(expr i)=vec[i]y enddef;
+def zval(expr i)=vec[i]z enddef;
+
+% vector (or point) equality (absolute version)
+def vec_eq_(expr i,j)=
+ ((xval(i)=xval(j)) and (yval(i)=yval(j)) and (zval(i)=zval(j)))
+enddef;
+
+% vector (or point) equality (local version)
+def vec_eq(expr i,j)=vec_eq_(pnt(i),pnt(j)) enddef;
+
+% vector inequality (absolute version)
+def vec_neq_(expr i,j)=(not vec_eq_(i,j)) enddef;
+
+% vector inequality (local version)
+def vec_neq(expr i,j)=(not vec_eq(i,j)) enddef;
+
+% definition of vector |i| by its coordinates (absolute version)
+def vec_def_(expr i,xi,yi,zi)= vec[i]x:=xi;vec[i]y:=yi;vec[i]z:=zi; enddef;
+
+% definition of vector |i| by its coordinates (local version)
+def vec_def(expr i,xi,yi,zi)= vec_def_(pnt(i),xi,yi,zi) enddef;
+
+% a point is stored as a vector (absolute version)
+let set_point_ = vec_def_;
+
+% a point is stored as a vector (local version)
+let set_point = vec_def;
+
+def set_point_vec_(expr i,v)=
+ set_point_(i,xval(v),yval(v),zval(v))
+enddef;
+
+def set_point_vec(expr i,v)=set_point_vec_(pnt(i),v) enddef;
+
+let vec_def_vec_=set_point_vec_;
+let vec_def_vec=set_point_vec;
+
+% vector sum: |vec[k]| $\leftarrow$ |vec[i]|$+$|vec[j]| (absolute version)
+def vec_sum_(expr k,i,j)=
+ vec[k]x:=vec[i]x+vec[j]x;
+ vec[k]y:=vec[i]y+vec[j]y;
+ vec[k]z:=vec[i]z+vec[j]z;
+enddef;
+
+% vector sum: |vec[k]| $\leftarrow$ |vec[i]|$+$|vec[j]| (local version)
+def vec_sum(expr k,i,j)=vec_sum_(pnt(k),pnt(i),pnt(j)) enddef;
+
+% vector translation: |vec[i]| $\leftarrow$ |vec[i]|$+$|vec[v]|
+def vec_translate_(expr i,v)=vec_sum_(i,i,v) enddef;
+
+% Here, the second parameter is absolute, because this is probably
+% the most common case.
+def vec_translate(expr i,v)=vec_translate_(pnt(i),v) enddef;
+
+% vector difference: |vec[k]| $\leftarrow$ |vec[i]|$-$|vec[j]|
+def vec_diff_(expr k,i,j)=
+ vec[k]x:=vec[i]x-vec[j]x;
+ vec[k]y:=vec[i]y-vec[j]y;
+ vec[k]z:=vec[i]z-vec[j]z;
+enddef;
+
+def vec_diff(expr k,i,j)=vec_diff_(pnt(k),pnt(i),pnt(j)) enddef;
+
+% dot product of |vec[i]| and |vec[j]|
+vardef vec_dprod_(expr i,j)=
+ (vec[i]x*vec[j]x+vec[i]y*vec[j]y+vec[i]z*vec[j]z)
+enddef;
+
+vardef vec_dprod(expr i,j)=vec_dprod_(pnt(i),pnt(j)) enddef;
+
+% modulus of |vec[i]|, absolute version
+% In the computation, we try to avoid overflows or underflows;
+% we perform a scaling in order to avoid losing too much
+% information in certain cases
+vardef vec_mod_(expr i)=
+ save prod,m_;
+ hide(
+ new_vec(v_a);
+ m_=max(abs(xval(i)),abs(yval(i)),abs(zval(i)));
+ if m_>0:vec_mult_(v_a,i,1/m_);else:vec_def_vec_(v_a,vec_null);fi;
+ prod=m_*sqrt(vec_dprod_(v_a,v_a));
+ free_vec(v_a);
+ )
+ prod
+enddef;
+
+% modulus of |vec[i]|, local version
+% If the return value must be compared to 0,
+% use |vec_eq| with |vec_null| instead.
+vardef vec_mod(expr i)= vec_mod_(pnt(i)) enddef;
+
+% unit vector |vec[i]| corresponding to vector |vec[j]|
+% only non-null vectors are changed
+def vec_unit_(expr i,j)=
+ if vec_mod_(j)>0: vec_mult_(i,j,1/vec_mod_(j));
+ else:vec_def_vec_(i,j);
+ fi;
+enddef;
+
+def vec_unit(expr i,j)=vec_unit_(pnt(i),pnt(j)) enddef;
+
+% vector product: |vec[k]| $\leftarrow$ |vec[i]| $\land$ |vec[j]|
+def vec_prod_(expr k,i,j)=
+ vec[k]x:=vec[i]y*vec[j]z-vec[i]z*vec[j]y;
+ vec[k]y:=vec[i]z*vec[j]x-vec[i]x*vec[j]z;
+ vec[k]z:=vec[i]x*vec[j]y-vec[i]y*vec[j]x;
+enddef;
+
+def vec_prod(expr k,i,j)=vec_prod_(pnt(k),pnt(i),pnt(j)) enddef;
+
+% scalar multiplication: |vec[j]| $\leftarrow$ |vec[i]*v| (absolute version)
+def vec_mult_(expr j,i,v)=
+ vec[j]x:=v*vec[i]x;vec[j]y:=v*vec[i]y;vec[j]z:=v*vec[i]z;
+enddef;
+
+% scalar multiplication: |vec[j]| $\leftarrow$ |vec[i]*v| (local version)
+def vec_mult(expr j,i,v)=vec_mult_(pnt(j),pnt(i),v) enddef;
+
+% middle of two points (absolute version)
+def mid_point_(expr k,i,j)= vec_sum_(k,i,j);vec_mult_(k,k,.5); enddef;
+
+% middle of two points (local version)
+def mid_point(expr k,i,j)= mid_point_(pnt(k),pnt(i),pnt(j)); enddef;
+
+%%\newpage
+%%\title{Vector rotation}
+% Rotation of |vec[v]| around |vec[axis]| by an angle |alpha|
+
+%% The vector $\vec{v}$ is first projected on the axis
+%% giving vectors $\vec{a}$ and $\vec{h}$:
+%%\figure{vect-fig.9}
+%% If we set
+%% $\vec{b}={\ora{axis}\over \left\Vert\vcenter{\ora{axis}}\right\Vert}$,
+%% the rotated vector $\vec{v'}$ is equal to $\vec{h}+\vec{f}$
+%% where $\vec{f}=\cos\alpha \cdot \vec{a} + \sin\alpha\cdot \vec{c}$.
+%% and $\vec{h}=(\vec{v}\cdot\vec{b})\vec{b}$
+%%\figure{vect-fig.10}
+
+% The rotation is independent of |vec[axis]|'s module.
+% |v| = old and new vector
+% |axis| = rotation axis
+% |alpha| = rotation angle
+%
+vardef vec_rotate_(expr v,axis,alpha)=
+ new_vec(v_a);new_vec(v_b);new_vec(v_c);
+ new_vec(v_d);new_vec(v_e);new_vec(v_f);
+ new_vec(v_g);new_vec(v_h);
+ vec_mult_(v_b,axis,1/vec_mod_(axis));
+ vec_mult_(v_h,v_b,vec_dprod_(v_b,v)); % projection of |v| on |axis|
+ vec_diff_(v_a,v,v_h);
+ vec_prod_(v_c,v_b,v_a);
+ vec_mult_(v_d,v_a,cosd(alpha));
+ vec_mult_(v_e,v_c,sind(alpha));
+ vec_sum_(v_f,v_d,v_e);
+ vec_sum_(v,v_f,v_h);
+ free_vec(v_h);free_vec(v_g);
+ free_vec(v_f);free_vec(v_e);free_vec(v_d);
+ free_vec(v_c);free_vec(v_b);free_vec(v_a);
+enddef;
+
+% The second parameter is left absolute because this is probably the most
+% common case.
+vardef vec_rotate(expr v,axis,alpha)=vec_rotate_(pnt(v),axis,alpha) enddef;
+
+%%\newpage
+%%\title{Operations on objects}
+% |iname| is the handler for an instance of an object of class |name|
+% |iname| must be a letter string
+% |vardef| is not used because at some point we give other names
+% to |assign_obj| with |let| and this cannot be done with |vardef|.
+% (see MFbook for details)
+def assign_obj(expr iname,name)=
+ begingroup
+ save tmpdef;
+ string tmpdef; % we need to add double quotes (char 34)
+ tmpdef="def " & iname & "_class=" & ditto & name & ditto & " enddef";
+ scantokens tmpdef;
+ def_obj(iname);
+ endgroup
+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="def_" & obj_class_(name) & "(" & ditto & name & ditto & ")";
+ tmpdef
+ endgroup
+enddef;
+
+% This macro puts an object back where it was right at the beginning,
+% or rather, where the |set| definition puts it (which may be different
+% than the initial position, in case it depends on parameters).
+% |iname| is the name of an object instance.
+vardef reset_obj(expr iname)=
+ save tmpdef;
+ string tmpdef;
+ define_current_point_offset_(iname);
+ tmpdef="set_" & obj_class_(iname) & "_points";
+ 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.
+% If the shape of the object has been changed since it was
+% created, these changes are lost.
+vardef put_obj(expr iname,pos,scale,psi,theta,phi)=
+ reset_obj(iname);scale_obj(iname,scale);
+ new_vec(v_x);new_vec(v_y);new_vec(v_z);
+ vec_def_vec_(v_x,vec_I);
+ vec_def_vec_(v_y,vec_J);
+ vec_def_vec_(v_z,vec_K);
+ rotate_obj_abs_pv(iname,point_null,v_z,psi);
+ vec_rotate_(v_x,v_z,psi);vec_rotate_(v_y,v_z,psi);
+ rotate_obj_abs_pv(iname,point_null,v_y,theta);
+ vec_rotate_(v_x,v_y,theta);vec_rotate_(v_z,v_y,theta);
+ rotate_obj_abs_pv(iname,point_null,v_x,phi);
+ vec_rotate_(v_y,v_x,phi);vec_rotate_(v_z,v_x,phi);
+ free_vec(v_z);free_vec(v_y);free_vec(v_x);
+ translate_obj(iname,pos);
+enddef;
+
+%%\newpage
+%%\title{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 |vec[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)=
+ define_current_point_offset_(name);
+ new_vec(v_a);
+ for i:=1 upto obj_points_(name):
+ vec_diff_(v_a,pnt(i),p);
+ vec_rotate_(v_a,v,a);
+ vec_sum_(pnt(i),v_a,p);
+ endfor;
+ free_vec(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 $\ora{pq}$. The angle of rotation is |a|.
+vardef rotate_obj_pp(expr name,p,q,a)=
+ define_current_point_offset_(name);
+ new_vec(v_a);new_vec(axis);
+ vec_diff_(axis,pnt(q),pnt(p));
+ for i:=1 upto obj_points_(name):
+ vec_diff_(v_a,pnt(i),pnt(p));
+ vec_rotate_(v_a,axis,a);
+ vec_sum_(pnt(i),v_a,pnt(p));
+ endfor;
+ free_vec(axis);free_vec(v_a);
+enddef;
+
+% Translation of an object instance |name| by a vector |vec[v]|.
+vardef translate_obj(expr name,v)=
+ define_current_point_offset_(name);
+ for i:=1 upto obj_points_(name):
+ vec_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):
+ vec_mult(i,i,v);
+ endfor;
+enddef;
+
+
+%%\newpage
+%%\title{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 $\alpha$ equal to the angle from |i| to |j|.
+%%\figure{vect-fig.11}
+vardef rotate_in_plane_(expr k,o,i,j)=
+ save cosalpha,sinalpha,alpha;
+ new_vec(v_a);new_vec(v_b);new_vec(v_c);
+ vec_diff_(v_a,i,o);vec_diff_(v_b,j,o);vec_prod_(v_c,v_a,v_b);
+ cosalpha=vec_dprod_(v_a,v_b)/vec_mod_(v_a)/vec_mod_(v_b);
+ sinalpha=sqrt(1-cosalpha**2);
+ alpha=angle((cosalpha,sinalpha));
+ vec_rotate_(v_b,v_c,alpha);
+ vec_sum_(k,o,v_b);
+ free_vec(v_c);free_vec(v_b);free_vec(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
+%% $\widehat{(\ora{om},\ora{mc})}=\alpha$
+%% This is useful to define regular polyhedra
+%%\figure{vect-fig.7}
+vardef new_face_point_(expr c,o,i,j,alpha)=
+ new_vec(v_a);new_vec(v_b);new_vec(v_c);new_vec(v_d);new_vec(v_e);
+ vec_diff_(v_a,i,o);vec_diff_(v_b,j,o);
+ vec_sum_(v_c,v_a,v_b);
+ vec_mult_(v_d,v_c,.5);
+ vec_diff_(v_e,i,j);
+ vec_sum_(c,v_d,o);
+ vec_rotate_(v_d,v_e,alpha);
+ vec_sum_(c,v_d,c);
+ free_vec(v_e);free_vec(v_d);free_vec(v_c);free_vec(v_b);free_vec(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;
+
+%%\newpage
+%%\title{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
+%%\figure{vect-fig.8}
+vardef project_point(expr p,m)=
+ save tmpalpha;
+ new_vec(v_a);new_vec(v_b);
+ if projection_type=2: % oblique
+ if point_in_plane_p_pl_(m)(projection_plane):
+ % |m| is on the projection plane
+ vec_diff_(v_a,m,ObliqueCenter_);
+ y[p]:=drawing_scale*vec_dprod_(v_a,ProjJ_);
+ x[p]:=drawing_scale*vec_dprod_(v_a,ProjK_);
+ else: % |m| is not on the projection plane
+ new_line_(l)(m,ObliqueCenter_);
+ vec_diff_(l2,l2,Obs);
+ vec_sum_(l2,l2,m);
+ % (the direction does not depend on Obs)
+ if def_inter_p_l_pl_(v_a)(l)(projection_plane):
+ vec_diff_(v_a,v_a,ObliqueCenter_);
+ y[p]:=drawing_scale*vec_dprod_(v_a,ProjJ_);
+ x[p]:=drawing_scale*vec_dprod_(v_a,ProjK_);
+ else: message "Point " & decimal m & " cannot be projected";
+ x[p]:=too_big_;y[p]=too_big_;
+ fi;
+ free_line(l);
+ fi;
+ else:
+ vec_diff_(v_b,m,Obs); % vector |Obs|-|m|
+ % |vec[v_a]| is |vec[v_b]| expressed in (|ObsI_|,|ObsJ_|,|ObsK_|)
+ % coordinates.
+ vec[v_a]x:=vec[IObsI_]x*vec[v_b]x
+ +vec[IObsJ_]x*vec[v_b]y+vec[IObsK_]x*vec[v_b]z;
+ vec[v_a]y:=vec[IObsI_]y*vec[v_b]x
+ +vec[IObsJ_]y*vec[v_b]y+vec[IObsK_]y*vec[v_b]z;
+ vec[v_a]z:=vec[IObsI_]z*vec[v_b]x
+ +vec[IObsJ_]z*vec[v_b]y+vec[IObsK_]z*vec[v_b]z;
+ if vec[v_a]x<Obs_dist: % then, point |m| is too close
+ message "Point " & decimal m & " too close -> not drawn";
+ x[p]:=too_big_;y[p]=too_big_;
+ else:
+ if (angle(vec[v_a]x,vec[v_a]z)>h_field/2)
+ or (angle(vec[v_a]x,vec[v_a]y)>v_field/2):
+ message "Point " & decimal m & " out of screen -> not drawn";
+ x[p]:=too_big_;y[p]=too_big_;
+ else:
+ if projection_type=0: % central perspective
+ tmpalpha:=Obs_dist/vec[v_a]x;
+ else:
+ tmpalpha:=1; % parallel
+ fi;
+ y[p]:=drawing_scale*tmpalpha*vec[v_a]y;
+ x[p]:=drawing_scale*tmpalpha*vec[v_a]z;
+ fi;
+ fi;
+ fi;
+ free_vec(v_b);free_vec(v_a);
+enddef;
+
+% At some point, we may need to do an oblique projection
+% of vectors |ObsK_| and |ObsI_| on a plane, and to normalize
+% and orthogonalize the projections (with the projection of |ObsK_|
+% keeping the same direction). This is done here,
+% where we take two vectors, a direction (line) and
+% a plane, and return two vectors. This function assumes
+% there is an intersection between line |l| and plane |p|.
+% We do not test it here.
+
+vardef project_vectors(expr va,vb)(expr k,i)(text l)(text p)=
+ save vc;new_vec(vc);
+ if proj_v_v_l_pl_(va,k)(l)(p): % |va| is the projection of vector |k|
+ else: message "THIS SHOULD NOT HAPPEN";
+ fi;
+ if proj_v_v_l_pl_(vb,i)(l)(p): % |vb| is the projection of vector |i|
+ else: message "THIS SHOULD NOT HAPPEN";
+ fi;
+ % now, we orthonormalize these vectors:
+ vec_prod_(vc,va,vb);
+ vec_unit_(va,va);vec_unit_(vc,vc);vec_prod_(vb,vc,va);
+ free_vec(vc);
+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;
+
+% Projection screen
+vardef show_projection_screen=
+ save dx,dy;
+ dx=Obs_dist*sind(h_field/2)/cosd(h_field/2);
+ dy=Obs_dist*sind(v_field/2)/cosd(v_field/2);
+ new_vec(pa);new_vec(pb);new_vec(pc);new_vec(pd);new_vec(op);
+ new_vec(w);new_vec(h);
+ vec_mult_(op,ObsI_,Obs_dist);vec_sum_(op,op,Obs); % center of screen
+ vec_mult_(w,ObsK_,dx);vec_mult_(h,ObsJ_,dy);
+ vec_sum_(pa,op,w);vec_sum_(pa,pa,h); % upper right corner
+ vec_mult_(w,w,-2);vec_mult_(h,h,-2);
+ vec_sum_(pb,pa,w);vec_sum_(pc,pb,h);vec_sum_(pd,pa,h);
+ message "Screen at corners:";
+ show_point("urcorner: ",pa);
+ show_point("ulcorner: ",pb);
+ show_point("llcorner: ",pc);
+ show_point("lrcorner: ",pd);
+ show_point("Obs:",Obs);
+ free_vec(h);free_vec(w);
+ free_vec(op);free_vec(pd);free_vec(pc);free_vec(pb);free_vec(pa);
+enddef;
+
+
+%%\newpage
+%%\title{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,overflow,i,j,k,nv;
+ path p;boolean overflow;
+ overflow=false;
+ forsuffixes $=vertices:
+ if z[ipnt_($)]=(too_big_,too_big_):overflow:=true; fi;
+ exitif overflow;
+ endfor;
+ if overflow: message "Face can not be drawn, due to overflow";
+ else:
+ p=forsuffixes $=vertices:z[ipnt_($)]--endfor cycle;
+ % we do now search for three distinct and non-aligned suffixes:
+ % usually, the first three suffixes do
+ new_vec(normal_vec);new_vec(v_a);new_vec(v_b);new_vec(v_c);
+ % first, we copy all the indexes in an array, so that
+ % it is easier to go through them
+ i=1; % num0 is not used
+ forsuffixes $=vertices:num[i]=$;i:=i+1;endfor;
+ nv=i-1;
+ for $:=1 upto nv:
+ for $$:=$+1 upto nv:
+ for $$$:=$$+1 upto nv:
+ vec_diff_(v_a,pnt(num[$$]),pnt(num[$]));
+ vec_diff_(v_b,pnt(num[$$$]),pnt(num[$$]));
+ vec_prod_(normal_vec,v_a,v_b);
+ exitif vec_neq_(normal_vec,vec_null);
+ % |vec_mod_| must not be used for such a test
+ endfor;
+ exitif vec_neq_(normal_vec,vec_null);
+ endfor;
+ exitif vec_neq_(normal_vec,vec_null);
+ endfor;
+ if projection_type=0: % perspective
+ vec_diff_(v_c,pnt(num1),Obs);
+ else: % parallel
+ vec_def_vec_(v_c,ObsI_);
+ fi;
+ if filled_faces:
+ if vec_dprod_(normal_vec,v_c)<0:
+ fill p withcolor col;drawcontour(p,contour_width,contour_color)();
+ else: % |draw p dashed evenly;| if this is done, you must ensure
+ % that hidden faces are (re)drawn at the end
+ fi;
+ else:
+ if vec_dprod_(normal_vec,v_c)<0:%visible
+ if not drawhidden:drawcontour(p,contour_width,contour_color)();fi;
+ else: % hidden
+ if drawhidden:
+ drawcontour(p,contour_width,contour_color)(dashed evenly);
+ fi;
+ fi;
+ fi;
+ free_vec(v_c);free_vec(v_b);free_vec(v_a);free_vec(normal_vec);
+ fi;
+enddef;
+
+% |p| is the path to draw (a face contour), |thickness| is the pen width
+% |col| is the color and |type| is a line modifier.
+def drawcontour(expr p,thickness,col)(text type)=
+ if draw_contours and (thickness>0):
+ pickup pencircle scaled thickness;
+ draw p withcolor background; % avoid strange overlapping dashes
+ draw p type withcolor col;
+ pickup pencircle scaled .4pt;
+ fi;
+enddef;
+
+%%\newpage
+% Variables for face handling. First, we have an array for lists of vertices
+% corresponding to faces.
+string face_points_[];% analogous to |vec| 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:="draw_face(" & face_points_[face(i)]
+ & ")(hexcolor(" & ditto & face_color_[face(i)] & ditto
+ & "),true)";scantokens tmpdef;
+ endfor;
+ % then, the visible faces:
+ for i:=1 upto obj_faces_(name):
+ tmpdef:="draw_face(" & face_points_[face(i)]
+ & ")(hexcolor(" & ditto & face_color_[face(i)] & ditto
+ & "),false)";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,vec_null);
+ project_point(2,vec_I);
+ project_point(3,vec_J);
+ project_point(4,vec_K);
+ if (z1<>(too_big_,too_big_)):
+ if (z2<>(too_big_,too_big_)):
+ drawarrow z1--z2 dashed evenly withcolor r;
+ fi;
+ if (z3<>(too_big_,too_big_)):
+ drawarrow z1--z3 dashed evenly withcolor g;
+ fi;
+ if (z4<>(too_big_,too_big_)):
+ drawarrow z1--z4 dashed evenly withcolor b;
+ fi;
+ fi;
+enddef;
+
+% Draw a polygonal line through the list of points
+% This implementation does not work if you call
+% |draw_lines(i,i+4)| because \MP{} adds parentheses around
+% the value of |i|.
+def draw_lines(text vertices)=
+ begingroup % so that we can |let| |draw_lines|
+ save j,num,np;
+ % first, we copy all the indexes in an array, so that
+ % it is easier to go through them
+ j=1;
+ for $=vertices:num[j]=$;j:=j+1;endfor;
+ np=j-1;
+ for j:=1 upto np-1:
+ draw z[ipnt_(num[j])]--z[ipnt_(num[j+1])];
+ endfor;
+ endgroup
+enddef;
+
+let draw_line=draw_lines;
+
+% Draw an arrow between points |i| and |j| of current object
+% This is used from the |draw| definition of an object.
+def draw_arrow(expr i,j)=
+ drawarrow z[ipnt_(i)]--z[ipnt_(j)];
+enddef;
+
+% Draw a line between points |i| of object |obja| and |j| of |objb|
+% This is used when outside an object (i.e., we can't presuppose
+% any object offset)
+vardef draw_line_inter(expr obja, i, objb, j)=
+ project_point(1,pnt_obj(obja,i));
+ project_point(2,pnt_obj(objb,j));
+ draw z1--z2;
+enddef;
+
+% Draw an arrow between points |i| of object |obja| and |j| of |objb|
+% This is used when outside an object (i.e., we can't presuppose
+% any object offset)
+vardef draw_arrow_inter(expr obja, i, objb, j)=
+ project_point(1,pnt_obj(obja,i));
+ project_point(2,pnt_obj(objb,j));
+ draw z1--z2;
+enddef;
+
+%%\newpage
+% 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_="def obj_name(expr i)=if i<1:";
+
+% |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_ & "elseif i<=" & decimal n & ":" & ditto
+ & name & ditto;
+ scantokens begingroup index_to_name_ & "fi;enddef;" 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;
+
+% |i| is a local face number of object |inst|
+% |rgbcolor| is a string representing a color in rgb hexadecimal
+def set_obj_face_color(expr inst,i,rgbcolor)=
+ face_color_[face_obj(inst,i)]:=rgbcolor;
+enddef;
+
+
+%%\newpage
+%%\title{Compute the vectors corresponding to the observer's viewpoint}
+% (vectors |ObsI_|,|ObsJ_| and |ObsK_| in the |vec_I|,|vec_J|,
+% |vec_K| reference; and vectors |IObsI_|,|IObsJ_| and |IObsK_|
+% which are |vec_I|,|vec_J|,|vec_K|
+% in the |ObsI_|,|ObsJ_|,|ObsK_| reference)
+%%\figure{vect-fig.16}
+%% (here, $\psi>0$, $\theta<0$ and $\phi>0$; moreover,
+%% $\vert\theta\vert \leq 90^\circ$)
+
+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_| $\land$ |ObsJ_|)
+ % The vectors are found by rotations of |vec_I|,|vec_J|,|vec_K|.
+ vec_def_vec_(ObsI_,vec_I);vec_def_vec_(ObsJ_,vec_J);
+ vec_def_vec_(ObsK_,vec_K);
+ vec_rotate_(ObsI_,ObsK_,psi);
+ vec_rotate_(ObsJ_,ObsK_,psi);% gives ($u$,$v$,$z$)
+ vec_rotate_(ObsI_,ObsJ_,theta);
+ vec_rotate_(ObsK_,ObsJ_,theta);% gives ($Obs_x$,$v$,$w$)
+ vec_rotate_(ObsJ_,ObsI_,phi);
+ vec_rotate_(ObsK_,ObsI_,phi);% gives ($Obs_x$,$Obs_y$,$Obs_z$)
+ % The passage matrix $P$ from |vec_I|,|vec_J|,|vec_K|
+ % to |ObsI_|,|ObsJ_|,|ObsK_| is the matrix
+ % composed of the vectors |ObsI_|,|ObsJ_| and |ObsK_| expressed
+ % in the base |vec_I|,|vec_J|,|vec_K|.
+ % We have $X=P X'$ where $X$ are the coordinates of a point
+ % in |vec_I|,|vec_J|,|vec_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.
+ vec_def_vec_(IObsI_,vec_I);vec_def_vec_(IObsJ_,vec_J);
+ vec_def_vec_(IObsK_,vec_K);
+ vec_rotate_(IObsK_,IObsI_,-phi);vec_rotate_(IObsJ_,IObsI_,-phi);
+ vec_rotate_(IObsK_,IObsJ_,-theta);vec_rotate_(IObsI_,IObsJ_,-theta);
+ vec_rotate_(IObsJ_,IObsK_,-psi);vec_rotate_(IObsI_,IObsK_,-psi);
+enddef;
+
+%%\newpage
+%%\title{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
+% |phi| = angle
+% This function also computes two vectors needed in case
+% of an oblique projection.
+vardef point_of_view_abs(expr target,phi)=
+ save psi,theta;
+ new_vec(v_a);
+ vec_diff_(v_a,target,Obs);
+ vec_mult_(v_a,v_a,1/vec_mod_(v_a));
+ psi=angle((vec[v_a]x,vec[v_a]y));
+ theta=-angle((vec[v_a]x++vec[v_a]y,vec[v_a]z));
+ compute_reference(psi,theta,phi);
+ if projection_type=2: % oblique
+ % we start by checking that at a minimum the three points defining
+ % the projection plane have different indexes; it doesn't mean
+ % the plane if well defined, but if two values are identical,
+ % the plane can't be well defined.
+ if ((projection_plane1<>projection_plane2) and
+ (projection_plane1<>projection_plane3) and
+ (projection_plane2<>projection_plane3)):
+ new_line_(l)(Obs,Obs);
+ vec_sum_(l2,ObsI_,Obs);
+ if def_inter_p_l_pl_(ObliqueCenter_)(l)(projection_plane):
+ project_vectors(ProjK_,ProjJ_)(ObsK_,ObsJ_)(l)(projection_plane);
+ % define the projection direction
+ set_line_(projection_direction)(Obs,ObliqueCenter_);
+ else:
+ message "Anomalous oblique projection:";
+ message " the observer is watching parallely to the plane";
+ fi;
+ free_line(l);
+ else:
+ message "Anomalous projection plane; did you define it?";
+ fi;
+ fi;
+ free_vec(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)=
+ new_vec(v_a);
+ define_current_point_offset_(name);% enables |pnt|
+ dist:=vec_mod_(v_a,pnt(n),Obs);
+ free_vec(v_a);
+enddef;
+
+%%\newpage
+%%\title{Vector and point allocation}
+% Allocation is done through a stack of vectors
+numeric last_vec_;
+last_vec_=0;
+
+% vector allocation
+% (this must not be a |vardef| because the vector |v| saved is not saved
+% in this macro, but in the calling context)
+def new_vec(text v)=
+ save v;
+ new_vec_(v);
+enddef;
+
+def new_vec_(text v)=
+ v:=incr(last_vec_);
+ %|message "Vector " & decimal (last_vec_+1) & " allocated";|
+enddef;
+
+let new_point = new_vec;
+let new_point_ = new_vec_;
+
+def new_points(text p)(expr n)=
+ save p;
+ numeric p[];
+ for i:=1 upto n:new_point_(p[i]);endfor;
+enddef;
+
+% Free a vector
+% A vector can only be freed safely when it was the last vector created.
+def free_vec(expr i)=
+ if i=last_vec_: last_vec_:=last_vec_-1;
+ %|message "Vector " & decimal i & " freed";|
+ else: errmessage("Vector " & decimal i & " can't be freed!");
+ fi;
+enddef;
+
+let free_point = free_vec;
+
+def free_points(text p)(expr n)=
+ for i:=n step-1 until 1:free_point(p[i]);endfor;
+enddef;
+
+%%\title{Debugging}
+
+def show_vec(expr t,i)=
+ message "Vector " & t & "="
+ & "(" & decimal vec[i]x & "," & decimal vec[i]y & ","
+ & decimal vec[i]z & ")";
+enddef;
+
+% One can write |show_point("2",pnt_obj("obj",2));|
+let show_point=show_vec;
+
+def show_pair(expr t,zz)=
+ message t & "=(" & decimal xpart(zz) & "," & decimal ypart(zz) & ")";
+enddef;
+
+%%\newpage
+%%\title{Access to object features}
+% |a| must be a string representing a class name, such as |"dodecahedron"|.
+% |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),"_points",name)
+enddef;
+
+def obj_faces_(expr name)=
+ obj_(obj_class_(name),"_faces",name)
+enddef;
+
+vardef obj_point_offset_(expr name)=
+ obj_(obj_class_(name),"_point_offset",name)
+enddef;
+
+vardef obj_face_offset_(expr name)=
+ obj_(obj_class_(name),"_face_offset",name)
+enddef;
+
+def obj_class_(expr name)=obj_(name,"_class","") enddef;
+
+%%\newpage
+def define_point_offset_(expr name,o)=
+ begingroup save n,tmpdef;
+ string n,tmpdef;
+ n=obj_class_(name) & "_point_offset" & name;
+ expandafter numeric scantokens n;
+ scantokens n:=last_point_offset_;
+ last_point_offset_:=last_point_offset_+o;
+ tmpdef="def " & obj_class_(name) & "_points" & name &
+ "=" & decimal o & " enddef";
+ scantokens tmpdef;
+ endgroup
+enddef;
+
+def define_face_offset_(expr name,o)=
+ begingroup save n,tmpdef;
+ string n,tmpdef;
+ n=obj_class_(name) & "_face_offset" & name;
+ expandafter numeric scantokens n;
+ scantokens n:=last_face_offset_;
+ last_face_offset_:=last_face_offset_+o;
+ tmpdef="def " & obj_class_(name) & "_faces" & name &
+ "=" & decimal o & " enddef";
+ 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;
+
+
+%%\newpage
+%%\title{Drawing an object}
+% |name| is an object instance
+vardef draw_obj(expr name)=
+ save tmpdef;
+ string tmpdef;
+ current_obj:=name;
+ tmpdef="draw_" & obj_class_(name);
+ project_obj(name);% compute screen coordinates
+ save overflow; boolean overflow; overflow=false;
+ for $:=1 upto obj_points_(name):
+ if z[ipnt_($)]=(too_big_,too_big_):overflow:=true;
+ x[ipnt_($)] := 10; % so that the figure can be drawn anyway
+ y[ipnt_($)] := 10;
+ % why can't I write z[ipnt_($)]:=(10,10); ?
+ fi;
+ exitif overflow;
+ endfor;
+ if overflow:
+ message "Figure has overflows";
+ message " (at least one point is not visible ";
+ message " and had to be drawn at a wrong place)";
+ fi;
+ scantokens tmpdef(name);
+enddef;
+
+%%\title{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 nvertices,last;
+ nvertices=0;
+ new_vec(v_a);vec_def_(v_a,0,0,0)
+ forsuffixes $=vertices:
+ vec_sum_(v_a,v_a,pnt($));
+ nvertices:=nvertices+1;
+ last:=$;
+ endfor;
+ vec_mult_(v_a,v_a,-1/nvertices);
+ translate_obj(name,v_a);% object centered on the origin
+ scale_obj(name,1/vec_mod(last));
+ free_vec(v_a);
+enddef;
+
+
+%%\newpage
+%%\title{General definitions}
+% Vector arrays
+numeric vec[]x,vec[]y,vec[]z;
+
+% Reference vectors $\vec{0}$, $\vec{\imath}$, $\vec{\jmath}$ and $\vec{k}$
+% and their definition
+new_vec(vec_null);new_vec(vec_I);new_vec(vec_J);new_vec(vec_K);
+vec_def_(vec_null,0,0,0);
+vec_def_(vec_I,1,0,0);vec_def_(vec_J,0,1,0);vec_def_(vec_K,0,0,1);
+numeric point_null;
+point_null=vec_null;
+
+% Observer
+new_point(Obs);
+% default value:
+set_point_(Obs,0,0,20);
+
+% Observer's vectors
+new_vec(ObsI_);new_vec(ObsJ_);new_vec(ObsK_);
+% default values:
+vec_def_vec_(ObsI_,vec_I);
+vec_def_vec_(ObsJ_,vec_J);
+vec_def_vec_(ObsK_,vec_K);
+
+new_vec(IObsI_);new_vec(IObsJ_);new_vec(IObsK_);
+
+% These vectors will be vectors of the projection plane,
+% in case of oblique projections:
+new_vec(ProjK_);new_vec(ProjJ_); % there is no |ProjI_|
+
+% This will be the center of the projection plane, in oblique projections
+new_point(ObliqueCenter_);
+
+
+% distance observer/plane (must be $>0$)
+numeric Obs_dist; % represents |Obs_dist| $\times$ |drawing_scale|
+% default value:
+Obs_dist=2; % means |Obs_dist| $\times$ |drawing_scale|
+
+% current object being drawn
+string current_obj;
+
+% kind of projection: 0 for linear (or central) perspective, 1 for parallel,
+% 2 for oblique projection
+% (default is 0)
+numeric projection_type;
+projection_type:=0;
+
+% Definition of a projection plane (only used in oblique projections)
+%
+new_plane_(projection_plane)(1,1,1); % the initial value is irrelevant
+
+% Definition of a projection direction (only used in oblique projections)
+new_line_(projection_direction)(1,1); % the initial value is irrelevant
+
+% this positions the observer at vector |p| (the point observed)
+% + |d| (distance) * (k-(i+j))
+def isometric_projection(expr i,j,k,p,d,phi)=
+ trimetric_projection(i,j,k,1,1,1,p,d,phi);
+enddef;
+
+% this positions the observer at vector |p| (the point observed)
+% + |d| (distance) * (ak-(i+j))
+def dimetric_projection(expr i,j,k,a,p,d,phi)=
+ trimetric_projection(i,j,k,1,1,a,p,d,phi);
+enddef;
+
+% this positions the observer at vector |p| (the point observed)
+% + |d| (distance) * (k-(i+j))
+% |a|, |b| and |c| are multiplicative factors to vectors |i|, |j| and |k|
+vardef trimetric_projection(expr i,j,k,a,b,c,p,d,phi)=
+ save v_a,v_b,v_c;
+ new_vec(v_a);new_vec(v_b);new_vec(v_c);
+ vec_mult_(v_a,i,a);vec_mult_(v_b,j,b);vec_mult_(v_c,k,c);
+ vec_sum_(Obs,v_a,v_b);
+ vec_diff_(Obs,v_c,Obs);
+ vec_mult_(Obs,Obs,d);
+ vec_sum_(Obs,Obs,p);
+ point_of_view_abs(p,phi);
+ projection_type:=1;
+ free_vec(v_c);free_vec(v_b);free_vec(v_a);
+enddef;
+
+% |hor| is an horizontal plane (in the sense that it will represent
+% the horizontal for the observer)
+% |p| is the point in space that the observer targets (center of screen)
+% |a| is an angle (45 degrees corresponds to cavalier drawing)
+% |b| is an angle (see examples defined below)
+% |d| is the distance of the observer
+vardef oblique_projection(text hor)(expr p,a,b,d)=
+ save _l,v_a,v_b,v_c,xxx_,obsJangle_;
+ new_vec(v_a);new_vec(v_b);new_vec(v_c);
+ % we first compute a horizontal line:
+ new_line_(_l)(1,1);
+ if def_inter_l_pl_pl(_l)(hor)(projection_plane):
+ vec_diff_(v_a,_l2,_l1); % horizontal vector
+ % then, we find a normal to the projection plane:
+ def_normal_p_(v_b)(projection_plane);
+ % complete the line and the vector by a third vector (=vertical)
+ vec_prod_(v_c,v_a,v_b);
+ % we make |v_a| a copy of |v_b| since we no longer need |v_b|
+ vec_def_vec_(v_a,v_b);
+ % we rotate |v_b| by an angle |a| around |v_c|
+ vec_rotate_(v_b,v_c,a);
+ % we rotate |v_b| by an angle |b| around |v_a|
+ vec_rotate_(v_b,v_a,b);
+ % we put the observer at the distance |d| of |p| in
+ % the direction of |v_b|:
+ vec_unit_(v_b,v_b);
+ vec_mult_(v_b,v_b,d);vec_sum_(Obs,p,v_b);
+ % We now have to make sure that point |p| and point |Obs|
+ % are on different sides of the projection plane. For this,
+ % we compute two dot products:
+ new_vec(v_d);new_vec(v_e);
+ vec_diff_(v_d,p,_l1);vec_diff_(v_e,Obs,_l1);
+ if vec_dprod_(v_d,v_a)*vec_dprod_(v_e,v_a)>=0:
+ % |p| and |Obs| are on the same side of the projection plane
+ % |Obs| needs to be recomputed.
+ vec_mult_(v_b,v_b,-1);
+ vec_sum_(Obs,p,v_b);
+ fi;
+ free_vec(v_e);free_vec(v_d);
+ projection_type:=2; % needs to be set before |point_of_view_abs|
+ point_of_view_abs(p,90); % this computes |ObliqueCenter_|
+ % and now, make sure the vectors defining the observer are right:
+ % Create the plane containing lines _l and projection_direction
+ % (defined by point_of_view_abs):
+ new_plane_(xxx_)(1,1,1);
+ def_plane_pl_l_l(xxx_)(_l)(projection_direction);
+ % Compute the angle of |ObsK_| with this plane:
+ obsJangle_=vangle_v_pl_(ObsK_)(xxx_);
+ % rotate |ObsJ_| and |ObsK_| by |obsJangle_| around |ObsI_|
+ vec_rotate_(ObsJ_,ObsI_,obsJangle_);
+ vec_rotate_(ObsK_,ObsI_,obsJangle_);
+ if abs(vangle_v_pl_(ObsK_)(xxx_))>1: % the rotation was done
+ % in the wrong direction
+ vec_rotate_(ObsJ_,ObsI_,-2obsJangle_);
+ vec_rotate_(ObsK_,ObsI_,-2obsJangle_);
+ fi;
+ % |vec_rotate_(ObsJ_,ObsI_,45);| % planometric test
+ % |vec_rotate_(ObsK_,ObsI_,45);| % planometric test
+ free_plane(xxx_);
+ % and now, |ProjJ_| and |ProjK_| must be recomputed:
+ project_vectors(ProjK_,ProjJ_)(ObsK_,ObsJ_)%
+ (projection_direction)(projection_plane);
+ else:
+ message "Error: the ``horizontal plane'' cannot be";
+ message " parallel to the projection plane.";
+ fi;
+ free_line(_l);
+ free_vec(v_c);free_vec(v_b);free_vec(v_a);
+enddef;
+
+% These two are the most common values for the third parameter
+% of |oblique_projection|
+numeric CAVALIER;CAVALIER=45;
+numeric CABINET;CABINET=angle((1,.5)); % atn(.5)
+
+% 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
+
+% 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;
+
+% This array relates an absolute object point number to the
+% absolute point number (that is, to the |vec| array).
+% The absolute object point number is the rank of a point
+% with respect to all object points. The absolute point number
+% considers in addition the extra points, such as |Obs|, which do
+% not belong to an object.
+% If |i| is an absolute object point number, |points_[i]|
+% is the absolute point number.
+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):new_point_(pnt(i));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;
+
+%%\newpage
+% 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)=
+ points_[n+obj_point_offset_(name)]
+ %hide(define_current_point_offset_(name);) pnt(n) % HAS SIDE EFFECTS
+enddef;
+
+% Absolute face number corresponding to local face |n|
+% in object instance |name|
+vardef face_obj(expr name,n)=
+ (n+obj_face_offset_(name))
+ %hide(define_current_face_offset_(name);) face(n) % HAS SIDE EFFECTS
+enddef;
+
+
+% Scale
+numeric drawing_scale;
+drawing_scale=2cm;
+
+% Color
+% This function is useful when a color is expressed in hexadecimal.
+% This does the opposite from |tohexcolor|
+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;
+
+% Convert a color triple into a hexadecimal color string.
+% |rv|, |gv| and |bv| are values between 0 and 1.
+% This does the opposite from |hexcolor|
+vardef tohexcolor(expr rv,gv,bv)=
+ save dig;numeric dig[];
+ hide(
+ dig2=floor(rv*255);dig1=floor((dig2)/16);dig2:=dig2-16*dig1;
+ dig4=floor(gv*255);dig3=floor((dig4)/16);dig4:=dig4-16*dig3;
+ dig6=floor(bv*255);dig5=floor((dig6)/16);dig6:=dig6-16*dig5;
+ for i:=1 upto 6:
+ if dig[i]<10:dig[i]:=dig[i]+48;
+ else:dig[i]:=dig[i]+87;
+ fi;
+ endfor;
+ )
+ char(dig1)&char(dig2)&char(dig3)&char(dig4)&char(dig5)&char(dig6)
+enddef;
+
+% Conversions
+
+% Returns a string encoding the integer |n| as follows:
+% if $n=10*a+b$ with $b<10$,
+% |alphabetize|(|n|)=|alphabetize|(|a|) |&| |char (65+b)|
+% For instance, alphabetize(3835) returns "DIDF"
+% This function is useful in places where digits are not allowed.
+def alphabetize(expr n)=
+ if (n>9):
+ alphabetize(floor(n/10)) & fi
+ char(65+n-10*floor(n/10))
+enddef;
+
+% Filling and contours
+boolean filled_faces,draw_contours;
+filled_faces=true;
+draw_contours=true;
+numeric contour_width; % thickness of contours
+contour_width=1pt;
+color contour_color; % face contours
+contour_color=black;
+
+% 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;
+
+endinput
diff --git a/Master/texmf-dist/metapost/mp3d/3danim.mp b/Master/texmf-dist/metapost/mp3d/3danim.mp
new file mode 100644
index 00000000000..3e49922a69f
--- /dev/null
+++ b/Master/texmf-dist/metapost/mp3d/3danim.mp
@@ -0,0 +1,269 @@
+%%\input epsf
+%%\def\newpage{\vfill\eject}
+%%\advance\vsize1in
+%%\let\ora\overrightarrow
+%%\def\title#1{\hrule\vskip1mm#1\par\vskip1mm\hrule\vskip5mm}
+%%\def\figure#1{\par\centerline{\epsfbox{#1}}}
+%%\title{{\bf 3DANIM.MP: STANDARD ANIMATION DEFINITIONS IN METAPOST}}
+
+%% version 1.34, 17 August 2003
+%% {\bf Denis Roegel} ({\tt roegel@loria.fr})
+%%
+%%This package provides standard animation definitions and
+%%must be used with the {\bf 3d} package.
+
+if known three_d_anim_version:
+ expandafter endinput % avoids loading this package twice
+fi;
+
+% First, we load the 3D package
+input 3d
+
+message "*** 3danim, v1.34 (c) D. Roegel, 17 August 2003 ***";
+numeric three_d_anim_version;
+three_d_anim_version=1.34;
+
+%%\newpage
+%%\title{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;
+
+extra_endfig:=extra_endfig & "compute_bbox;";
+
+boolean show_animation_parameters;
+show_animation_parameters=false;
+% The paper height is purely virtual; it is a paper size
+% such that the full animation can rest on it.
+% It must be compatible with the |PAPERSIZE| option
+% of |gs| in |write_script|.
+% If the size (i.e. the bounding box) of the animation is too large,
+% you may have to change the |paper_height| and the |gs| parameter.
+% The ratio between the paper height and width is assumed to be sqrt(2).
+numeric paper_height;
+paper_height=29.7; % A4 paper height in cm
+numeric output_res;
+output_res=36; % default output resolution of bitmap: 36 dots per inch
+
+% 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,phbp,pwbp;
+ w=xmax_-xmin_;h=ymax_-ymin_;
+ if show_animation_parameters:
+ message "animation bbox: (llx=" & decimal round(xmin_)
+ & ",lly=" & decimal round(ymin_)
+ & ",w=" & decimal round(w) & ",h=" & decimal round(h) & ")";
+ fi;
+ % the lower left corner is put at position (20,20)
+ trx=20-xmin_;
+ try=20-ymin_;
+ if show_animation_parameters:
+ message "translate parameters: "
+ & decimal round(trx) & " " & decimal round(try);
+ fi;
+ xmin_:=xmin_+trx;ymin_:=ymin_+try;
+ delta=10; % extra space
+ phbp=paper_height/2.54*72; % paper height in PostScript points
+ pwbp=phbp/sqrt(2); % paper width in PostScript points
+ % (assuming a sqrt(2) ratio)
+ pnmx=round(xmin_*(output_res/72)-delta);
+ pnmy=round((paper_height/2.54*72-ymin_-h)*(output_res/72)-delta);
+ pnmw=round(w*(output_res/72)+2*delta);
+ pnmh=round(h*(output_res/72)+2*delta);
+ if show_animation_parameters:
+ message "pnmcut parameters (with -r" & decimal output_res & "): "
+ & decimal pnmx & " " & decimal pnmy & " "
+ & decimal pnmw & " " & decimal pnmh;
+ fi;
+ if (pnmx<0) or (pnmy<0) or (pnmx+pnmw>=pwbp) or (pnmy+pnmh>=phbp):
+ message "!! paper size overflow: you may need to reduce the size of";
+ message "!! the animation, or to change the paper size";
+ fi;
+ write_script(round(trx),round(try),
+ pnmx,pnmy,pnmw,pnmh,output_res,jobname,"create_animation.sh");
+enddef;
+
+%%\newpage
+%%\title{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("#! /bin/sh");
+ write_to_file("");
+ write_to_file("/bin/rm -f "&output&".log");
+ write_to_file("for i in `ls "&output&".*| grep '"&output&".[0-9]'`;do");
+ if false: "endfor" fi % indentation hack for meta-mode.el
+ write_to_file("echo $i");
+ write_to_file("echo '=============='");
+ s:="awk < $i '{print} /^%%Page: /{print "&ditto;
+ s:=s&decimal trx&" "&decimal try&" translate\n"&ditto&"}' > $i.ps";
+ write_to_file(s);
+ % ghostscript PostScript into ppm
+ % (the paper size must be compatible with the definition of |paper_height|)
+ s:="gs -sDEVICE=ppmraw -sPAPERSIZE=a4 -dNOPAUSE ";
+ s:=s&"-r"&decimal res &" -sOutputFile=$i.ppm -q -- $i.ps";
+ write_to_file(s);
+ write_to_file("/bin/rm -f $i.ps");
+ % possible alternative:
+ % |s:="mogrify -compress -crop " & decimal(w) & "x" & decimal(h);|
+ % |s:=s&"+"& decimal(xmin) &"+"&decimal(ymin);|
+ % |s:=s&" -colors 32 -format gif $i.ppm";|
+ s:="ppmquant 32 $i.ppm | pnmcut "& decimal(xmin) &" "&decimal(ymin);
+ s:=s&" "&decimal(w)&" "&decimal(h) &" | ";
+ s:=s&"ppmtogif > `expr $i.ppm : '\(.*\)ppm'`gif";
+ write_to_file(s);
+ write_to_file("/bin/rm -f $i.ppm");
+ write_to_file("done");
+ write_to_file("/bin/rm -f "&output&".gif");
+ s:="gifmerge -10 -l1000 ";
+ s:=s&output&".*.gif > "&output&".gif";
+ write_to_file(s);
+ write_to_file("/bin/rm -f "&output&".*.gif");
+ write_to_file(EOF);% end of file
+enddef;
+
+% These definitions produce {\it one\/} image of some kind.
+
+% In the standard animations, the observer follows a circle, shown below:
+%%\figure{vect-fig.17}
+
+% 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,vec_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,vec_I,ang);
+ rotate_obj_pv(name_b,13,vec_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;
+
+%%\newpage
+% 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);
+ new_vec(v_a);
+ vec_def_(v_a,.03*cosd(-a*ang+90),.03*sind(-a*ang+90),0);
+ translate_obj(name_c,v_a);
+ free_vec(v_a);
+ rotate_obj_pv(name_a,1,vec_I,ang);
+ rotate_obj_pv(name_b,13,vec_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,vec_I,ang);
+ rotate_obj_pv(name_b,13,vec_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;
+
+
+%%\newpage
+% 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("obj",name);
+ for i:=imin upto imax:one_image("obj",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("obja",name_a);assign_obj("objb",name_b);
+ translate_obj("objb",vec_K);translate_obj("objb",vec_K);
+ for i:=imin upto imax:
+ one_image_two_objects("obja","objb",i+index,i,10,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("obja",name_a);assign_obj("objb",name_b);
+ assign_obj("objc",name_c);
+ scale_obj("objb",.7);
+ new_vec(v_a);
+ vec_def_vec_(v_a,vec_K);vec_mult_(v_a,v_a,4);put_obj("objb",v_a,1,0,0,0);
+ free_vec(v_a);
+ scale_obj("objc",.5);
+ translate_obj("objc",vec_K);translate_obj("objc",vec_K);
+ for i:=imin upto imax:
+ one_image_three_objects("obja","objb","objc",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("obja",name);assign_obj("objb",name);
+ translate_obj("objb",vec_K);translate_obj("objb",vec_K);
+ for i:=imin upto imax:
+ one_image_two_identical_objects("obja","objb",i+index,i,10,ang);
+ endfor;
+ show_animation_bbox;
+enddef;
+
+endinput
+
diff --git a/Master/texmf-dist/metapost/mp3d/3dgeom.mp b/Master/texmf-dist/metapost/mp3d/3dgeom.mp
new file mode 100644
index 00000000000..22a30315b65
--- /dev/null
+++ b/Master/texmf-dist/metapost/mp3d/3dgeom.mp
@@ -0,0 +1,931 @@
+%%\input epsf
+%%\def\newpage{\vfill\eject}
+%%\advance\vsize1in
+%%\let\ora\overrightarrow
+%%\def\title#1{\hrule\vskip1mm#1\par\vskip1mm\hrule\vskip5mm}
+%%\def\figure#1{\par\centerline{\epsfbox{#1}}}
+%%\title{{\bf 3DGEOM.MP: 3D GEOMETRY IN METAPOST}}
+
+%% version 1.34, 17 August 2003
+%% {\bf Denis Roegel} ({\tt roegel@loria.fr})
+
+% This package provides useful definitions for geometrical drawings.
+% It contains functions dealing with lines, planes, etc.
+
+if known three_d_geom_version:
+ expandafter endinput % avoids loading this package twice
+fi;
+
+% First, we load the 3D package
+input 3d
+% and some utilities
+input 3dutil
+
+message "*** 3dgeom, v1.34 (c) D. Roegel 17 August 2003 ***";
+numeric three_d_geom_version; three_d_geom_version:=1.34;
+
+% WARNING:
+% Known bugs: unnecessary overflows can occur, especially when
+% computing the intersection of two planes.
+
+
+% Among other things, this file defines so-called ``structures.''
+% These structures are different from the ``objects'' manipulated
+% by the main 3d package. For some explanations, see the article
+%
+% Denis Roegel: La géométrie dans l'espace avec METAPOST,
+% Cahiers GUTenberg number 39-40, 2001, pages 107-138.
+% (in French, conference proceedings of GUT2001)
+%
+%
+%
+% Future versions of this module will consider the following structures,
+% not all of which are currently implemented:
+%
+% structure name standard abreviation
+% point p
+% line l
+% plane pl
+% circle c
+% triangle tr
+% sphere s
+% cone co
+% cylinder cy
+% tetrahedron te
+%
+% These names are considered reserved and should not be used for classes.
+%
+% The left column names are used when defining a structure with |def|,
+% |set| or freeing it with |free|.
+%
+% When a function using parameters of these types is defined,
+% the abreviations of the types are part of the function name.
+% For instance, the function giving the intersection between a
+% line and a plane is named |def_inter_l_pl|.
+%
+% Functions computing intersections should be named |def_inter|
+% and should be followed by the resulting type. For instance,
+% the intersection of two lines is |def_inter_p_l_l|,
+% the intersection of two planes is |def_inter_l_pl_pl|
+%
+% Functions computing inscriptions (like a circle inscribed
+% in a triangle) should be named |def_ins|.
+% For instance, |def_ins_c_tr|.
+%
+% Functions computing circumscriptions (like a circle circumscribing
+% a triangle) should be named |def_circums|.
+% For instance, |def_circums_c_tr|.
+%
+% Functions computing exinscriptions (like a circle exinscribed
+% in a triangle) should be named |def_exins|.
+% For instance, |def_exins_c_tr|.
+%
+% Functions computing tangencies (like a tangent to a circle)
+% should be named |def_tang|.
+% For instance, |def_tang_l_c|.
+%
+% Functions computing orthogonal planes, lines, etc. should
+% be named |def_orth|.
+%
+% All these functions can have more parameters than what the name
+% implies.
+%
+% These rules are guidelines, not a standard. If you have some idea
+% on naming conventions, please let me know at roegel@loria.fr.
+%
+% Possibly, more thought should be given in
+% order to distinguish pseudo-objects like ``circle''
+% from the other objects of 3d.mp (like the polyhedra, etc.).
+
+
+% Structures can be allocated, set and freed.
+
+% Our first structure is the line. A line is defined by two points.
+% This is not an object in the usual sense of the 3d package.
+% It is just made of two points.
+% |l| is the line name: it must be different from already known variables
+% |i| and |j| are point numbers
+% (absolute version)
+def new_line_(text l)(expr i,j)=
+ new_points(l)(2);
+ set_line_(l)(i,j);
+enddef;
+
+% The following version takes local point numbers instead of absolute ones.
+def new_line(text l)(expr i,j)=new_line_(l)(pnt(i),pnt(j)) enddef;
+
+% This is used to set a line:
+% (absolute version)
+def set_line_(text l)(expr i,j)=
+ vec_def_vec_(l1,i);
+ vec_def_vec_(l2,j); % l[2]=l[1]+1 (this is assumed elsewhere,
+ % so should never change)
+enddef;
+
+% (local version)
+def set_line(text l)(expr i,j)=set_line_(l)(pnt(i),pnt(j)) enddef;
+
+def free_line(text l)=
+ free_points(l)(2);
+enddef;
+
+% A circle |c| of center |i|, radius |r| and in plane |p|.
+% We store the center as a point, and (r,p[1]) in another point.
+def new_circle(text c)(expr i,r)(text p)=
+ new_points(c)(2);
+ vec_def_vec(c1,i);
+ vec_def(c2,r,p1,0);
+enddef;
+
+% should |set_circle| be defined?
+
+def free_circle(text c)=
+ free_points(c)(2);
+enddef;
+
+% Planes are similar to lines. A plane is just a triple of points.
+% (absolute version)
+def new_plane_(text p)(expr i,j,k)=
+ new_points(p)(3);
+ set_plane_(p)(i,j,k);
+enddef;
+
+% (local version)
+def new_plane(text p)(expr i,j,k)=new_plane_(p)(pnt(i),pnt(j),pnt(k)) enddef;
+
+% (absolute version)
+def set_plane_(text p)(expr i,j,k)=
+ vec_def_vec_(p1,i);
+ vec_def_vec_(p2,j); % p[2]=p[1]+1 (this is assumed elsewhere,
+ % so should never change)
+ vec_def_vec_(p3,k); % p[3]=p[3]+1 (this is assumed elsewhere,
+ % so should never change)
+enddef;
+
+% (local version)
+def set_plane(text p)(expr i,j,k)=set_plane_(p)(pnt(i),pnt(j),pnt(k)) enddef;
+
+def free_plane(text p)=
+ free_points(p)(3);
+enddef;
+
+% Spheres are not yet used, but here is how they will be allocated and freed.
+
+% A sphere is defined with a center |c| and a radius |r|.
+% We store it using two points.
+def new_sphere(text s)(expr c,r)=
+ new_points(s)(2);
+ vec_def_vec(s1,c);
+ vec_def(s2,r,0,0);
+enddef;
+
+% Should |set_sphere| be defined?
+
+def free_sphere(text s)=
+ free_points(s)(2);
+enddef;
+
+% Lines and planes may be used locally or globally to define
+% new points or new lines.
+
+% In order to define a line which is given by a point and a vector,
+% compute a second point before defining the line.
+% In order to define a line which is given by two planes,
+% define the planes and compute the intersection.
+
+% If a plane is given by a parametric equation (1 point, 2 vectors),
+% compute two additional points and define the plane.
+% If a plane is given by an equation ax+by+cz+d=0, compute three
+% points and define the plane.
+
+% Currently, plane equations are not handled separately.
+
+% Projection of a vector |j| on a plane |p|, along a line |l|.
+% The projection, if it exists, is vector |i|.
+% Returns |true| is there is a projection, and |false| if there is none.
+vardef proj_v_v_l_pl_(expr i,j)(text l)(text p)=
+ save pa,pb,int; boolean int;
+ hide(
+ new_point(pa);new_point(pb);
+ % we project two points: the origin, and origin+v(j):
+ if def_proj_pl_(pa)(p)(point_null)(l):
+ if def_proj_pl_(pb)(p)(j)(l):
+ vec_diff_(i,pb,pa);
+ int=true;
+ else:
+ message "Second point can not be projected";
+ int=false;
+ fi;
+ else: int=false;
+ message "Origin can not be projected";
+ fi;
+ free_point(pb);free_point(pa);
+ )
+ int
+enddef;
+
+% The next function checks if a point is part of a plane.
+% Returns |true| is point |i| is in the plane |p|.
+vardef point_in_plane_p_pl_(expr i)(text p)=
+ save v_a;boolean res;
+ hide(
+ new_vec(v_a);new_vec(v_b);
+ def_normal_p_(v_a)(p);
+ vec_diff_(v_b,p1,i);
+ if vec_dprod_(v_a,v_b)=0: res=true;else: res=false;fi;
+ free_vec(v_b);free_vec(v_a);
+ )
+ res
+enddef;
+
+% The next function finds the angle of a vector with respect to a plane.
+% Returns the angle of a vector |v| with respect to a plane |p|.
+vardef vangle_v_pl_(expr v)(text p)=
+ save v_a,an_;
+ hide(
+ new_vec(v_a);
+ % we compute a vector normal to the plane:
+ def_normal_p_(v_a)(p);
+ an_=90-vangle_v_v_(v,v_a);
+ free_vec(v_a);
+ )
+ an_
+enddef;
+
+% Compute the angle between two vectors
+% The angle is always between 0 and 180,
+% since this is the best one can do with two vectors.
+% If we had a third vector, we could be more accurate.
+vardef vangle_v_v_(expr va,vb)=
+ save cosa_,sina_;
+ hide(
+ cosa_=vec_dprod_(va,vb)/vec_mod_(va)/vec_mod_(vb);
+ if cosa_>1: % sometimes, this happens with rounding errors
+ sina_=0;
+ else:
+ sina_= 1 +-+ cosa_; % sqrt(1-cosa_**2)
+ fi;
+ )
+ angle((cosa_,sina_))
+enddef;
+
+% Define a plane with two lines:
+def def_plane_pl_l_l(text p)(text l)(text m)=
+ set_plane_(p)(l1,l2,1); % the last value is irrelevant
+ vec_diff_(p3,m2,m1);vec_sum_(p3,p3,l1);
+enddef;
+
+% Define the plane orthogonal to a line and going through a point
+% (not necessarily belonging to the plane):
+% the plane must already have been defined
+% |p|=plane, |l|=line, |i|=point
+%...
+% (absolute version)
+vardef def_orth_pl_l_p_(text p)(text l)(expr i)=
+ new_vec(va);new_vec(vb);new_vec(vc);new_vec(h);
+ vec_def_vec_(p1,i); % this is the first point of the plane
+ vec_diff_(va,l2,l1);
+ vec_def_vec_(vb,i);
+ if abs(xval(va))<absmin(yval(va),zval(va)):
+ vec_sum_(vb,vb,vec_I);
+ elseif abs(yval(va))<absmin(xval(va),zval(va)):
+ vec_sum_(vb,vb,vec_J);
+ else:vec_sum_(vb,vb,vec_K);
+ fi;
+ % now, |vb| is a point not on the line and not too close to it
+ % we compute a vertical to the line
+ def_vert_l_(h,vb)(l);
+ vec_diff_(vb,vb,h);vec_unit_(vb,vb);
+ vec_sum_(p2,vb,p1);
+ % |p[2]| is now a point of the plane
+ % a third point is obtained by cross product
+ vec_prod_(vc,va,vb);vec_unit_(vc,vc);vec_sum_(p3,vc,p1);
+ free_vec(h);free_vec(vc);free_vec(vb);free_vec(va);
+enddef;
+
+% (local version)
+vardef def_orth_pl_l_p(text p)(text l)(expr i)=
+ def_orth_pl_l_p_(p)(l)(pnt(i))
+enddef;
+
+% Line orthogonal to a plane and going through a point
+% (not necessarily belonging to the plane);
+% from the three points defining the plane, compute a normal,
+% and add it to the point, this gives a second point,
+% and make a line out of it
+% (absolute version)
+vardef def_orth_l_pl_p_(text l)(text p)(expr i)=
+ new_vec(va);new_vec(vb);
+ vec_def_vec_(l1,i);
+ def_normal_p_(l2)(p);
+ vec_sum_(l2,l2,l1);
+ free_vec(vb);free_vec(va);
+enddef;
+
+% (local version)
+vardef def_orth_l_pl_p(text l)(text p)(expr i)=
+ def_orth_l_pl_p_(l)(p)(pnt(i))
+enddef;
+
+% Unitary vector normal to a plane.
+% |v| is a vector that must have been defined
+% (absolute version)
+vardef def_normal_p_(expr v)(text p)=
+ new_vec(va);new_vec(vb);
+ vec_diff_(va,p2,p1);vec_diff_(vb,p3,p1);vec_prod_(v,va,vb);
+ vec_unit_(v,v);
+ free_vec(vb);free_vec(va);
+enddef;
+
+% Unitary vector normal to a plane (local version)
+vardef def_normal_p(expr v)(text p)=def_normal_p_(pnt(v))(p) enddef;
+
+% The following two functions are old versions of the
+% line/plane intersection. They are not used anymore.
+%
+% Intersection line/plane
+% Point |i| is the intersection
+% The return value is |true| if the intersection is a point,
+% |false| otherwise
+% (absolute version)
+vardef old_def_inter_p_l_pl_(expr i)(text l)(text p)=
+ save d,t,int;boolean int;
+ hide(
+ new_vec(va);new_vec(vb);new_vec(vc);
+ % first, we compute a vector normal to the plane
+ vec_diff_(va,p2,p1);
+ vec_diff_(vb,p3,p2);
+ vec_prod_(vc,va,vb);
+ % we want the plane equation as ax+by+cz+d=0
+ % the normal vector gives us (a,b,c)
+ % d is then easy to compute
+ d=-xval(vc)*xval(p1)-yval(vc)*yval(p1)-zval(vc)*zval(p1);
+ vec_diff_(i,l2,l1);
+ if vec_dprod_(i,vc)=0: % the line is parallel to the plane
+ int:=false;
+ else:
+ int:=true;
+ t=-(d+xval(vc)*xval(l1)+yval(vc)*yval(l1)+zval(vc)*zval(l1))
+ /vec_dprod_(i,vc);
+ vec_mult_(i,i,t);vec_sum_(i,i,l1);
+ fi;
+ free_vec(vc);free_vec(vb);free_vec(va);
+ )
+ int
+enddef;
+
+% same (local version)
+vardef old_def_inter_p_l_pl(expr i)(text l)(text p)=
+ def_inter_p_l_pl_(pnt(i))(l)(p)
+enddef;
+
+% Intersection line/plane (absolute version)
+% Point |i| is the intersection.
+% The return value is |true| if the intersection is a point,
+% |false| otherwise
+vardef def_inter_p_l_pl_(expr i)(text l)(text p)=
+ save int;boolean int;
+ hide(
+ new_points(loc)(3);
+ vec_diff_(loc1,p2,p1);vec_diff_(loc2,p3,p1);vec_prod_(loc3,loc1,loc2);
+ vec_diff_(loc1,p1,l1);vec_diff_(loc2,l2,l1);
+ if vec_dprod_(loc2,loc3)<>0:
+ vec_mult_(loc2,loc2,vec_dprod_(loc1,loc3)/vec_dprod_(loc2,loc3));
+ vec_sum_(i,l1,loc2);
+ int:=true;
+ % Remark: in order to prove that point |i| is on the plane, it
+ % suffices to compute vec(ci).(vec(cd) /\ vec(ce))
+ % =(-vec(ac)+vec(ai)).(vec(cd) /\ vec(ce))
+ % =-vec(ac).(vec(cd) /\ vec(ce))
+ % +(vec(ab).(vec(cd) /\ vec(ce))) vec(ac).(vec(cd) /\ vec(ce))
+ % ----------------------------
+ % vec(ab).(vec(cd) /\ vec(ce))
+ % =0
+ else: % the line is parallel to the plane
+ int:=false;
+ fi;
+ free_points(loc)(3);
+ )
+ int
+enddef;
+
+% Intersection line/plane (local version)
+vardef def_inter_p_l_pl(expr i)(text l)(text p)=
+ def_inter_p_l_pl_(pnt(i))(l)(p)
+enddef;
+
+% The following function is used in |def_inter_l_pl_pl|.
+% We could simplify it by breaking it in two.
+vardef def_inter_l_pl_pl_base_case_(text l)(expr pa,pb,pc)(text q)=
+ save trial;
+ new_line_(trial)(pa,pb);
+ if def_inter_p_l_pl_(l1)(trial)(q):
+ else: % there is no intersection or the intersection is the line
+ vec_def_vec_(trial1,pa);
+ mid_point_(trial2,pb,pc);
+ if def_inter_p_l_pl_(l1)(trial)(q):
+ else:
+ message "THIS SHOULD NOT HAPPEN, PLEASE REPORT THIS PROBLEM";
+ fi;
+ fi;
+ set_line_(trial)(pa,pc);
+ if def_inter_p_l_pl_(l2)(trial)(q):
+ else: % there is no intersection or the intersection is the line
+ vec_def_vec_(trial1,pa);
+ mid_point_(trial2,pb,pc);
+ if def_inter_p_l_pl_(l2)(trial)(q):
+ else:
+ message "THIS SHOULD NOT HAPPEN, PLEASE REPORT THIS PROBLEM";
+ fi;
+ fi;
+ free_line(trial);
+enddef;
+
+% Intersection of two planes.
+% TO DO: this function is not yet robust enough, because
+% unnecessary overflows can occur.
+% A boolean is set if there is no intersection.
+% The line |l| must already have been defined.
+vardef def_inter_l_pl_pl(text l)(text p)(text q)=
+ save trial,da,db,dc,int;boolean int;
+ hide(
+ % we first search the point of p1, p2, p3 which is the farthest
+ % from q;
+ da=dist_pl_(p1)(q);db=dist_pl_(p2)(q);dc=dist_pl_(p3)(q);
+ if (da=db) and (db=dc): % the two planes are parallel
+ int:=false;
+ else:
+ int:=true;
+ if (da>=db) and (da>=dc):
+ def_inter_l_pl_pl_base_case_(l)(p1,p2,p3)(q);
+ elseif (db>=da) and (db>=dc):
+ def_inter_l_pl_pl_base_case_(l)(p2,p1,p3)(q);
+ else:
+ def_inter_l_pl_pl_base_case_(l)(p3,p1,p2)(q);
+ fi;
+ fi;
+ )
+ int
+enddef;
+
+% Visual intersection between lines (jk) and (lm).
+% The computed intersection lies on (jk).
+% Returns true if there is an intersection, false otherwise.
+% (absolute version)
+vardef def_visual_inter_(expr i)(expr j,k,l,m)=
+ save pla,plb,la,lb,d,int;boolean int;
+ hide(
+ new_plane_(pla)(Obs,l,m);new_plane_(plb)(Obs,j,k);
+ new_line_(la)(0,0);new_line_(lb)(j,k);
+ if def_inter_l_pl_pl(la)(pla)(plb):
+ int:=true;
+ % |d| is the closest distance between lines |la| and |lb|
+ % We don't use |d| here, and are only interested in point |i|.
+ d=def_inter_p_l_l_(i)(la)(lb);
+ else:
+ int:=false;
+ fi;
+ free_line(lb);free_line(la);free_plane(plb);free_plane(pla);
+ ) int
+enddef;
+
+% same (local version)
+vardef def_visual_inter(expr i)(expr j,k,l,m)=
+ def_visual_inter_(pnt(i),pnt(j),pnt(k),pnt(l),pnt(m))
+enddef;
+
+% Point of a line at a given distance from a given point.
+% |i| = new point |d|=distance |j|=point |l|=line
+% $|d|>0$ or $|d|<0$ give two different points.
+% If there is an intersection, the function returns |true|;
+% otherwise it returns |false|.
+% (absolute version)
+vardef def_point_at_(expr i)(expr d)(expr j)(text l)=
+ save dj,ld,int;boolean int;
+ hide(
+ new_point(h);new_point(hc);
+ def_vert_l_(h,j)(l);
+ vec_diff_(hc,j,h);
+ if d*d-vec_dprod_(hc,hc)>=0: int:=true;
+ ld=sign(d)*sqrt(d*d-vec_dprod_(hc,hc));
+ vec_diff_(i,l1,l2);
+ vec_unit_(i,i);
+ vec_mult_(i,i,ld);
+ vec_sum_(i,i,h);
+ else: int:=false;
+ fi;
+ free_point(hc);
+ free_point(h);
+ )
+ int
+enddef;
+
+% same (local version)
+vardef def_point_at(expr i)(expr d)(expr j)(text l)=
+ def_point_at_(pnt(i))(d)(pnt(j))(l)
+enddef;
+
+% Define a vertical of a line.
+% Point |i| is obtained as the intersection of a vertical
+% starting from point |j| and reaching the line |l|.
+vardef def_vert_l_(expr i,j)(text l)=
+ new_points(loc)(3);
+ vec_diff_(loc1,j,l1);vec_diff_(loc2,l2,l1);
+ vec_mult_(loc3,loc2,vec_dprod_(loc1,loc2)/vec_dprod_(loc2,loc2));
+ vec_sum_(i,loc3,l1);
+ free_points(loc)(3);
+enddef;
+
+% Define a vertical. (local version)
+vardef def_vert_l(expr i,j)(text l)=
+ def_vert_l_(pnt(i),pnt(j))(l);
+enddef;
+
+% Vertical falling on a plane.
+% Point |j| falls on plane |p| at point |i| (absolute version)
+vardef def_vert_pl_(expr i)(expr j)(text p)=
+ save d;
+ new_vec(va);new_vec(vb);
+ def_normal_p_(va)(p);
+ vec_diff_(vb,j,p1);
+ d=-vec_dprod_(vb,va);
+ vec_mult_(va,va,d);
+ vec_sum_(vb,vb,va);
+ vec_sum_(i,p1,vb);
+ free_vec(vb);free_vec(va);
+enddef;
+
+% same (local version)
+vardef def_vert_pl(expr i)(expr j)(text p)=
+ def_vert_pl_(pnt(i))(pnt(j))(p)
+enddef;
+
+% Distance to a plane.
+% (absolute version)
+vardef dist_pl_(expr i)(text p)=
+ save d;
+ hide(
+ new_vec(va);
+ def_vert_pl_(va)(i)(p);
+ vec_diff_(va,va,i);
+ d=vec_mod_(va);
+ free_vec(va);
+ )
+ d
+enddef;
+
+% (local version)
+def dist_pl(expr i)(text p)=dist_pl_(pnt(i))(p) enddef;
+
+% Projections on planes or lines, according to a direction.
+% This one is very hazardous: use epsilon
+% Find point |i| on |l| from point |j| using direction |d|
+
+def def_proj_l_(expr i)(text l)(expr j)(text d)=
+ NOT YET IMPLEMENTED
+enddef;
+
+def def_proj_l(expr i)(text l)(expr j)(text d)=
+ def_proj_l_(pnt(i))(l)(pnt(j))(d)
+enddef;
+
+% Find point |i| on |p| from point |j| using direction |d|.
+vardef def_proj_pl_(expr i)(text p)(expr j)(text d)=
+ save l_,int; boolean int;
+ hide(
+ % we compute the intersection between line (|j|+|d|) and plane |p|
+ new_line_(l_)(1,1); % we must take a name that cannot
+ % conflict with the text replacement of |d|
+ vec_diff_(l_2,d2,d1);vec_sum_(l_2,l_2,j);
+ vec_def_vec_(l_1,j);
+ if def_inter_p_l_pl_(i)(l_)(p):int=true;
+ else: int=false;
+ fi;
+ free_line(l_);
+ )
+ int
+enddef;
+
+def def_proj_pl(expr i)(text p)(expr j)(text d)=
+ def_proj_pl_(pnt(i))(p)(pnt(j))(d)
+enddef;
+
+% Central projection on a plane.
+def def_cproj_pl_(expr i)(text p)(expr j)(expr k)=
+% use |def_proj_p|
+ NOT YET IMPLEMENTED
+enddef;
+
+% Central projection on a plane.
+def def_cproj_pl(expr i)(text p)(expr j)(expr k)=
+ def_cproj_pl_(pnt(i))(p)(pnt(j))(pnt(k))
+enddef;
+
+
+% Intersection of two lines (hazardous).
+% Due to rounding errors, two lines that should intersect
+% may not do so in reality. Therefore,
+% we compute the point which is the middle of the two
+% closest points between the lines and return the distance
+% between the two lines. If the lines are parallel (possibly
+% identical), we return -1.
+vardef def_inter_p_l_l_(expr i)(text l)(text m)=
+ save ga,gb,gc,gd,ge,gf,t,u,d,mx;
+ hide(
+ new_point(va);new_point(vb);new_point(vc);new_point(h);new_point(k);
+ vec_diff_(va,m1,l1);
+ vec_diff_(vb,l2,l1);
+ vec_diff_(vc,m2,m1);
+ ga=vec_dprod_(vc,vb);gb=-vec_dprod_(vb,vb);
+ gc=vec_dprod_(va,vb);gd=vec_dprod_(vc,vc);
+ ge=-ga;gf=vec_dprod_(va,vc);
+ % compute the max of ga,gb,...
+ mx:=absmax(ga,gb);mx:=absmax(mx,gc);mx:=absmax(mx,gd);mx:=absmax(mx,ge);
+ mx:=absmax(mx,gf);
+ ga:=ga/mx;gb:=gb/mx;gc:=gc/mx;gd:=gd/mx;ge:=ge/mx;gf:=gf/mx;
+ if ga*ge=gb*gd: % the lines are parallel
+ % we return -1
+ d=-1;
+ else:
+ t=(gc*gd-ga*gf)/(ga*ge-gb*gd);u=(gb*gf-gc*ge)/(ga*ge-gb*gd);
+ vec_diff_(h,l2,l1);vec_mult_(h,h,t);vec_sum_(h,h,l1);
+ vec_diff_(k,m2,m1);vec_mult_(k,k,u);vec_sum_(k,k,m1);
+ % |h| and |k| are now the closest points
+ % we set |i| to the middle of |h| and |k| and return the distance |hk|
+ mid_point_(i,h,k);
+ vec_diff_(h,h,k);d=vec_mod_(h);
+ fi;
+ free_point(k);free_point(h);free_point(vc);free_point(vb);free_point(va);
+ )
+ d
+enddef;
+
+def def_inter_p_l_l(expr i)(text l)(text m)=
+ def_inter_p_l_l_(pnt(i))(l)(m)
+enddef;
+
+% Find point |i| symmetric of point |j| with respect to point |k|
+def def_sym_(expr i)(expr j)(expr k)=
+ NOT YET IMPLEMENTED
+enddef;
+
+def def_sym(expr i)(expr j)(expr k)=
+ def_sym_(pnt(i))(pnt(j))(pnt(k))
+enddef;
+
+% Find point |i| symmetric of point |j| with respect to plane |p|
+def def_sym_pl_(expr i)(expr j)(text p)=
+ NOT YET IMPLEMENTED
+enddef;
+
+def def_sym_pl(expr i)(expr j)(text p)=
+ def_sym_pl_(pnt(i))(pnt(j))(p)
+enddef;
+
+% Find point |i| symmetric of point |j| with respect to line |l|.
+% That's a mere 180 degrees rotation around the line.
+def def_sym_l_(expr i)(expr j)(text l)=
+ NOT YET IMPLEMENTED
+enddef;
+
+def def_sym_l(expr i)(expr j)(text l)=
+ def_sym_l_(pnt(i))(pnt(j))(l)
+enddef;
+
+
+% Intersection circle/line (hazardous).
+% If some intersection does not exist, |infty| is put for its values
+def def_inter_p_p_c_l_(expr i,j)(text c)(text l)=
+ NOT YET IMPLEMENTED
+enddef;
+
+def def_inter_p_p_c_l(expr i,j)(text c)(text l)=
+ def_inter_p_p_c_l_(pnt(i),pnt(j))(c)(l)
+enddef;
+
+% circle/plane
+% A similar coding will distinguish the four cases:
+% one point, two points, the full circle, nothing
+def def_inter_p_p_c_pl_(expr i,j)(text c)(text p)=
+ NOT YET IMPLEMENTED
+enddef;
+
+def def_inter_p_p_c_pl(expr i,j)(text c)(text p)=
+ def_inter_p_p_c_pl_(pnt(i),pnt(j))(c)(p)
+enddef;
+
+% circle/circle
+% A similar coding will distinguish the four cases:
+% one point, two points, the full circle, nothing
+def def_inter_p_p_c_c_(expr i,j)(text ca)(text cb)=
+ NOT YET IMPLEMENTED
+enddef;
+
+def def_inter_p_p_c_c(expr i,j)(text ca)(text cb)=
+ def_inter_p_p_c_c_(pnt(i),pnt(j))(ca)(cb)
+enddef;
+
+% Computation of tangent lines and planes.
+
+% Tangent line to a circle at a given point.
+def def_tang_l_c_p_(text l)(text c)(expr i)=
+ NOT YET IMPLEMENTED
+enddef;
+
+def def_tang_l_c_p(text l)(text c)(expr i)=
+ def_tang_l_c_p_(l)(c)(pnt(i))
+enddef;
+
+% Tangent plane to a sphere at a given point.
+def def_tang_pl_s_p_(text p)(text s)(expr i)=
+ NOT YET IMPLEMENTED
+enddef;
+
+def def_tang_pl_s_p(text p)(text s)(expr i)=
+ def_tang_pl_s_p_(p)(s)(pnt(i))
+enddef;
+
+% Sphere defined by four non-coplanar points.
+def def_sphere_through_(text s)(expr i,j,k,l)=
+ NOT YET IMPLEMENTED
+enddef;
+
+def def_sphere_through(text s)(expr i,j,k,l)=
+ def_sphere_through_(s)(pnt(i),pnt(j),pnt(k),pnt(l))
+enddef;
+
+% Line going through a point and parallel to another line.
+def def_parallel_l_p_pl_(text l)(expr i)(text m)=
+ NOT YET IMPLEMENTED
+enddef;
+
+def def_parallel_l_p_pl(text l)(expr i)(text m)=
+ def_parallel_l_p_pl_(l)(pnt(i))(m)
+enddef;
+
+% Plane going through a point and parallel to another plane.
+def def_parallel_pl_p_pl_(text p)(expr i)(text q)=
+ NOT YET IMPLEMENTED
+enddef;
+
+def def_parallel_pl_p_pl(text p)(expr i)(text q)=
+ def_parallel_pl_p_pl_(p)(pnt(i))(q)
+enddef;
+
+def def_rectangle_one_side_(expr p)(text l)(text pa)(text pb)(text pc)=
+ if def_inter_l_pl_pl(l)(pb)(pc):
+ else:
+ message "YOUR PLANES ARE NOT WELL SPECIFIED 1";
+ fi;
+ if def_inter_p_l_pl_(p)(l)(pa):
+ else:
+ message "YOUR PLANES ARE NOT WELL SPECIFIED 2";
+ fi;
+enddef;
+
+% A rectangle (for instance representing a plane) can be defined
+% from five planes; the rectangle is made of four points (corners)
+% |pa| is the plane containing the rectangle
+vardef def_rectangle_pl_pl_pl_pl_pl_(expr ca,cb,cc,cd)
+ (text pa)(text pb)(text pc)(text pd)(text pe)=
+ save l;
+ new_line_(l)(1,1);
+ def_rectangle_one_side_(ca)(l)(pa)(pb)(pc);
+ def_rectangle_one_side_(cb)(l)(pa)(pc)(pd);
+ def_rectangle_one_side_(cc)(l)(pa)(pd)(pe);
+ def_rectangle_one_side_(cd)(l)(pa)(pe)(pb);
+ free_line(l);
+enddef;
+
+% Instead of using four additional planes, one can also use eight points:
+% the order of the point is important.
+vardef def_rectangle_pl_(expr ca,cb,cc,cd)
+ (text pa)(expr pta,ptb,ptc,ptd,pte,ptf,ptg,pth)=
+ save pb,pc,pd,pe;
+ % we create the four additionnal planes
+ new_plane_(pb)(pta,ptb,pte);new_plane_(pc)(ptb,ptc,ptf);
+ new_plane_(pd)(ptc,ptd,ptg);new_plane_(pe)(ptd,pta,pth);
+ def_rectangle_pl_pl_pl_pl_pl_(ca,cb,cc,cd)(pa)(pb)(pc)(pd)(pe);
+ free_plane(pe);free_plane(pd);free_plane(pc);free_plane(pb);
+enddef;
+
+def draw_rectangle(expr i,j,k,l)=
+ draw_line(i,j);draw_line(j,k);draw_line(k,l);draw_line(l,i);
+enddef;
+
+numeric mark_h,mark_l;mark_h=2mm;mark_l=1mm;
+
+def draw_one_mark(expr p,a)=
+ draw (p+unitvector(dir(a))*mark_h/2)--(p-unitvector(dir(a))*mark_h/2);
+enddef;
+
+% Draw |n| marks between points |i| and |j|.
+% |i| and |j| are local points and there is no absolute version
+% since this is a drawing function.
+vardef draw_equal_marks(expr i,j,n)=
+ save a,k,l,start;
+ a=angle(z[ipnt_(j)]-z[ipnt_(i)])+90;
+ l=(x[ipnt_(j)]-x[ipnt_(i)])++(y[ipnt_(j)]-y[ipnt_(i)]);
+ if n=1:
+ draw_one_mark(.5[z[ipnt_(i)],z[ipnt_(j)]],a);
+ elseif n>1:
+ start=0.5-(n-1)*mark_l/(2*l);
+ for k:=0 upto n-1:
+ draw_one_mark((start+k*mark_l/l)[z[ipnt_(i)],z[ipnt_(j)]],a);
+ endfor;
+ else: message "parameter " & decimal n & " should be positive";
+ fi;
+enddef;
+
+numeric square_angle_size;
+square_angle_size=0.2;
+
+% (absolute version)
+def def_right_angle_(expr pi,pj,pk,i,j,k)=
+ vec_diff_(pj,j,i);vec_diff_(pk,k,i);
+ if vec_mod_(pj)>0:
+ vec_mult_(pj,pj,square_angle_size/vec_mod_(pj));
+ fi;
+ if vec_mod_(pk)>0:
+ vec_mult_(pk,pk,square_angle_size/vec_mod_(pk));
+ fi;
+ vec_sum_(pi,i,pj);vec_sum_(pi,pi,pk);
+ vec_sum_(pj,pj,i);vec_sum_(pk,pk,i);
+enddef;
+
+% (local version)
+def def_right_angle(expr pi,pj,pk,i,j,k)=
+ def_right_angle_(pnt(pi),pnt(pj),pnt(pk),pnt(i),pnt(j),pnt(k));
+enddef;
+
+% Right angle on a plane projection.
+% Similar to |def_right_angle_|.
+% This also defines the vertical projection as |vp|.
+vardef def_right_angle_p_(expr pi,pj,pk,vp)(expr i)(text p)=
+ def_vert_pl_(vp)(i)(p);
+ new_vec(va);
+ vec_diff_(va,p1,p2);
+ vec_sum_(va,va,vp); % va is now a second point on the plane,
+ % different from the projection
+ def_right_angle_(pi,pj,pk,vp,va,i);
+ free_vec(va);
+enddef;
+
+def draw_right_angle(expr pi,pj,pk)=
+ draw z[ipnt_(pj)]--z[ipnt_(pi)]--z[ipnt_(pk)];
+enddef;
+
+def draw_double_right_angle(expr pi,pj,pk,pl)=
+ draw z[ipnt_(pj)]--z[ipnt_(pi)]--z[ipnt_(pk)]--z[ipnt_(pl)]--cycle;
+enddef;
+
+% |draw_line| with extra drawing in either directions
+def draw_line_extra(expr i,j)(expr exi,exj)=
+ draw exi[z[ipnt_(i)],z[ipnt_(j)]]--exj[z[ipnt_(i)],z[ipnt_(j)]];
+enddef;
+
+% defines point |i| at position |t| on segment |a|-|b| (absolute version)
+def set_extra_point_(expr i,a,b,t)=
+ vec_diff_(i,b,a);vec_mult_(i,i,t);vec_sum_(i,i,a);
+enddef;
+
+% defines point |i| at position |t| on segment |a|-|b| (local version)
+def set_extra_point(expr i,a,b,t)=
+ set_extra_point_(pnt(i),pnt(a),pnt(b),t);
+enddef;
+
+% labels with local points
+vardef thelabel_obj@#(expr s,n) =
+ thelabel.@#(s,z[ipnt_(n)])
+enddef;
+
+def label_obj = draw thelabel_obj enddef;
+
+% The plane |p| (which must have been initialized) is defined
+% as the screen plane. This is useful for computing vanishing points
+def def_screen_pl(text p)=
+ vec_mult_(p1,ObsI_,Obs_dist);vec_sum_(p1,p1,Obs); % center of screen
+ vec_sum_(p2,p1,ObsJ_);vec_sum_(p3,p1,ObsK_);
+enddef;
+
+% |i| is the resulting point, |l| defines a line in space,
+% |s| is the screen plane
+% Returns |true| is there is a vanishing point, otherwise |false|.
+vardef def_vanishing_point_p_l_pl_(expr i)(text l)(text s)=
+ save vp;boolean vp;
+ hide(
+ new_vec(v);
+ vec_diff_(v,l2,l1);vec_sum_(v,Obs,v);
+ new_line_(obsl)(Obs,v);
+ if def_inter_p_l_pl_(i)(obsl)(s):vp=true;else:vp=false;fi;
+ free_line(obsl);
+ free_vec(v);
+ )
+ vp
+enddef;
+
+def def_vanishing_point_p_l_pl(expr i)(text l)(text s)=
+ def_vanishing_point_p_l_pl_(pnt(i))(l)(s)
+enddef;
+
+endinput
diff --git a/Master/texmf-dist/metapost/mp3d/3dpoly.mp b/Master/texmf-dist/metapost/mp3d/3dpoly.mp
new file mode 100644
index 00000000000..5d3b8f5f2fa
--- /dev/null
+++ b/Master/texmf-dist/metapost/mp3d/3dpoly.mp
@@ -0,0 +1,294 @@
+%%\input epsf
+%%\def\newpage{\vfill\eject}
+%%\def\vc#1{$\vcenter{#1}$}
+%%\def\title#1{\hrule\vskip1mm#1\par\vskip1mm\hrule\vskip5mm}
+%%\def\figure#1{\par\centerline{\epsfbox{#1}}}
+%%\title{{\bf 3DPOLY.MP: DEFINITION OF REGULAR CONVEX POLYHEDRA IN METAPOST}}
+
+%% version 1.34, 17 August 2003
+%% {\bf Denis Roegel} ({\tt roegel@loria.fr})
+%%
+%%This package defines the five regular convex polyhedra classes and
+%%must be used with the {\bf 3d} package.
+
+if known three_d_poly_version:
+ expandafter endinput % avoids loading this package twice
+fi;
+
+message "*** 3dpoly, v1.34 (c) D. Roegel, 17 August 2003 ***";
+numeric three_d_poly_version;
+three_d_poly_version=1.34;
+
+%%Each object definition is parameterized by
+%%an instance identification.
+
+%%All polyhedra are normalized (centered on the origin and inscriptible
+%%in a sphere of radius 1).
+
+%%The polyhedra are defined with a certain number of points and faces.
+%%The points comprise both the vertices and the center of the object.
+
+%%It might be useful to remind the Euler formula linking
+%%the number of faces $f$,
+%%vertices $v$ and edges $e$ in polyhedra with no hole in dimension 3:
+%%$f+v-e=2$. This can be verified in the following table:
+
+%%\centerline{\vbox{\halign{\quad#\hfil\quad&&\quad\hfil#\quad\cr
+%% \omit \hfil type\hfil&\omit\hfil faces\hfil
+%% &\omit\hfil vertices\hfil&\omit\hfil edges\hfil\cr
+%% tetrahedron&4&4&6\cr
+%% cube&6&8&12\cr
+%% octahedron&8&6&12\cr
+%% dodecahedron&12&20&30\cr
+%% icosahedron&20&12&30\cr}}}
+
+%%For more details, see the documentation enclosed in the distribution.
+
+%%\newpage\title{TETRAHEDRON (4 faces)\kern1cm\epsfbox{tetra.ps}}
+
+%% Construction of a tetrahedron:
+
+%% One face is \vc{\epsfbox{vect-fig.1}} with $h=c\sqrt{3}/2$
+
+%% The angle $\alpha$ between two faces is given by
+%% $\sin(\alpha/2)={\displaystyle c/2\over \displaystyle c\sqrt{3}/2}=1/\sqrt3$
+
+%% \figure{vect-fig.2}
+
+def set_tetrahedron_points(expr inst)=
+ set_point(1)(0,0,0);set_point(2)(1,0,0);
+ set_point(3)(cosd(60),sind(60),0);% $60=360/6$
+ sinan=1/sqrt(3);cosan=sqrt(1-sinan**2);
+ an=180-2*angle((cosan,sinan));
+ new_face_point(4,1,2,3,an);
+ normalize_obj(inst)(1,2,3,4);
+ set_point(5)(0,0,0);% center of tetrahedron
+enddef;
+
+vardef def_tetrahedron(expr inst)=
+ new_obj_points(inst,5);% 4 vertices and the center
+ new_obj_faces(inst,4);% 4 faces in the definition
+ set_tetrahedron_points(inst);
+ set_obj_face(1,"1,2,4","b4fefe");
+ set_obj_face(2,"2,3,4","b49bc0");
+ set_obj_face(3,"1,4,3","b4c8fe");
+ set_obj_face(4,"1,3,2","b4fe40");
+enddef;
+
+%%\newpage\title{CUBE (HEXAHEDRON) (6 faces)\kern1cm\epsfbox{cube.ps}}
+
+%%The definition of a cube is straightforward.
+
+def set_cube_points(expr inst)=
+ set_point(1)(0,0,0);
+ set_point(2)(1,0,0);
+ set_point(3)(1,1,0);
+ set_point(4)(0,1,0);
+ set_point(5)(0,0,1);
+ set_point(6)(1,0,1);
+ set_point(7)(1,1,1);
+ set_point(8)(0,1,1);
+ normalize_obj(inst)(1,7);% 1 and 7 are opposite vertices
+ set_point(9)(0,0,0);% center of cube
+enddef;
+
+vardef def_cube(expr inst)=
+ new_obj_points(inst,9);% 8 vertices and the center
+ new_obj_faces(inst,6);% 6 faces in the definition
+ set_cube_points(inst);
+ set_obj_face(1,"1,2,6,5","b4fe40");
+ set_obj_face(2,"2,3,7,6","45d040");
+ set_obj_face(3,"4,8,7,3","45a114");
+ set_obj_face(4,"1,5,8,4","45a1d4");
+ set_obj_face(5,"5,6,7,8","4569d4");
+ set_obj_face(6,"4,3,2,1","112da1");
+enddef;
+
+
+%%\newpage\title{OCTAHEDRON (8 faces)\kern1cm\epsfbox{octa.ps}}
+
+%% A section of a pyramid is: \vc{\epsfbox{vect-fig.3}}
+
+%% The height of this pyramid is given
+%% by $H^2={(c\sqrt3/2)}^2-{(c/2)}^2=c^2/2$, hence $H=c/\sqrt2$
+
+def set_octahedron_points(expr inst)=
+ set_point(1)(0,0,0);
+ set_point(2)(1,0,0);
+ set_point(3)(1,1,0);
+ set_point(4)(0,1,0);
+ set_point(5)(.5,.5,1/sqrt(2));
+ set_point(6)(.5,.5,-1/sqrt(2));
+ normalize_obj(inst)(5,6);% 5 and 6 are opposite vertices
+ set_point(7)(0,0,0);% center of octahedron
+enddef;
+
+vardef def_octahedron(expr inst)=
+ new_obj_points(inst,7);% 6 vertices and the center
+ new_obj_faces(inst,8);% 8 faces in the definition
+ set_octahedron_points(inst);
+ set_obj_face(1,"1,2,5","b4fefe");
+ set_obj_face(2,"2,3,5","45d040");
+ set_obj_face(3,"3,4,5","4569d4");
+ set_obj_face(4,"4,1,5","b49bc0");
+ set_obj_face(5,"6,1,4","45a1d4");
+ set_obj_face(6,"6,2,1","b4c8fe");
+ set_obj_face(7,"6,3,2","b49b49");
+ set_obj_face(8,"6,4,3","112da1");
+enddef;
+
+
+%%\newpage\title{DODECAHEDRON (12 faces)\kern1cm\epsfbox{dodeca.ps}}
+
+%% Two adjacent faces (pentagons) of the dodecahedron are as follows:
+%%\figure{vect-fig.13}
+%% The angle between the faces is the angle between $\overrightarrow{OA}$
+%% and $\overrightarrow{OB}$.
+%% $OA=OB=d_3$ as defined in \figure{vect-fig.14}
+%% $d_1=2r\sin(\pi/5)$, $d_2=2r\sin(2\pi/5)$
+%% and $d_3=d_1\sin(2\pi/5)=2r\sin(\pi/5)\sin(2\pi/5)$
+%%
+%% The angle $\alpha$ is defined by the following conditions:
+%% \figure{vect-fig.15}
+%% Thus $d_2^2=d_3^2+d_3^2-2d_3^2\cos\alpha$, which leads to
+%% $\cos\alpha=1-{\displaystyle d_2^2\over\displaystyle2d_3^2}
+%% =-{\displaystyle\cos(2\pi/5)\over\displaystyle2\sin^2(\pi/5)}$
+%%\newpage
+
+def set_dodecahedron_points(expr inst)=
+ new_points(fc)(10);% face centers
+ set_point_(fc1,0,0,0);set_point(1)(1,0,0);
+ set_point(2)(cosd(72),sind(72),0);% 72=360/5
+ rotate_in_plane(3,fc1,1,2);
+ rotate_in_plane(4,fc1,2,3);
+ rotate_in_plane(5,fc1,3,4);
+ cosan=-cosd(72)/(2*sind(36)*sind(36));sinan=sqrt(1-cosan**2);
+ an=180-angle((cosan,sinan));
+ new_abs_face_point(fc2,fc1,1,2,an);
+ new_abs_face_point(fc3,fc1,2,3,an);
+ new_abs_face_point(fc4,fc1,3,4,an);
+ new_abs_face_point(fc5,fc1,4,5,an);
+ new_abs_face_point(fc6,fc1,5,1,an);
+ rotate_in_plane(6,fc2,2,1);
+ rotate_in_plane(7,fc2,1,6);
+ rotate_in_plane(8,fc2,6,7);
+ rotate_in_plane(9,fc3,2,8);
+ rotate_in_plane(10,fc3,8,9);
+ rotate_in_plane(11,fc4,3,10);
+ rotate_in_plane(12,fc4,10,11);
+ rotate_in_plane(13,fc5,4,12);
+ rotate_in_plane(14,fc5,12,13);
+ rotate_in_plane(15,fc6,5,14);
+ new_abs_face_point(fc7,fc2,6,7,an);
+ new_abs_face_point(fc8,fc3,8,9,an);
+ new_abs_face_point(fc9,fc4,10,11,an);
+ new_abs_face_point(fc10,fc5,12,13,an);
+ rotate_in_plane(16,fc7,6,15);
+ rotate_in_plane(17,fc7,15,16);
+ rotate_in_plane(18,fc8,7,17);
+ rotate_in_plane(19,fc9,9,18);
+ rotate_in_plane(20,fc10,11,19);
+ normalize_obj(inst)(2,20);% opposite vertices
+ set_point(21)(0,0,0);% center of dodecahedron
+ free_points(fc)(10);
+enddef;
+%%\newpage
+vardef def_dodecahedron(expr inst)=
+ save cosan,sinan,an;
+ new_obj_points(inst,21);% 21 points in the definition
+ new_obj_faces(inst,12);% 12 faces in the definition
+ set_dodecahedron_points(inst);
+ set_obj_face(1,"5,4,3,2,1","ff0fa1");
+ set_obj_face(2,"8,7,6,1,2","b40000");
+ set_obj_face(3,"10,9,8,2,3","b49b49");
+ set_obj_face(4,"12,11,10,3,4","b49bc0");
+ set_obj_face(5,"14,13,12,4,5","b4c8fe");
+ set_obj_face(6,"6,15,14,5,1","b4fefe");
+ set_obj_face(7,"6,7,17,16,15","b4fe40");
+ set_obj_face(8,"8,9,18,17,7","45d040");
+ set_obj_face(9,"10,11,19,18,9","45a114");
+ set_obj_face(10,"11,12,13,20,19","45a1d4");
+ set_obj_face(11,"14,15,16,20,13","4569d4");
+ set_obj_face(12,"16,17,18,19,20","112da1");
+enddef;
+
+%%\newpage\title{ICOSAHEDRON (20 faces)\kern1cm\epsfbox{icosa.ps}}
+
+%%Two faces of an icosahedron are linked in the following way:
+
+%% \vc{\epsfbox{vect-fig.4}} where $h=c\sqrt3/2$
+
+%% $d$ is also a diagonal in a pentagon: \vc{\epsfbox{vect-fig.5}}
+
+%% $c=2r\sin(\pi/5)$,
+%% $d=2r\sin(2\pi/5)=2c\cos(\pi/5)$
+
+%% The angle between two faces is computed as follows:
+%% \vc{\epsfbox{vect-fig.6}}
+
+%% We have ${(2\cos(\pi/5))}^2=2{(\sqrt3/2)}^2-2{(\sqrt3/2)}^2\cos\alpha$
+
+%% Thus $cos\alpha=1-{8\over3}\cos^2(\pi/5)$
+%%\newpage
+
+def set_icosahedron_points(expr inst)=
+ set_point(1)(0,0,0);set_point(2)(1,0,0);
+ set_point(3)(cosd(60),sind(60),0);% 60=360/6
+ cosan=1-8/3*cosd(36)*cosd(36);sinan=sqrt(1-cosan**2);
+ an=180-angle((cosan,sinan));
+ new_face_point(4,1,2,3,an);
+ new_face_point(5,2,3,1,an);
+ new_face_point(6,3,1,2,an);
+ new_face_point(7,2,4,3,an);
+ new_face_point(8,3,5,1,an);
+ new_face_point(9,1,6,2,an);
+ new_face_point(10,3,4,7,an);
+ new_face_point(11,3,7,5,an);
+ new_face_point(12,1,8,6,an);
+ normalize_obj(inst)(1,10);% opposite vertices
+ set_point(13)(0,0,0);% center of icosahedron
+enddef;
+
+vardef def_icosahedron(expr inst)=
+ save cosan,sinan,an;
+ new_obj_points(inst,13);% 12 vertices and the center
+ new_obj_faces(inst,20);% 20 faces in the definition
+ set_icosahedron_points(inst);
+ set_obj_face(1,"3,2,1","b40000");
+ set_obj_face(2,"2,3,4","ff0fa1");
+ set_obj_face(3,"3,7,4","b49b49");
+ set_obj_face(4,"3,5,7","b49bc0");
+ set_obj_face(5,"3,1,5","b4c8fe");
+ set_obj_face(6,"1,8,5","b4fefe");
+ set_obj_face(7,"1,6,8","b4fe40");
+ set_obj_face(8,"1,2,6","45d040");
+ set_obj_face(9,"2,9,6","45a114");
+ set_obj_face(10,"2,4,9","45a1d4");
+ set_obj_face(11,"9,4,10","4569d4");
+ set_obj_face(12,"4,7,10","112da1");
+ set_obj_face(13,"7,5,11","b4fefe");
+ set_obj_face(14,"5,8,11","b49bc0");
+ set_obj_face(15,"8,6,12","45a114");
+ set_obj_face(16,"6,9,12","b49b49");
+ set_obj_face(17,"8,12,11","b40000");
+ set_obj_face(18,"7,11,10","45a1d4");
+ set_obj_face(19,"12,10,11","b4c8fe");
+ set_obj_face(20,"9,10,12","ff0fa1");
+enddef;
+
+%%\newpage\title{General draw functions}
+
+def draw_polyhedron(expr name)=
+ draw_faces(name);
+enddef;
+
+let draw_tetrahedron=draw_polyhedron;
+let draw_cube=draw_polyhedron;
+let draw_octahedron=draw_polyhedron;
+let draw_dodecahedron=draw_polyhedron;
+let draw_icosahedron=draw_polyhedron;
+
+let new_poly=assign_obj;
+
+endinput
+
diff --git a/Master/texmf-dist/metapost/mp3d/3dutil.mp b/Master/texmf-dist/metapost/mp3d/3dutil.mp
new file mode 100644
index 00000000000..e83f5f18c05
--- /dev/null
+++ b/Master/texmf-dist/metapost/mp3d/3dutil.mp
@@ -0,0 +1,37 @@
+%%\input epsf
+%%\def\newpage{\vfill\eject}
+%%\advance\vsize1in
+%%\let\ora\overrightarrow
+%%\def\title#1{\hrule\vskip1mm#1\par\vskip1mm\hrule\vskip5mm}
+%%\def\figure#1{\par\centerline{\epsfbox{#1}}}
+%%\title{{\bf 3DUTIL.MP: 3D UTILITIES IN METAPOST}}
+
+%% version 1.34, 17 August 2003
+%% {\bf Denis Roegel} ({\tt roegel@loria.fr})
+
+% This package provides some basic functions used in certain
+% files of the 3d package.
+
+if known three_d_util_version:
+ expandafter endinput % avoids loading this package twice
+fi;
+
+message "*** 3dutil, v1.34 (c) D. Roegel 17 August 2003 ***";
+numeric three_d_util_version; three_d_util_version:=1.34;
+
+% maybe this exists?
+def sign(expr n)=
+ (if n>=0: 1 else: -1 fi)
+enddef;
+
+%
+def absmin(expr v,w)=
+ (if abs(v)<=abs(w):abs(v) else: abs(w) fi)
+enddef;
+
+%
+def absmax(expr v,w)=
+ (if abs(v)>abs(w):abs(v) else: abs(w) fi)
+enddef;
+
+endinput
diff --git a/Master/texmf-dist/metapost/mp3d/animpoly.mp b/Master/texmf-dist/metapost/mp3d/animpoly.mp
new file mode 100644
index 00000000000..99300950053
--- /dev/null
+++ b/Master/texmf-dist/metapost/mp3d/animpoly.mp
@@ -0,0 +1,35 @@
+% Animation of polyhedra
+% Example of use of 3d and poly packages
+% METAPOST
+% Denis Roegel, 17 August 2003
+
+% package 3d
+input 3d
+
+% polyhedra definitions
+input 3dpoly
+
+% animations
+input 3danim
+
+% set scale
+drawing_scale:=8cm;
+
+%
+filled_faces:=true; % default value
+%filled_faces:=false;
+
+%
+%show_animation_parameters:=true; % default value is |false|
+
+% animations
+%animate_object("tetrahedron",1,100,100);
+%animate_object("cube",1,100,100);
+%animate_object("octahedron",1,100,100);
+%animate_object("dodecahedron",1,100,100);
+%animate_object("icosahedron",1,100,100);
+%animate_two_objects("dodecahedron","icosahedron",1,100,100);
+animate_three_objects("dodecahedron","icosahedron","octahedron",1,100,100);
+%animate_two_identical_objects("dodecahedron",1,100,100);
+
+end
diff --git a/Master/texmf-dist/metapost/mp3d/create_animation.sh b/Master/texmf-dist/metapost/mp3d/create_animation.sh
new file mode 100644
index 00000000000..e99ca7383ca
--- /dev/null
+++ b/Master/texmf-dist/metapost/mp3d/create_animation.sh
@@ -0,0 +1,15 @@
+#! /bin/sh
+
+/bin/rm -f animpoly.log
+for i in `ls animpoly.*| grep 'animpoly.[0-9]'`;do
+echo $i
+echo '=============='
+awk < $i '{print} /^%%Page: /{print "142 123 translate\n"}' > $i.ps
+gs -sDEVICE=ppmraw -sPAPERSIZE=a4 -dNOPAUSE -r36 -sOutputFile=$i.ppm -q -- $i.ps
+/bin/rm -f $i.ps
+ppmquant 32 $i.ppm | pnmcut 0 114 141 307 | ppmtogif > `expr $i.ppm : '\(.*\)ppm'`gif
+/bin/rm -f $i.ppm
+done
+/bin/rm -f animpoly.gif
+gifmerge -10 -l1000 animpoly.*.gif > animpoly.gif
+/bin/rm -f animpoly.*.gif