diff options
Diffstat (limited to 'Master/texmf-dist/metapost/mp3d/3dgeom.mp')
-rw-r--r-- | Master/texmf-dist/metapost/mp3d/3dgeom.mp | 931 |
1 files changed, 931 insertions, 0 deletions
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 |