-5 4 setyrange
-6 5 setxrange

-5 8 4 SetCamPos
0 0 0 SetCamView
50 setxunit
GetCamPos setlightsrc
2 setlightintensity
2 setlinejoin

%% syntaxe : M I k hompoint3d --> M' tq vect(IM') = k vect(IM)
/hompoint3d {
3 dict begin
   /k exch def
   /I defpoint3d
   /M defpoint3d
   I M vecteur3d
   k mulv3d
   I addv3d
end
} def

%% syntaxe : solid coeff i s@lidaffineface
/solidchanfreine {
10 dict begin
   /coeff exch def
   /solid exch def
   /result newsolid def
   solid issolid not {
      (Erreur : mauvais type d argument dans solidchanfreine) ==
      quit
   } if
   /n solid solidnombresommets def
   /nf solid solidnombrefaces def

   %% ajout des faces reduites
   0 1 nf 1 sub {
      /i exch def
      /Fsommets solid i solidgetsommetsface def
      /Findex solid i solidgetface def
      /ns Fsommets length 3 idiv def
      /couleurfaceorigine solid i solidgetfcolor def
      Fsommets isobarycentre3d /G defpoint3d
      %% on ajoute les nouveaux sommets
      /Sindex [] def
      0 1 ns 1 sub {
         /j exch def
         /Sindex [ Sindex aload pop
            Fsommets j getp3d /M defpoint3d
            result M G coeff hompoint3d solidaddsommet
         ] store
      } for
      %% Sindex contient les indices des nouveaux sommets
      result Sindex solidaddface
   } for

etape1 {
   %% ajout des faces rectangulaires entre faces d'origines adjacentes
   %% pour chaque face de depart
   0 1 nf 2 sub {
      /i exch def
      /F solid i solidgetface def
      /Fres result i solidgetface def
      %% pour chaque arete de la face
      0 1 F length 1 sub {
         /j exch def
         /trouve false def
         /indice1 F j get def
         /indice2 F j 1 add F length mod get def
         /a1 j def
         /a2 j 1  add F length mod def
         %% on regarde toutes les autres faces
         i 1 add 1 nf 1 sub {
            /k exch def
            /Ftest solid k solidgetface def
            indice1 Ftest in {pop true} {false} ifelse
            indice2 Ftest in {pop true} {false} ifelse
            and {
               /indiceFadj k def
               indice1 Ftest in pop /k1 exch def
               indice2 Ftest in pop /k2 exch def
               /trouve true def
            exit
            } if
         } for
         trouve {
            /Fadj solid indiceFadj solidgetface def
            result [
               Fres a1 get
               result indiceFadj solidgetface k1 get
               result indiceFadj solidgetface k2 get
               Fres a2 get
            ] solidaddface
         } if
      } for
   } for
} if

etape2 {
   %% pour chaque face
   0 1 nf 2 sub {
      /i exch def
      /F solid i solidgetface def
      %% et pour chaque sommet de cette face
      0 1 F length 1 sub {
         /j exch def
         /k F j get def
         solid k solidfacesadjsommet /adj exch def
         %% adj est le tableau des indices des faces adjacentes
         %% au sommet d'indice k
         %% rque : toutes les faces d'indice strict inferieur a i
         %% sont deja traitees
         %% Pour chaque face adjacente, on repere l'indice du sommet concerne dans
         %% la face
         adj min i lt not {
            /indadj [] def
            0 1 adj length 1 sub {
               /m exch def
               k solid adj m get solidgetface in {
                  /ok exch def
                  /indadj [indadj aload pop ok] store
               } if
            } for
         
            /aajouter [
               0 1 adj length 1 sub {
                  /m exch def
                  result adj m get solidgetface indadj m get get
               } for
            ] def

	    %% la table des sommets
	    [0 1 aajouter length 1 sub {
	       /m exch def
	       result aajouter m get solidgetsommet
	    } for]
	    solid k solidgetsommet %% le point indiquant la direction de la normale
	    ordonnepoints3d
	    /indicestries exch def

	    result [
	       0 1 indicestries length 1 sub {
	          /m exch def
		  aajouter indicestries m get get
	       } for
	    ] solidaddface
	 } if
      } for
   } for
} if   
   result
end
} def


%% syntaxe : array1 M ordonnepoints3d --> array2
%% array1 = tableau de points 3d coplanaires (plan P)
%% M = point3d indiquant la direction de la normale a P
%% array2 = les indices des points de depart, ranges dans le 
%% sens trigo par rapport a la normale
/ordonnepoints3d {
5 dict begin
   /M defpoint3d
   /table exch def
   table isobarycentre3d /G defpoint3d
   %% calcul de la normale
   table 0 getp3d /ptref defpoint3d
   table 1 getp3d /A defpoint3d
   G ptref vecteur3d
   G A vecteur3d
   vectprod3d /vecteurnormal defpoint3d
   vecteurnormal G M vecteur3d scalprod3d 0 lt {
      vecteurnormal -1 mulv3d /vecteurnormal defpoint3d
   } if
   %% la table des angles
   table duparray exch pop
   {1 dict begin
      /M defpoint3d
      G ptref vecteur3d
      G M vecteur3d
      vecteurnormal angle3doriente
   end} papply3d
   [0 1 table length 3 idiv 1 sub {} for]
   exch doublebubblesort pop
end
} def

%% syntaxe : solid i solidfacesadjsommet --> array
%% array est le tableau des indices des faces adjacentes au
%% sommet d indice i
/solidfacesadjsommet {
10 dict begin
   /indicesommet exch def
   /solid exch def
   /result [] def
   %% pour chaque face
   0 1 solid solidnombrefaces 1 sub {
      /i exch def
      /F solid i solidgetface def
      indicesommet F in {
         pop
	 /result [result aload pop i] store
      } if
   } for
   result
end
} def

/k .4 def

/etape1 true def
/etape2 false def

3 newcube
dup .7  solidchanfreine
%dup .3 solidtronque
%dup .8 solidaffine
%dup .55 solidchanfreine
%dup [.5 .5] solidputhuecolors
dup videsolid
dup (rouge) (jaune) solidputcolors
dup drawsolid**
%solidshowsommets