%%\input epsf
%%\def\newpage{\vfill\eject}
%%\advance\vsize1in
%%\let\ora\overrightarrow
%%\def\title#1{\hrule\vskip1mm#1\par\vskip1mm\hrule\vskip5mm}
%%\def\figure#1{\par\centerline{\epsfbox{#1}}}
%%\title{{\bf 3D.MP: 3-DIMENSIONAL REPRESENTATIONS IN METAPOST}}

%% version 1.34, 17 August 2003
%% {\bf Denis Roegel} ({\tt roegel@loria.fr}) 

%% This package provides definitions enabling the manipulation
%% and animation of 3-dimensional objects.
%% Such objects can be included in a \TeX{} file or used on web pages
%% for instance. See the documentation enclosed in the distribution for
%% more details.

%% Thanks to John Hobby and Ulrik Vieth for helpful hints.

%% PROJECTS FOR THE FUTURE:

%%   $-$ take light sources into account and show shadows and darker faces

%%   $-$ handle overlapping of objects ({\it obj\_name\/} can be used when
%%     going through all faces)

if known three_d_version: 
  expandafter endinput % avoids loading this package twice
fi;

message "*** 3d,          v1.34 (c) D. Roegel, 17 August 2003 ***";
numeric three_d_version;
three_d_version=1.34;

% This package needs |3dgeom| in a few places. |3dgeom| also loads |3d|
% but that's not a problem.
%
input 3dgeom; 

%%\newpage
%%\title{Vector operations}

% components of vector |i|
def xval(expr i)=vec[i]x enddef;
def yval(expr i)=vec[i]y enddef;
def zval(expr i)=vec[i]z enddef;

% vector (or point) equality (absolute version)
def vec_eq_(expr i,j)=
  ((xval(i)=xval(j)) and (yval(i)=yval(j)) and (zval(i)=zval(j)))
enddef;

% vector (or point) equality (local version)
def vec_eq(expr i,j)=vec_eq_(pnt(i),pnt(j)) enddef;

% vector inequality (absolute version)
def vec_neq_(expr i,j)=(not vec_eq_(i,j)) enddef;

% vector inequality (local version)
def vec_neq(expr i,j)=(not vec_eq(i,j)) enddef;

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

% definition of vector |i| by its coordinates (local version)
def vec_def(expr i,xi,yi,zi)= vec_def_(pnt(i),xi,yi,zi) enddef;

% a point is stored as a vector (absolute version)
let set_point_ = vec_def_;

% a point is stored as a vector (local version)
let set_point = vec_def;

def set_point_vec_(expr i,v)=
  set_point_(i,xval(v),yval(v),zval(v))
enddef;

def set_point_vec(expr i,v)=set_point_vec_(pnt(i),v) enddef;

let vec_def_vec_=set_point_vec_;
let vec_def_vec=set_point_vec;

% vector sum: |vec[k]| $\leftarrow$ |vec[i]|$+$|vec[j]| (absolute version)
def vec_sum_(expr k,i,j)=
  vec[k]x:=vec[i]x+vec[j]x;
  vec[k]y:=vec[i]y+vec[j]y;
  vec[k]z:=vec[i]z+vec[j]z;
enddef;

% vector sum: |vec[k]| $\leftarrow$ |vec[i]|$+$|vec[j]| (local version)
def vec_sum(expr k,i,j)=vec_sum_(pnt(k),pnt(i),pnt(j)) enddef;

% vector translation: |vec[i]| $\leftarrow$ |vec[i]|$+$|vec[v]|
def vec_translate_(expr i,v)=vec_sum_(i,i,v) enddef;

% Here, the second parameter is absolute, because this is probably
% the most common case.
def vec_translate(expr i,v)=vec_translate_(pnt(i),v) enddef;

% vector difference: |vec[k]| $\leftarrow$ |vec[i]|$-$|vec[j]|
def vec_diff_(expr k,i,j)=
  vec[k]x:=vec[i]x-vec[j]x;
  vec[k]y:=vec[i]y-vec[j]y;
  vec[k]z:=vec[i]z-vec[j]z;
enddef;

def vec_diff(expr k,i,j)=vec_diff_(pnt(k),pnt(i),pnt(j)) enddef;

% dot product of |vec[i]| and |vec[j]|
vardef vec_dprod_(expr i,j)=
  (vec[i]x*vec[j]x+vec[i]y*vec[j]y+vec[i]z*vec[j]z)
enddef;

vardef vec_dprod(expr i,j)=vec_dprod_(pnt(i),pnt(j)) enddef;

% modulus of |vec[i]|, absolute version
% In the computation, we try to avoid overflows or underflows;
% we perform a scaling in order to avoid losing too much
% information in certain cases
vardef vec_mod_(expr i)=
  save prod,m_;
  hide(
    new_vec(v_a);
    m_=max(abs(xval(i)),abs(yval(i)),abs(zval(i)));
    if m_>0:vec_mult_(v_a,i,1/m_);else:vec_def_vec_(v_a,vec_null);fi;
    prod=m_*sqrt(vec_dprod_(v_a,v_a));
    free_vec(v_a);
  )
  prod      
enddef;

% modulus of |vec[i]|, local version
% If the return value must be compared to 0,
% use |vec_eq| with |vec_null| instead.
vardef vec_mod(expr i)= vec_mod_(pnt(i)) enddef;

% unit vector |vec[i]| corresponding to vector |vec[j]|
% only non-null vectors are changed
def vec_unit_(expr i,j)=
  if vec_mod_(j)>0: vec_mult_(i,j,1/vec_mod_(j));
  else:vec_def_vec_(i,j);
  fi;
enddef;

def vec_unit(expr i,j)=vec_unit_(pnt(i),pnt(j)) enddef;

% vector product: |vec[k]| $\leftarrow$ |vec[i]| $\land$ |vec[j]|
def vec_prod_(expr k,i,j)=
  vec[k]x:=vec[i]y*vec[j]z-vec[i]z*vec[j]y;
  vec[k]y:=vec[i]z*vec[j]x-vec[i]x*vec[j]z;
  vec[k]z:=vec[i]x*vec[j]y-vec[i]y*vec[j]x;
enddef;

def vec_prod(expr k,i,j)=vec_prod_(pnt(k),pnt(i),pnt(j)) enddef;

% scalar multiplication: |vec[j]| $\leftarrow$ |vec[i]*v| (absolute version)
def vec_mult_(expr j,i,v)=
  vec[j]x:=v*vec[i]x;vec[j]y:=v*vec[i]y;vec[j]z:=v*vec[i]z;
