From ac3c55a3216b5988f0e48ba9414ddb059f19a699 Mon Sep 17 00:00:00 2001 From: Karl Berry Date: Wed, 11 Jan 2006 22:58:36 +0000 Subject: trunk/Master/texmf-dist/metapost git-svn-id: svn://tug.org/texlive/trunk@104 c570f23f-e606-0410-a88d-b1316a301751 --- Master/texmf-dist/metapost/mp3d/3d.mp | 1259 ++++++++++++++++++++ Master/texmf-dist/metapost/mp3d/3danim.mp | 269 +++++ Master/texmf-dist/metapost/mp3d/3dgeom.mp | 931 +++++++++++++++ Master/texmf-dist/metapost/mp3d/3dpoly.mp | 294 +++++ Master/texmf-dist/metapost/mp3d/3dutil.mp | 37 + Master/texmf-dist/metapost/mp3d/animpoly.mp | 35 + .../texmf-dist/metapost/mp3d/create_animation.sh | 15 + 7 files changed, 2840 insertions(+) create mode 100644 Master/texmf-dist/metapost/mp3d/3d.mp create mode 100644 Master/texmf-dist/metapost/mp3d/3danim.mp create mode 100644 Master/texmf-dist/metapost/mp3d/3dgeom.mp create mode 100644 Master/texmf-dist/metapost/mp3d/3dpoly.mp create mode 100644 Master/texmf-dist/metapost/mp3d/3dutil.mp create mode 100644 Master/texmf-dist/metapost/mp3d/animpoly.mp create mode 100644 Master/texmf-dist/metapost/mp3d/create_animation.sh (limited to 'Master/texmf-dist/metapost/mp3d') 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 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))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 -- cgit v1.2.3