enddef;

% scalar multiplication: |vec[j]| $\leftarrow$ |vec[i]*v| (local version)
def vec_mult(expr j,i,v)=vec_mult_(pnt(j),pnt(i),v) enddef;

% middle of two points (absolute version)
def mid_point_(expr k,i,j)= vec_sum_(k,i,j);vec_mult_(k,k,.5); enddef;

% middle of two points (local version)
def mid_point(expr k,i,j)= mid_point_(pnt(k),pnt(i),pnt(j)); enddef;

%%\newpage
%%\title{Vector rotation}
% Rotation of |vec[v]| around |vec[axis]| by an angle |alpha|

%% The vector $\vec{v}$ is first projected on the axis
%% giving vectors $\vec{a}$ and $\vec{h}$:
%%\figure{vect-fig.9}
%% If we set 
%% $\vec{b}={\ora{axis}\over \left\Vert\vcenter{\ora{axis}}\right\Vert}$,
%% the rotated vector $\vec{v'}$ is equal to $\vec{h}+\vec{f}$
%% where $\vec{f}=\cos\alpha \cdot \vec{a} + \sin\alpha\cdot \vec{c}$.
%% and $\vec{h}=(\vec{v}\cdot\vec{b})\vec{b}$ 
%%\figure{vect-fig.10}

% The rotation is independent of |vec[axis]|'s module.
% |v| = old and new vector
% |axis| = rotation axis
% |alpha| = rotation angle
%
vardef vec_rotate_(expr v,axis,alpha)=
  new_vec(v_a);new_vec(v_b);new_vec(v_c);
  new_vec(v_d);new_vec(v_e);new_vec(v_f);
  new_vec(v_g);new_vec(v_h);
  vec_mult_(v_b,axis,1/vec_mod_(axis));
  vec_mult_(v_h,v_b,vec_dprod_(v_b,v)); % projection of |v| on |axis|
  vec_diff_(v_a,v,v_h);
  vec_prod_(v_c,v_b,v_a);
  vec_mult_(v_d,v_a,cosd(alpha));
  vec_mult_(v_e,v_c,sind(alpha));
  vec_sum_(v_f,v_d,v_e);
  vec_sum_(v,v_f,v_h);
  free_vec(v_h);free_vec(v_g);
  free_vec(v_f);free_vec(v_e);free_vec(v_d);
  free_vec(v_c);free_vec(v_b);free_vec(v_a);
enddef;

% The second parameter is left absolute because this is probably the most
% common case.
vardef vec_rotate(expr v,axis,alpha)=vec_rotate_(pnt(v),axis,alpha) enddef;

%%\newpage
%%\title{Operations on objects}
% |iname| is the handler for an instance of an object of class |name|
% |iname| must be a letter string
% |vardef| is not used because at some point we give other names
% to |assign_obj| with |let| and this cannot be done with |vardef|.
% (see MFbook for details)
def assign_obj(expr iname,name)=
  begingroup 
  save tmpdef;
  string tmpdef; % we need to add double quotes (char 34)
  tmpdef="def " & iname & "_class=" & ditto & name & ditto & " enddef";
  scantokens tmpdef;
  def_obj(iname);
  endgroup
enddef;

% |name| is the the name of an object instance
% It must be made only of letters (or underscores), but no digits.
def def_obj(expr name)=
  scantokens begingroup 
    save tmpdef;string tmpdef;
    tmpdef="def_" & obj_class_(name) & "(" & ditto & name & ditto & ")";
    tmpdef
  endgroup
enddef;

% This macro puts an object back where it was right at the beginning,
% or rather, where the |set| definition puts it (which may be different
% than the initial position, in case it depends on parameters).
% |iname| is the name of an object instance.
vardef reset_obj(expr iname)=
  save tmpdef;
  string tmpdef;
  define_current_point_offset_(iname);
  tmpdef="set_" & obj_class_(iname) & "_points";
  scantokens tmpdef(iname);
enddef;

% Put an object at position given by |pos| (a vector) and
% with orientations given by angles |psi|, |theta|, |phi|.
% The object is scaled by |scale|.
% |iname| is the name of an object instance.
% If the shape of the object has been changed since it was
% created, these changes are lost.
vardef put_obj(expr iname,pos,scale,psi,theta,phi)=
  reset_obj(iname);scale_obj(iname,scale);
  new_vec(v_x);new_vec(v_y);new_vec(v_z);
  vec_def_vec_(v_x,vec_I);
  vec_def_vec_(v_y,vec_J);
  vec_def_vec_(v_z,vec_K);
  rotate_obj_abs_pv(iname,point_null,v_z,psi);
  vec_rotate_(v_x,v_z,psi);vec_rotate_(v_y,v_z,psi);
  rotate_obj_abs_pv(iname,point_null,v_y,theta);
  vec_rotate_(v_x,v_y,theta);vec_rotate_(v_z,v_y,theta);
  rotate_obj_abs_pv(iname,point_null,v_x,phi);
  vec_rotate_(v_y,v_x,phi);vec_rotate_(v_z,v_x,phi);
  free_vec(v_z);free_vec(v_y);free_vec(v_x);    
  translate_obj(iname,pos);
enddef;

%%\newpage
%%\title{Rotation, translation and scaling of objects}
% Rotation of an object instance |name| around an axis 
% going through a point |p| (local to the object)
% and directed by vector |vec[v]|. The angle of rotation is |a|.
vardef rotate_obj_pv(expr name,p,v,a)=
  define_current_point_offset_(name);
  rotate_obj_abs_pv(name,pnt(p),v,a);
enddef;

vardef rotate_obj_abs_pv(expr name,p,v,a)=
  define_current_point_offset_(name);
  new_vec(v_a);
  for i:=1 upto obj_points_(name):
    vec_diff_(v_a,pnt(i),p);
    vec_rotate_(v_a,v,a);
    vec_sum_(pnt(i),v_a,p);
  endfor;
  free_vec(v_a);
enddef;

% Rotation of an object instance |name| around an axis 
% going through a point |p| (local to the object)
% and directed by vector $\ora{pq}$. The angle of rotation is |a|.
vardef rotate_obj_pp(expr name,p,q,a)=
  define_current_point_offset_(name);
  new_vec(v_a);new_vec(axis);
  vec_diff_(axis,pnt(q),pnt(p));
  for i:=1 upto obj_points_(name):
    vec_diff_(v_a,pnt(i),pnt(p));
    vec_rotate_(v_a,axis,a);
    vec_sum_(pnt(i),v_a,pnt(p));
  endfor;
  free_vec(axis);free_vec(v_a);
enddef;

% Translation of an object instance |name| by a vector |vec[v]|.
vardef translate_obj(expr name,v)=
  define_current_point_offset_(name);
  for i:=1 upto obj_points_(name):
    vec_sum_(pnt(i),pnt(i),v);
  endfor;
enddef;

% Scalar multiplication of an object instance |name| by a scalar |v|.
vardef scale_obj(expr name,v)=
  define_current_point_offset_(name);
  for i:=1 upto obj_points_(name):
    vec_mult(i,i,v);
  endfor;
enddef;


%%\newpage
%%\title{Functions to build new points in space}
% Rotation in a plane: this is useful to define a regular polygon.
% |k| is a new point obtained from point |j| by rotation around |o|
% by a angle $\alpha$ equal to the angle from |i| to |j|.
%%\figure{vect-fig.11}
vardef rotate_in_plane_(expr k,o,i,j)=
  save cosalpha,sinalpha,alpha;
  new_vec(v_a);new_vec(v_b);new_vec(v_c);
  vec_diff_(v_a,i,o);vec_diff_(v_b,j,o);vec_prod_(v_c,v_a,v_b);
  cosalpha=vec_dprod_(v_a,v_b)/vec_mod_(v_a)/vec_mod_(v_b);
  sinalpha=sqrt(1-cosalpha**2);
  alpha=angle((cosalpha,sinalpha));
  vec_rotate_(v_b,v_c,alpha);
  vec_sum_(k,o,v_b);
  free_vec(v_c);free_vec(v_b);free_vec(v_a);
enddef;

vardef rotate_in_plane(expr k,o,i,j)=
  rotate_in_plane_(pnt(k),o,pnt(i),pnt(j)) 
enddef;

% Build a point on a adjacent face.
%% The middle $m$ of points $i$ and $j$ is such that
%% $\widehat{(\ora{om},\ora{mc})}=\alpha$ 
%% This is useful to define regular polyhedra
%%\figure{vect-fig.7}
vardef new_face_point_(expr c,o,i,j,alpha)=
  new_vec(v_a);new_vec(v_b);new_vec(v_c);new_vec(v_d);new_vec(v_e);
  vec_diff_(v_a,i,o);vec_diff_(v_b,j,o);
  vec_sum_(v_c,v_a,v_b);
  vec_mult_(v_d,v_c,.5);
  vec_diff_(v_e,i,j);
  vec_sum_(c,v_d,o);
  vec_rotate_(v_d,v_e,alpha);
  vec_sum_(c,v_d,c);
  free_vec(v_e);free_vec(v_d);free_vec(v_c);free_vec(v_b);free_vec(v_a);
enddef;

vardef new_face_point(expr c,o,i,j,alpha)=
  new_face_point_(pnt(c),pnt(o),pnt(i),pnt(j),alpha)
enddef;

vardef new_abs_face_point(expr c,o,i,j,alpha)=
  new_face_point_(c,o,pnt(i),pnt(j),alpha)
enddef;

%%\newpage
%%\title{Computation of the projection of a point on the ``screen''}
% |p| is the projection of |m|
% |m| = point in space (3 coordinates)
% |p| = point of the intersection plane 
%%\figure{vect-fig.8}
vardef project_point(expr p,m)=
  save tmpalpha;
  new_vec(v_a);new_vec(v_b);
  if projection_type=2: % oblique
    if point_in_plane_p_pl_(m)(projection_plane):
      % |m| is on the projection plane
      vec_diff_(v_a,m,ObliqueCenter_);
      y[p]:=drawing_scale*vec_dprod_(v_a,ProjJ_);
      x[p]:=drawing_scale*vec_dprod_(v_a,ProjK_);
    else: % |m| is not on the projection plane
      new_line_(l)(m,ObliqueCenter_);
      vec_diff_(l2,l2,Obs);
      vec_sum_(l2,l2,m);
      % (the direction does not depend on Obs)
      if def_inter_p_l_pl_(v_a)(l)(projection_plane):
        vec_diff_(v_a,v_a,ObliqueCenter_);
        y[p]:=drawing_scale*vec_dprod_(v_a,ProjJ_);
        x[p]:=drawing_scale*vec_dprod_(v_a,ProjK_);
      else: message "Point " & decimal m & " cannot be projected";
        x[p]:=too_big_;y[p]=too_big_;
      fi;
      free_line(l);
    fi;
  else:
    vec_diff_(v_b,m,Obs); % vector |Obs|-|m|
      % |vec[v_a]| is |vec[v_b]| expressed in (|ObsI_|,|ObsJ_|,|ObsK_|)
      % coordinates.
    vec[v_a]x:=vec[IObsI_]x*vec[v_b]x
    +vec[IObsJ_]x*vec[v_b]y+vec[IObsK_]x*vec[v_b]z;
    vec[v_a]y:=vec[IObsI_]y*vec[v_b]x
    +vec[IObsJ_]y*vec[v_b]y+vec[IObsK_]y*vec[v_b]z;
    vec[v_a]z:=vec[IObsI_]z*vec[v_b]x
    +vec[IObsJ_]z*vec[v_b]y+vec[IObsK_]z*vec[v_b]z;
    if vec[v_a]x<Obs_dist: % then, point |m| is too close
      message "Point " & decimal m & " too close -> not drawn";
      x[p]:=too_big_;y[p]=too_big_;
    else:
      if (angle(vec[v_a]x,vec[v_a]z)>h_field/2)
        or (angle(vec[v_a]x,vec[v_a]y)>v_field/2):
        message "Point " & decimal m & " out of screen -> not drawn";
        x[p]:=too_big_;y[p]=too_big_;
      else:
        if projection_type=0: % central perspective
	  tmpalpha:=Obs_dist/vec[v_a]x;
        else:
	  tmpalpha:=1; % parallel
        fi;
        y[p]:=drawing_scale*tmpalpha*vec[v_a]y;
        x[p]:=drawing_scale*tmpalpha*vec[v_a]z;
      fi;
    fi;
  fi;
  free_vec(v_b);free_vec(v_a);
enddef;

% At some point, we may need to do an oblique projection
% of vectors |ObsK_| and |ObsI_| on a plane, and to normalize
% and orthogonalize the projections (with the projection of |ObsK_|
% keeping the same direction). This is done here,
% where we take two vectors, a direction (line) and
% a plane, and return two vectors. This function assumes
% there is an intersection between line |l| and plane |p|.
% We do not test it here.

vardef project_vectors(expr va,vb)(expr k,i)(text l)(text p)=
  save vc;new_vec(vc);
  if proj_v_v_l_pl_(va,k)(l)(p): % |va| is the projection of vector |k|
  else: message "THIS SHOULD NOT HAPPEN";
  fi;
  if proj_v_v_l_pl_(vb,i)(l)(p): % |vb| is the projection of vector |i|
  else: message "THIS SHOULD NOT HAPPEN";
  fi;
  % now, we orthonormalize these vectors:
  vec_prod_(vc,va,vb);
  vec_unit_(va,va);vec_unit_(vc,vc);vec_prod_(vb,vc,va);
  free_vec(vc);
enddef;

% Object projection
% This is a mere iteration on |project_point|
def project_obj(expr name)=
  define_current_point_offset_(name);
  for i:=1 upto obj_points_(name):
    project_point(ipnt_(i),pnt(i));endfor;
enddef;

% Projection screen
vardef show_projection_screen=
  save dx,dy;
  dx=Obs_dist*sind(h_field/2)/cosd(h_field/2);
  dy=Obs_dist*sind(v_field/2)/cosd(v_field/2);
  new_vec(pa);new_vec(pb);new_vec(pc);new_vec(pd);new_vec(op);
  new_vec(w);new_vec(h);
  vec_mult_(op,ObsI_,Obs_dist);vec_sum_(op,op,Obs); % center of screen
  vec_mult_(w,ObsK_,dx);vec_mult_(h,ObsJ_,dy);
  vec_sum_(pa,op,w);vec_sum_(pa,pa,h); % upper right corner
  vec_mult_(w,w,-2);vec_mult_(h,h,-2);
  vec_sum_(pb,pa,w);vec_sum_(pc,pb,h);vec_sum_(pd,pa,h);
  message "Screen at corners:";
  show_point("urcorner: ",pa);
  show_point("ulcorner: ",pb);
  show_point("llcorner: ",pc);
  show_point("lrcorner: ",pd);
  show_point("Obs:",Obs);
  free_vec(h);free_vec(w);
  free_vec(op);free_vec(pd);free_vec(pc);free_vec(pb);free_vec(pa);
enddef;


%%\newpage
%%\title{Draw one face, hiding it if it is hidden}
% The order of the vertices determines what is the visible side
% of the face. The order must be clockwise when the face is seen.
% |drawhidden| is a boolean; if |true| only hidden faces are drawn; if |false|,
% only visible faces are drawn. Therefore, |draw_face| is called twice
% by |draw_faces|.
vardef draw_face(text vertices)(expr col,drawhidden)=
  save p,num,overflow,i,j,k,nv;
  path p;boolean overflow;
  overflow=false;
  forsuffixes $=vertices:
    if z[ipnt_($)]=(too_big_,too_big_):overflow:=true; fi;
    exitif overflow;
  endfor;
  if overflow: message "Face can not be drawn, due to overflow";
  else:
    p=forsuffixes $=vertices:z[ipnt_($)]--endfor cycle;
    % we do now search for three distinct and non-aligned suffixes:
    % usually, the first three suffixes do
    new_vec(normal_vec);new_vec(v_a);new_vec(v_b);new_vec(v_c);
    % first, we copy all the indexes in an array, so that
    % it is easier to go through them
    i=1; % num0 is not used
    forsuffixes $=vertices:num[i]=$;i:=i+1;endfor;
    nv=i-1;
    for $:=1 upto nv:
      for $$:=$+1 upto nv:
	for $$$:=$$+1 upto nv:
	  vec_diff_(v_a,pnt(num[$$]),pnt(num[$]));
          vec_diff_(v_b,pnt(num[$$$]),pnt(num[$$]));
          vec_prod_(normal_vec,v_a,v_b);
          exitif vec_neq_(normal_vec,vec_null);
	      % |vec_mod_| must not be used for such a test
	endfor;
	exitif vec_neq_(normal_vec,vec_null);
      endfor;
      exitif vec_neq_(normal_vec,vec_null);
    endfor;
    if projection_type=0: % perspective
      vec_diff_(v_c,pnt(num1),Obs);
    else: % parallel
      vec_def_vec_(v_c,ObsI_); 
    fi;
    if filled_faces:
      if vec_dprod_(normal_vec,v_c)<0:
        fill p withcolor col;drawcontour(p,contour_width,contour_color)();
      else: % |draw p dashed evenly;| if this is done, you must ensure
              % that hidden faces are (re)drawn at the end    
      fi;
    else:
      if vec_dprod_(normal_vec,v_c)<0:%visible
        if not drawhidden:drawcontour(p,contour_width,contour_color)();fi;
      else: % hidden
        if drawhidden:
	  drawcontour(p,contour_width,contour_color)(dashed evenly);
        fi;
      fi;
    fi;
    free_vec(v_c);free_vec(v_b);free_vec(v_a);free_vec(normal_vec);
  fi;
enddef;

% |p| is the path to draw (a face contour), |thickness| is the pen width
% |col| is the color and |type| is a line modifier.
def drawcontour(expr p,thickness,col)(text type)=
  if draw_contours and (thickness>0):
    pickup pencircle scaled thickness;
    draw p withcolor background; % avoid strange overlapping dashes
    draw p type withcolor col;
    pickup pencircle scaled .4pt;
  fi;
enddef;

%%\newpage
% Variables for face handling. First, we have an array for lists of vertices
% corresponding to faces. 
string face_points_[];% analogous to |vec| arrays

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

% |name| is the name of an object instance
vardef draw_faces(expr name)=
  save tmpdef;string tmpdef;
  define_current_face_offset_(name);
    % first the hidden faces (dashes must be drawn first):
  for i:=1 upto obj_faces_(name):
    tmpdef:="draw_face(" & face_points_[face(i)] 
      & ")(hexcolor(" & ditto & face_color_[face(i)] & ditto 
      & "),true)";scantokens tmpdef;
  endfor;
    % then, the visible faces:
  for i:=1 upto obj_faces_(name):
    tmpdef:="draw_face(" & face_points_[face(i)] 
      & ")(hexcolor(" & ditto & face_color_[face(i)] & ditto 
      & "),false)";scantokens tmpdef;
  endfor;
enddef;

% Draw point |n| of object instance |name|
vardef draw_point(expr name,n)=
  define_current_point_offset_(name);
  project_point(ipnt_(n),pnt(n));
  if z[ipnt_(n)] <> (too_big_,too_big_):
    pickup pencircle scaled 5pt;
    drawdot(z[ipnt_(n)]);
    pickup pencircle scaled .4pt;
  fi;
enddef;

vardef draw_axes(expr r,g,b)=
  project_point(1,vec_null);
  project_point(2,vec_I);
  project_point(3,vec_J);
  project_point(4,vec_K);
  if (z1<>(too_big_,too_big_)):
    if (z2<>(too_big_,too_big_)):
      drawarrow z1--z2 dashed evenly withcolor r;
    fi;
    if (z3<>(too_big_,too_big_)):
      drawarrow z1--z3 dashed evenly withcolor g;
    fi;
    if (z4<>(too_big_,too_big_)):
      drawarrow z1--z4 dashed evenly withcolor b;
    fi;
  fi;
enddef;

% Draw a polygonal line through the list of points
% This implementation does not work if you call
% |draw_lines(i,i+4)| because \MP{} adds parentheses around
% the value of |i|.
def draw_lines(text vertices)=
  begingroup % so that we can |let| |draw_lines|
  save j,num,np;
  % first, we copy all the indexes in an array, so that
  % it is easier to go through them
  j=1;
  for $=vertices:num[j]=$;j:=j+1;endfor;
  np=j-1;
  for j:=1 upto np-1:
    draw z[ipnt_(num[j])]--z[ipnt_(num[j+1])];
  endfor;
  endgroup
enddef;

let draw_line=draw_lines;

% Draw an arrow between points |i| and |j| of current object
% This is used from the |draw| definition of an object.
def draw_arrow(expr i,j)=
  drawarrow z[ipnt_(i)]--z[ipnt_(j)];
enddef;

% Draw a line between points |i| of object |obja| and |j| of |objb|
% This is used when outside an object (i.e., we can't presuppose
% any object offset)
vardef draw_line_inter(expr obja, i, objb, j)=
  project_point(1,pnt_obj(obja,i));
  project_point(2,pnt_obj(objb,j));
  draw z1--z2;
enddef;

% Draw an arrow between points |i| of object |obja| and |j| of |objb|
% This is used when outside an object (i.e., we can't presuppose
% any object offset)
vardef draw_arrow_inter(expr obja, i, objb, j)=
  project_point(1,pnt_obj(obja,i));
  project_point(2,pnt_obj(objb,j));
  draw z1--z2;
enddef;

%%\newpage
% Definition of a macro |obj_name| returning an object name 
% when given an absolute
% face number. This definition is built incrementally through a string, 
% everytime a new object is defined.
% |obj_name| is defined by |redefine_obj_name_|.

% Initial definition
string index_to_name_;
index_to_name_="def obj_name(expr i)=if i<1:";

% |name| is the name of an object instance
% |n| is the absolute index of its last face
def redefine_obj_name_(expr name,n)=
  index_to_name_:=index_to_name_ & "elseif i<=" & decimal n & ":" & ditto
      & name & ditto;
  scantokens begingroup index_to_name_ & "fi;enddef;" endgroup;
enddef;

% |i| is an absolute face number
% |vertices| is a string representing a list of vertices
% |rgbcolor| is a string representing a color in rgb hexadecimal
def set_face(expr i,vertices,rgbcolor)=
  face_points_[i]:=vertices;face_color_[i]:=rgbcolor;
enddef;

% |i| is a local face number
% |vertices| is a string representing a list of vertices
% |rgbcolor| is a string representing a color in rgb hexadecimal
def set_obj_face(expr i,vertices,rgbcolor)=set_face(face(i),vertices,rgbcolor)
enddef;

% |i| is a local face number of object |inst|
% |rgbcolor| is a string representing a color in rgb hexadecimal
def set_obj_face_color(expr inst,i,rgbcolor)=
  face_color_[face_obj(inst,i)]:=rgbcolor;
enddef;


%%\newpage
%%\title{Compute the vectors corresponding to the observer's viewpoint}
% (vectors |ObsI_|,|ObsJ_| and |ObsK_| in the |vec_I|,|vec_J|,
% |vec_K| reference; and vectors |IObsI_|,|IObsJ_| and |IObsK_| 
% which are |vec_I|,|vec_J|,|vec_K| 
% in the |ObsI_|,|ObsJ_|,|ObsK_| reference)
%%\figure{vect-fig.16}
%% (here, $\psi>0$, $\theta<0$ and $\phi>0$; moreover, 
%% $\vert\theta\vert \leq 90^\circ$)

def compute_reference(expr psi,theta,phi)=
   % |ObsI_| defines the direction of observation; 
   % |ObsJ_| and |ObsK_| the orientation
   % (but one of these two vectors is enough,
   % since |ObsK_| = |ObsI_| $\land$ |ObsJ_|)
   % The vectors are found by rotations of |vec_I|,|vec_J|,|vec_K|.
  vec_def_vec_(ObsI_,vec_I);vec_def_vec_(ObsJ_,vec_J);
  vec_def_vec_(ObsK_,vec_K);
  vec_rotate_(ObsI_,ObsK_,psi);
  vec_rotate_(ObsJ_,ObsK_,psi);% gives ($u$,$v$,$z$)
  vec_rotate_(ObsI_,ObsJ_,theta);
  vec_rotate_(ObsK_,ObsJ_,theta);% gives ($Obs_x$,$v$,$w$)
  vec_rotate_(ObsJ_,ObsI_,phi);
  vec_rotate_(ObsK_,ObsI_,phi);% gives ($Obs_x$,$Obs_y$,$Obs_z$)
   % The passage matrix $P$ from |vec_I|,|vec_J|,|vec_K| 
   % to |ObsI_|,|ObsJ_|,|ObsK_| is the matrix
   % composed of the vectors |ObsI_|,|ObsJ_| and |ObsK_| expressed 
   % in the base |vec_I|,|vec_J|,|vec_K|.
   % We have $X=P X'$ where $X$ are the coordinates of a point 
   % in |vec_I|,|vec_J|,|vec_K|
   % and $X'$ the coordinates of the same point in |ObsI_|,|ObsJ_|,|ObsK_|.
   % In order to get $P^{-1}$, it suffices to build vectors using
   % the previous rotations in the inverse order.
  vec_def_vec_(IObsI_,vec_I);vec_def_vec_(IObsJ_,vec_J);
  vec_def_vec_(IObsK_,vec_K);
  vec_rotate_(IObsK_,IObsI_,-phi);vec_rotate_(IObsJ_,IObsI_,-phi);
  vec_rotate_(IObsK_,IObsJ_,-theta);vec_rotate_(IObsI_,IObsJ_,-theta);
  vec_rotate_(IObsJ_,IObsK_,-psi);vec_rotate_(IObsI_,IObsK_,-psi);
enddef;

%%\newpage
%%\title{Point of view}
% This macro computes the three angles necessary for |compute_reference|
% |name| =  name of an instance of an object 
% |target| = target point (local to object |name|)
% |phi| = angle
vardef point_of_view_obj(expr name,target,phi)=
  define_current_point_offset_(name);% enables |pnt|
  point_of_view_abs(pnt(target),phi);
enddef;

% Compute absolute perspective. |target| is an absolute point number
% |phi| = angle
% This function also computes two vectors needed in case
% of an oblique projection.
vardef point_of_view_abs(expr target,phi)=
  save psi,theta;
  new_vec(v_a);
  vec_diff_(v_a,target,Obs);
  vec_mult_(v_a,v_a,1/vec_mod_(v_a));
  psi=angle((vec[v_a]x,vec[v_a]y));
  theta=-angle((vec[v_a]x++vec[v_a]y,vec[v_a]z));
  compute_reference(psi,theta,phi);
  if projection_type=2: % oblique
    % we start by checking that at a minimum the three points defining
    % the projection plane have different indexes; it doesn't mean
    % the plane if well defined, but if two values are identical,
    % the plane can't be well defined.
    if ((projection_plane1<>projection_plane2) and
      (projection_plane1<>projection_plane3) and
      (projection_plane2<>projection_plane3)):
      new_line_(l)(Obs,Obs);
      vec_sum_(l2,ObsI_,Obs);
      if def_inter_p_l_pl_(ObliqueCenter_)(l)(projection_plane):
        project_vectors(ProjK_,ProjJ_)(ObsK_,ObsJ_)(l)(projection_plane);
	% define the projection direction
	set_line_(projection_direction)(Obs,ObliqueCenter_);
      else:
	message "Anomalous oblique projection:";
	message "  the observer is watching parallely to the plane";
      fi;
      free_line(l);
    else:
      message "Anomalous projection plane; did you define it?";
    fi;
  fi;
  free_vec(v_a);
enddef;


% Distance between the observer and point |n| of object |name|
% Result is put in |dist|
vardef obs_distance(text dist)(expr name,n)=
  new_vec(v_a);
  define_current_point_offset_(name);% enables |pnt|
  dist:=vec_mod_(v_a,pnt(n),Obs);
  free_vec(v_a);
enddef;

%%\newpage
%%\title{Vector and point allocation}
% Allocation is done through a stack of vectors
numeric last_vec_;
last_vec_=0;

% vector allocation
% (this must not be a |vardef| because the vector |v| saved is not saved
% in this macro, but in the calling context)
def new_vec(text v)=
  save v;
  new_vec_(v);
enddef;

def new_vec_(text v)=
  v:=incr(last_vec_);
       %|message "Vector " & decimal (last_vec_+1) & " allocated";|
enddef;

let new_point = new_vec;
let new_point_ = new_vec_;

def new_points(text p)(expr n)=
  save p;
  numeric p[];
  for i:=1 upto n:new_point_(p[i]);endfor;
enddef;

% Free a vector
% A vector can only be freed safely when it was the last vector created.
def free_vec(expr i)=
  if i=last_vec_: last_vec_:=last_vec_-1;
     %|message "Vector " & decimal i & " freed";|
  else: errmessage("Vector " & decimal i & " can't be freed!");
  fi;
enddef;

let free_point = free_vec;

def free_points(text p)(expr n)=
  for i:=n step-1 until 1:free_point(p[i]);endfor;
enddef;

%%\title{Debugging}

def show_vec(expr t,i)=
  message "Vector " & t & "="
  & "(" & decimal vec[i]x & "," & decimal vec[i]y & ","
    & decimal vec[i]z & ")";
enddef;

% One can write |show_point("2",pnt_obj("obj",2));|
let show_point=show_vec;

def show_pair(expr t,zz)=
  message t & "=(" & decimal xpart(zz) & "," & decimal ypart(zz) & ")";
enddef;

%%\newpage
%%\title{Access to object features}
% |a| must be a string representing a class name, such as |"dodecahedron"|.
% |b| is the tail of a macro name.

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

def obj_points_(expr name)=
  obj_(obj_class_(name),"_points",name)
enddef;

def obj_faces_(expr name)=
  obj_(obj_class_(name),"_faces",name)
enddef;

vardef obj_point_offset_(expr name)=
  obj_(obj_class_(name),"_point_offset",name)
enddef;

vardef obj_face_offset_(expr name)=
  obj_(obj_class_(name),"_face_offset",name)
enddef;

def obj_class_(expr name)=obj_(name,"_class","") enddef;

%%\newpage
def define_point_offset_(expr name,o)=
  begingroup save n,tmpdef;
    string n,tmpdef;
    n=obj_class_(name) & "_point_offset" & name;
    expandafter numeric scantokens n;
    scantokens n:=last_point_offset_;
    last_point_offset_:=last_point_offset_+o;
    tmpdef="def " & obj_class_(name) & "_points" & name & 
      "=" & decimal o & " enddef";
    scantokens tmpdef;
  endgroup
enddef;

def define_face_offset_(expr name,o)=
  begingroup save n,tmpdef;
    string n,tmpdef;
    n=obj_class_(name) & "_face_offset" & name;
    expandafter numeric scantokens n;
    scantokens n:=last_face_offset_;
    last_face_offset_:=last_face_offset_+o;
    tmpdef="def " & obj_class_(name) & "_faces" & name & 
      "=" & decimal o & " enddef";
    scantokens tmpdef;
  endgroup
enddef;

def define_current_point_offset_(expr name)=
  save current_point_offset_;
  numeric current_point_offset_;
  current_point_offset_:=obj_point_offset_(name);
enddef;

def define_current_face_offset_(expr name)=
  save current_face_offset_;
  numeric current_face_offset_;
  current_face_offset_:=obj_face_offset_(name);
enddef;


%%\newpage
%%\title{Drawing an object}
% |name| is an object instance
vardef draw_obj(expr name)=
  save tmpdef;
  string tmpdef;
  current_obj:=name; 
  tmpdef="draw_" & obj_class_(name);
  project_obj(name);% compute screen coordinates
  save overflow; boolean overflow; overflow=false;
  for $:=1 upto obj_points_(name): 
    if z[ipnt_($)]=(too_big_,too_big_):overflow:=true;
      x[ipnt_($)] := 10; % so that the figure can be drawn anyway
      y[ipnt_($)] := 10;
      % why can't I write z[ipnt_($)]:=(10,10); ?
    fi;
    exitif overflow;
  endfor;
  if overflow:
    message "Figure has overflows";
    message "  (at least one point is not visible ";
    message "   and had to be drawn at a wrong place)";
  fi;
  scantokens tmpdef(name);
enddef;

%%\title{Normalization of an object}
% This macro translates an object so that a list of vertices is centered
% on the origin, and the last vertex is put on a sphere whose radius is 1.
% |name| is the name of the object and |vertices| is a list
% of points whose barycenter will define the center of the object.
% (|vertices| need not be the list of all vertices)
vardef normalize_obj(expr name)(text vertices)=
  save nvertices,last;
  nvertices=0;
  new_vec(v_a);vec_def_(v_a,0,0,0)
  forsuffixes $=vertices:
    vec_sum_(v_a,v_a,pnt($));
    nvertices:=nvertices+1;
    last:=$;
  endfor;
  vec_mult_(v_a,v_a,-1/nvertices);
  translate_obj(name,v_a);% object centered on the origin
  scale_obj(name,1/vec_mod(last));
  free_vec(v_a);
enddef;


%%\newpage
%%\title{General definitions}
% Vector arrays
numeric vec[]x,vec[]y,vec[]z;

% Reference vectors $\vec{0}$, $\vec{\imath}$, $\vec{\jmath}$ and $\vec{k}$ 
% and their definition
new_vec(vec_null);new_vec(vec_I);new_vec(vec_J);new_vec(vec_K);
vec_def_(vec_null,0,0,0);
vec_def_(vec_I,1,0,0);vec_def_(vec_J,0,1,0);vec_def_(vec_K,0,0,1);
numeric point_null;
point_null=vec_null;

% Observer 
new_point(Obs);
% default value:
set_point_(Obs,0,0,20); 

% Observer's vectors
new_vec(ObsI_);new_vec(ObsJ_);new_vec(ObsK_);
% default values: 
vec_def_vec_(ObsI_,vec_I);
vec_def_vec_(ObsJ_,vec_J);
vec_def_vec_(ObsK_,vec_K);

new_vec(IObsI_);new_vec(IObsJ_);new_vec(IObsK_);

% These vectors will be vectors of the projection plane,
% in case of oblique projections:
new_vec(ProjK_);new_vec(ProjJ_); % there is no |ProjI_|

% This will be the center of the projection plane, in oblique projections
new_point(ObliqueCenter_);


% distance observer/plane (must be $>0$)
numeric Obs_dist; % represents |Obs_dist| $\times$ |drawing_scale|
% default value:
Obs_dist=2; % means |Obs_dist| $\times$ |drawing_scale|

% current object being drawn
string current_obj;

% kind of projection: 0 for linear (or central) perspective, 1 for parallel,
% 2 for oblique projection
% (default is 0)
numeric projection_type;
projection_type:=0;

% Definition of a projection plane (only used in oblique projections)
%
new_plane_(projection_plane)(1,1,1); % the initial value is irrelevant

% Definition of a projection direction (only used in oblique projections)
new_line_(projection_direction)(1,1);  % the initial value is irrelevant

% this positions the observer at vector |p| (the point observed)
% + |d| (distance) * (k-(i+j))
def isometric_projection(expr i,j,k,p,d,phi)=
  trimetric_projection(i,j,k,1,1,1,p,d,phi);
enddef;

% this positions the observer at vector |p| (the point observed)
% + |d| (distance) * (ak-(i+j))
def dimetric_projection(expr i,j,k,a,p,d,phi)=
  trimetric_projection(i,j,k,1,1,a,p,d,phi);
enddef;

% this positions the observer at vector |p| (the point observed)
% + |d| (distance) * (k-(i+j))
% |a|, |b| and |c| are multiplicative factors to vectors |i|, |j| and |k|
vardef trimetric_projection(expr i,j,k,a,b,c,p,d,phi)=
  save v_a,v_b,v_c;
  new_vec(v_a);new_vec(v_b);new_vec(v_c);
  vec_mult_(v_a,i,a);vec_mult_(v_b,j,b);vec_mult_(v_c,k,c);
  vec_sum_(Obs,v_a,v_b);
  vec_diff_(Obs,v_c,Obs);
  vec_mult_(Obs,Obs,d);
  vec_sum_(Obs,Obs,p);
  point_of_view_abs(p,phi);
  projection_type:=1;
  free_vec(v_c);free_vec(v_b);free_vec(v_a);
enddef;

% |hor| is an horizontal plane (in the sense that it will represent
% the horizontal for the observer)
% |p| is the point in space that the observer targets (center of screen)
% |a| is an angle (45 degrees corresponds to cavalier drawing)
% |b| is an angle (see examples defined below)
% |d| is the distance of the observer
vardef oblique_projection(text hor)(expr p,a,b,d)=
  save _l,v_a,v_b,v_c,xxx_,obsJangle_;
  new_vec(v_a);new_vec(v_b);new_vec(v_c);
  % we first compute a horizontal line:
  new_line_(_l)(1,1);
  if def_inter_l_pl_pl(_l)(hor)(projection_plane):
    vec_diff_(v_a,_l2,_l1); % horizontal vector
    % then, we find a normal to the projection plane:
    def_normal_p_(v_b)(projection_plane);
    % complete the line and the vector by a third vector (=vertical)
    vec_prod_(v_c,v_a,v_b);
    % we make |v_a| a copy of |v_b| since we no longer need |v_b|
    vec_def_vec_(v_a,v_b);
    % we rotate |v_b| by an angle |a| around |v_c|
    vec_rotate_(v_b,v_c,a);
    % we rotate |v_b| by an angle |b| around |v_a|
    vec_rotate_(v_b,v_a,b);
    % we put the observer at the distance |d| of |p| in
    % the direction of |v_b|:
    vec_unit_(v_b,v_b);
    vec_mult_(v_b,v_b,d);vec_sum_(Obs,p,v_b);
    % We now have to make sure that point |p| and point |Obs|
    % are on different sides of the projection plane. For this,
    % we compute two dot products:
    new_vec(v_d);new_vec(v_e);
    vec_diff_(v_d,p,_l1);vec_diff_(v_e,Obs,_l1);
    if vec_dprod_(v_d,v_a)*vec_dprod_(v_e,v_a)>=0:
      % |p| and |Obs| are on the same side of the projection plane
      % |Obs| needs to be recomputed.
      vec_mult_(v_b,v_b,-1);
      vec_sum_(Obs,p,v_b);
    fi;
    free_vec(v_e);free_vec(v_d);
    projection_type:=2;    % needs to be set before |point_of_view_abs|
    point_of_view_abs(p,90); % this computes |ObliqueCenter_|
    % and now, make sure the vectors defining the observer are right:
    % Create the plane containing lines _l and projection_direction
    %  (defined by point_of_view_abs):
    new_plane_(xxx_)(1,1,1);
    def_plane_pl_l_l(xxx_)(_l)(projection_direction);
    % Compute the angle of |ObsK_| with this plane:
    obsJangle_=vangle_v_pl_(ObsK_)(xxx_);
    % rotate |ObsJ_| and |ObsK_| by |obsJangle_| around |ObsI_|
    vec_rotate_(ObsJ_,ObsI_,obsJangle_);
    vec_rotate_(ObsK_,ObsI_,obsJangle_);
    if abs(vangle_v_pl_(ObsK_)(xxx_))>1: % the rotation was done
                                         % in the wrong direction
      vec_rotate_(ObsJ_,ObsI_,-2obsJangle_);
      vec_rotate_(ObsK_,ObsI_,-2obsJangle_);
    fi;
    % |vec_rotate_(ObsJ_,ObsI_,45);| % planometric test
    % |vec_rotate_(ObsK_,ObsI_,45);| % planometric test
    free_plane(xxx_);
    % and now, |ProjJ_| and |ProjK_| must be recomputed:
    project_vectors(ProjK_,ProjJ_)(ObsK_,ObsJ_)%
               (projection_direction)(projection_plane);
  else:
    message "Error: the ``horizontal plane'' cannot be";
    message "  parallel to the projection plane.";
  fi;
  free_line(_l);
  free_vec(v_c);free_vec(v_b);free_vec(v_a);
enddef;

% These two are the most common values for the third parameter
% of |oblique_projection|
numeric CAVALIER;CAVALIER=45;
numeric CABINET;CABINET=angle((1,.5)); % atn(.5)

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

% Observer's orientation, defined by three angles
numeric Obs_psi,Obs_theta,Obs_phi; 
% default value:
Obs_psi=0;Obs_theta=90;Obs_phi=0;

% This array relates an absolute object point number to the
% absolute point number (that is, to the |vec| array).
% The absolute object point number is the rank of a point
% with respect to all object points. The absolute point number
% considers in addition the extra points, such as |Obs|, which do
% not belong to an object.
% If |i| is an absolute object point number, |points_[i]|
% is the absolute point number.
numeric points_[]; 

% |name| is the name of an object instance
% |npoints| is its number of defining points
def new_obj_points(expr name,npoints)= 
  define_point_offset_(name,npoints);define_current_point_offset_(name);
  for i:=1 upto obj_points_(name):new_point_(pnt(i));endfor;
enddef;

% |name| is the name of an object instance
% |nfaces| is its number of defining faces
def new_obj_faces(expr name,nfaces)= 
  define_face_offset_(name,nfaces);define_current_face_offset_(name);
  redefine_obj_name_(name,current_face_offset_+nfaces);
enddef;

%%\newpage
% Absolute point number corresponding to object point number |i|
% This macro must only be used within the function defining an object
% (such as |def_cube|) or the function drawing an object (such as
% |draw_cube|).
def ipnt_(expr i)=i+current_point_offset_ enddef;
def pnt(expr i)=points_[ipnt_(i)] enddef;

def face(expr i)=(i+current_face_offset_) enddef;

% Absolute point number corresponding to local point |n|
% in object instance |name|
vardef pnt_obj(expr name,n)=
  points_[n+obj_point_offset_(name)]
  %hide(define_current_point_offset_(name);) pnt(n) % HAS SIDE EFFECTS
enddef;

% Absolute face number corresponding to local face |n|
% in object instance |name|
vardef face_obj(expr name,n)=
  (n+obj_face_offset_(name))
  %hide(define_current_face_offset_(name);) face(n) % HAS SIDE EFFECTS
enddef;


% Scale
numeric drawing_scale;
drawing_scale=2cm;

% Color
% This function is useful when a color is expressed in hexadecimal.
% This does the opposite from |tohexcolor|
def hexcolor(expr s)=
  (hex(substring (0,2) of s)/255,hex(substring (2,4) of s)/255,
    hex(substring (4,6) of s)/255)
enddef;

% Convert a color triple into a hexadecimal color string.
% |rv|, |gv| and |bv| are values between 0 and 1.
% This does the opposite from |hexcolor|
vardef tohexcolor(expr rv,gv,bv)=
  save dig;numeric dig[];
  hide(
    dig2=floor(rv*255);dig1=floor((dig2)/16);dig2:=dig2-16*dig1;
    dig4=floor(gv*255);dig3=floor((dig4)/16);dig4:=dig4-16*dig3;
    dig6=floor(bv*255);dig5=floor((dig6)/16);dig6:=dig6-16*dig5;
    for i:=1 upto 6:
      if dig[i]<10:dig[i]:=dig[i]+48;
      else:dig[i]:=dig[i]+87;
      fi;
    endfor;
  )
  char(dig1)&char(dig2)&char(dig3)&char(dig4)&char(dig5)&char(dig6)
enddef;

% Conversions

% Returns a string encoding the integer |n| as follows:
% if $n=10*a+b$ with $b<10$,
%   |alphabetize|(|n|)=|alphabetize|(|a|) |&| |char (65+b)|
% For instance, alphabetize(3835) returns "DIDF"
% This function is useful in places where digits are not allowed.
def alphabetize(expr n)=
  if (n>9):
    alphabetize(floor(n/10)) & fi
  char(65+n-10*floor(n/10))
enddef;

% Filling and contours
boolean filled_faces,draw_contours;
filled_faces=true;
draw_contours=true;
numeric contour_width; % thickness of contours
contour_width=1pt;
color contour_color; % face contours
contour_color=black;

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


% Object offset (the points defining an object are arranged
% in a single array, and the objects are easier to manipulate
% if the point numbers are divided into a number and an offset).
numeric last_point_offset_,last_face_offset_;
last_point_offset_=0;last_face_offset_=0;

endinput