mp-solid.mp

   1 %% =============================================================================
   2 %% Fichier : mp-solid.mp
   3 %% =============================================================================
   4 %% Le fichier spatial d'Anthony Phan m'a permis de mettre certaines choses au
   5 %% clair. ;-)
   6 %% 14/08/2008
   7
   8 prologues:=2;
   9
  10 %% Constantes ------------------------------------------------------------------
  11 u  := 1cm;
  12 pi := 3.141592654;
  13 c  := 57.29578;                               % conversion d'un radian en degrés
  14
  15 color rouge, vert, bleu, jaune, noir, blanc, orange, rose, violet, ciel,
  16    cielfonce, orangevif, gris, marron;
  17 rouge     = (1,0,0);
  18 bleu      = (0,0,1);
  19 noir      = (0,0,0);
  20 blanc     = (1,1,1);
  21 orange    = (1,0.5,0);
  22 violet    = blanc-vert;
  23 rose      = (1,0.7,0.7);
  24 cielfonce = 0.9*(0.25,1,1);
  25 ciel      = bleu+vert;
  26 orangevif = (1,0.25,0.1);
  27 vert      = (0,1,0);
  28 jaune     = blanc-bleu;
  29 gris      = 0.8*white;
  30
  31 input format;
  32 input marith;
  33 input sarith;
  34 %input donymodule;
  35
  36 input objets;
  37
  38 color Sommet[];
  39
  40 %% Anthony Phan
  41 vardef Norm primary z =
  42   abs (abs(Xpart z, Ypart z), Zpart z)
  43 enddef;
  44
  45 let Xpart = redpart;
  46 let Ypart = greenpart;
  47 let Zpart = bluepart;
  48 %
  49
  50 string typerepre,pointilles;
  51 typerepre := "persp";
  52
  53 vardef Initialisation(expr r,t,p,d)=
  54   Rho:=r;
  55   Theta:=t;
  56   Phi:=p;
  57   DE:=d;
  58   Aux1:=sind(Theta);
  59   Aux2:=sind(Phi);
  60   Aux3:=cosd(Theta);
  61   Aux4:=cosd(Phi);
  62   Aux5:=Aux3*Aux2;
  63   Aux6:=Aux1*Aux2;
  64   Aux7:=Aux3*Aux4;
  65   Aux8:=Aux1*Aux4;
  66   pointilles:="oui";
  67   Lumiere:=Oeil;
  68 enddef;
  69
  70 vardef Oeil=(Rho*Aux7,Rho*Aux8,Rho*Aux2)
  71 enddef;
  72
  73 vardef sin(expr t) = sind(c*t) enddef;
  74
  75 vardef cos(expr t) = cosd(c*t) enddef;
  76
  77 vardef tan(expr t) = sin(t)/cos(t) enddef;
  78
  79 vardef exp(expr x) = mexp(256)**x enddef;
  80
  81 vardef Exp primary x = mexp(256)**x enddef;
  82
  83 vardef ln(expr t) = mlog(t)/256 enddef;
  84
  85 vardef log(expr t) = ln(t)/ln(10) enddef;
  86
  87 vardef ch(expr x)=(exp(x)+exp (-x))/2 enddef;
  88
  89 vardef sh(expr x)=(exp(x)-exp(-x))/2 enddef;
  90
  91 vardef th(expr x)=sh(x)/ch(x) enddef;
  92
  93 vardef arcsin(expr x)=%Définition mathématique en radian
  94   pi*angle((sqrt(1-x**2),x))/180
  95 enddef;
  96
  97 vardef arccos(expr x)=%Définition mathématique en radian
  98   pi*angle((x,sqrt(1-x**2)))/180
  99 enddef;
 100
 101 vardef arctan(expr x)=arcsin(x/(1++x))
 102 enddef;
 103
 104 numeric satu,lum;
 105 satu:=0.45;
 106 lum:=1;
 107
 108 %Coordonnées dans le repère Oeil
 109 vardef GCoord(expr N)=
 110   (-Xpart(N)*Aux1+Ypart(N)*Aux3,-Xpart(N)*Aux5-Ypart(N)*Aux6+Zpart(N)*Aux4,-Xpart(N)*Aux7-Ypart(N)*Aux8-Zpart(N)*Aux2+Rho)
 111 enddef;
 112
 113 unit:=1;%pour les mises à l'échelle :) Merci pst-solides3d
 114
 115 vardef Projette(expr M)=
 116   %if typerepre="proj":
 117   %  unit*(DE*(Xpart(GCoord(M)/Zpart(GCoord(M))),(Ypart(GCoord(M))/Zpart(GCoord(M)))))
 118   %elseif typerepre="persp":
 119     unit*(DE*(Xpart(GCoord(M)),Ypart(GCoord(M))))
 120   %fi
 121 enddef;
 122
 123 vardef TraceAxes=
 124   color Origine,Unitex,Unitey,Unitez;
 125   Origine=(0,0,0);
 126   Unitex=(5,0,0);
 127   Unitey=(0,5,0);
 128   Unitez=(0,0,5);
 129   drawoptions(dashed dashpattern(on 12bp off 6bp on 3bp off 6bp));
 130   drawarrow Projette(Origine)--Projette(Unitex) withcolor blue;
 131   drawarrow Projette(Origine)--Projette(Unitey) withcolor vert;
 132   drawarrow Projette(Origine)--Projette(Unitez);
 133   drawoptions();
 134 enddef;
 135
 136 vardef TraceAxesD(expr xd,yd,zd)=
 137   drawoptions(dashed dashpattern(on12bp off6bp on3bp off6bp));;
 138   drawarrow Projette((0,0,0))--Projette((xd,0,0));
 139   drawarrow Projette((0,0,0))--Projette((0,yd,0));
 140   drawarrow Projette((0,0,0))--Projette((0,0,zd));
 141   drawoptions();
 142   label(btex $x$ etex,Projette((xd+0.3,0,0)));
 143   label(btex $y$ etex,Projette((0,yd+0.3,0)));
 144   label(btex $z$ etex,Projette((0,0,zd+0.3)));
 145 enddef;
 146
 147 vardef TraceAxesDD(expr xd,yd,zd,xf,yf,zf)=
 148   drawoptions(dashed evenly);
 149   draw Projette((0,0,0))--Projette((xd,0,0));
 150   draw Projette((0,0,0))--Projette((0,yd,0));
 151   draw Projette((0,0,0))--Projette((0,0,zd));
 152   drawoptions();
 153   drawarrow Projette((xd,0,0))--Projette((xf,0,0));
 154   drawarrow Projette((0,yd,0))--Projette((0,yf,0));
 155   drawarrow Projette((0,0,zd))--Projette((0,0,zf));
 156   label(btex $x$ etex,Projette((xf+0.3,0,0)));
 157   label(btex $y$ etex,Projette((0,yf+0.3,0)));
 158   label(btex $z$ etex,Projette((0,0,zf+0.3)));
 159 enddef;
 160
 161 primarydef u Vectprod v =
 162   (Ypart u * Zpart v - Zpart u * Ypart v,
 163     Zpart u * Xpart v - Xpart u * Zpart v,
 164     Xpart u * Ypart v - Ypart u * Xpart v)
 165 enddef;
 166
 167 vardef Normal(expr vecun,vecde,vectr)=
 168   save aa;
 169   color aa;
 170   P1:=redpart(vecde-vecun);
 171   P2:=greenpart(vecde-vecun);
 172   P3:=bluepart(vecde-vecun);
 173   Q1:=redpart(vectr-vecun);
 174   Q2:=greenpart(vectr-vecun);
 175   Q3:=bluepart(vectr-vecun);
 176   aa=(P2*Q3-Q2*P3,P3*Q1-Q3*P1,P1*Q2-Q1*P2);
 177   aa
 178 enddef;
 179
 180 vardef ProduitScalaire(expr wec,mor)=
 181   %Mexp(Mlog redpart(wec) Mmul Mlog redpart(mor))+Mexp(Mlog greenpart(wec) Mmul Mlog greenpart(mor))+Mexp(Mlog bluepart(wec) Mmul Mlog bluepart(mor))
 182   Xpart(wec)*Xpart(mor)+Ypart(wec)*Ypart(mor)+Zpart(wec)*Zpart(mor)
 183 enddef;
 184 %pour les rotations et translations
 185
 186 vardef RotX(expr ptx)=
 187   (Xpart(ptx),Ypart(ptx)*cosd(angx)-Zpart(ptx)*sind(angx),Ypart(ptx)*sind(angx)+Zpart(ptx)*cosd(angx))
 188 enddef;
 189
 190 vardef RotY(expr ptx)=
 191   (Xpart(ptx)*cosd(angy)+Zpart(ptx)*sind(angy),Ypart(ptx),-Xpart(ptx)*sind(angy)+Zpart(ptx)*cosd(angy))
 192 enddef;
 193
 194 vardef RotZ(expr ptx)=
 195   (Xpart(ptx)*cosd(angz)-Ypart(ptx)*sind(angz),Xpart(ptx)*sind(angz)+Ypart(ptx)*cosd(angz),Zpart(ptx))
 196 enddef;
 197
 198 vardef RotXYZ(expr ptx)=
 199   RotZ(RotY(RotX(ptx)))
 200 enddef;
 201
 202 angx:=0;
 203 angy:=0;
 204 angz:=0;
 205 color TR;
 206 TR:=(0,0,0);
 207 %pour les tubes
 208 vardef VT(expr t)=Fp(t)/Norm(Fp(t))
 209 enddef;
 210
 211 vardef VN(expr t)=Fd(t)/Norm(Fd(t))
 212 enddef;
 213
 214 vardef VBN(expr t)=VT(t) Vectprod VN(t)
 215 enddef;
 216
 217 path feuillet;
 218 numeric _tfig,_nfig;
 219 pair coinbg,coinbd,coinhd,coinhg;
 220 _nfig:=0;
 221
 222 def feuille(expr xa,ya,xb,yb) =
 223   feuillet := (xa,ya)--(xa,yb)--(xb,yb)--(xb,ya)--cycle;
 224   coinbg := (xa,ya);
 225   coinbd := (xb,ya);
 226   coinhd := (xb,yb);
 227   coinhg := (xa,yb);
 228   %modifié le 29.09.04
 229   z.so=(xpart(coinbg/1cm),ypart(coinbg/1cm));
 230   z.ne=(xpart(coinhd/1cm),ypart(coinhd/1cm));
 231   %fin modification
 232   extra_endfig := "clip currentpicture to feuillet;" & extra_endfig;
 233 enddef;
 234
 235 def figureespace(expr xa,ya,xb,yb) =
 236   _nfig:=_nfig+1;
 237   beginfig(_nfig);
 238     feuille(xa,ya,xb,yb);
 239     _tfig:= if (xb-xa)>(yb-ya): xb-xa else: yb-ya fi;
 240         _tfig:=2*_tfig;
 241 enddef;
 242
 243 def finespace=
 244 endfig;
 245 enddef;
 246
 247 def QS(expr ndeb,nfin)=
 248   begingroup
 249     save v,m,x;
 250     if ndeb<nfin:
 251       v:=ALT[cpt[ndeb]];
 252       m:=ndeb;
 253       for i=(ndeb+1) upto nfin:
 254         if ALT[cpt[i]]<v:
 255           m:=m+1;
 256           x:=cpt[m];cpt[m]:=cpt[i];cpt[i]:=x;
 257         fi
 258       endfor;
 259       x:=cpt[m];cpt[m]:=cpt[ndeb];cpt[ndeb]:=x;
 260       QS(ndeb,m-1);
 261       QS(m+1,nfin);
 262     fi
 263   endgroup
 264 enddef;
 265
 266 boolean Pointilles;
 267 Pointilles:=false;
 268
 269 boolean Vue[];
 270 color Fc[].iso;
 271 boolean Creux;
 272 Creux:=false;
 273
 274 vardef DessineObjetNew=
 275   numeric ALT[];
 276   %on détermine les zmax dans le repère spatial de l'écran
 277   for k=1 upto apj:
 278     Fc[k].iso:=(0,0,0) for l=1 upto ns[k][0]:
 279       +Fc[k][l]
 280     endfor;
 281     Fc[k].iso:=Fc[k].iso/ns[k][0];
 282   endfor;
 283   for k=1 upto apj:
 284     %ALT[k]=-Zpart(GCoord(Fc[k].iso));
 285     zmax:=infinity;
 286     for l=1 upto ns[k][0]:
 287       t:=-Zpart(GCoord(Fc[k][l]));
 288       if t<zmax:
 289         zmax:=t;
 290       fi;
 291     endfor;
 292     ALT[k]=zmax;
 293     cpt[k]:=k;
 294   endfor;
 295   %On trie suivant les zmax
 296   QS(1,apj);
 297   %On dessine
 298   if Pointilles:
 299     for k=1 upto apj:
 300       draw for l=1 upto ns[cpt[k]][0]:
 301         Projette(Fc[cpt[k]][l])--
 302       endfor
 303       cycle if Vue[cpt[k]]=false: dashed evenly fi;
 304     endfor;
 305   else:
 306     for k=1 upto apj:
 307       if Creux=false:
 308         if Vue[cpt[k]]:
 309           fill for l=1 upto ns[cpt[k]][0]:
 310             Projette(Fc[cpt[k]][l])--
 311           endfor
 312           cycle withcolor
 313           if arcenciel:
 314             lumin(cpt[k])*Hsvtorgb(((cpt[k]/apj)*360,satu,lum))
 315           else:
 316             lumin(cpt[k])*outcolor
 317           fi;
 318           if traits:
 319             draw for l=1 upto ns[cpt[k]][0]:
 320               Projette(Fc[cpt[k]][l])--
 321             endfor
 322             cycle;
 323           fi;
 324         fi;
 325       else:
 326         fill for l=1 upto ns[cpt[k]][0]:
 327           Projette(Fc[cpt[k]][l])--
 328         endfor
 329         cycle withcolor
 330         if Vue[cpt[k]]:
 331           if arcenciel:
 332             lumin(cpt[k])*Hsvtorgb(((cpt[k]/apj)*360,satu,lum))
 333           else:
 334             lumin(cpt[k])*outcolor
 335           fi
 336         else:
 337           lumin(cpt[k])*incolor
 338         fi;
 339         if traits:
 340           draw for l=1 upto ns[cpt[k]][0]:
 341             Projette(Fc[cpt[k]][l])--
 342           endfor
 343           cycle;
 344         fi;
 345       fi;
 346     endfor;
 347   fi;
 348 enddef;
 349
 350 color Fc[][];%représente les sommets des différentes faces;
 351
 352 invnormale=1;%Parfois dans la lecture des fichiers OFF, il est nécessaire d'opposer la normale pour obtenir un bon intérieur-extérieur
 353 boolean OFF,OBJ;%pour utiliser la lumière "correctement"
 354 OFF:=false;
 355 OBJ:=false;
 356
 357 vardef LectureOFF(expr nomsolide)=
 358   OFF:=true;
 359   %Détermination du nombre de sommets et de faces.
 360   string s_;
 361   s_=readfrom nomsolide;
 362   string ss[];
 363   if s_<>EOF:
 364     ss1 := loptok s_;
 365     t_ := if ss1="%": 0 else: 1 fi;
 366     forever:
 367       ss[incr t_] := loptok s_;
 368       exitif ss[t_]="";
 369     endfor
 370   else: false
 371   fi;
 372   NbS:=round(Mexp Mlog_str ss1);
 373   NF:=round(Mexp Mlog_str ss2);
 374   message("Il y a "&decimal(NbS)&" sommets.");
 375   message("Il y a "&decimal(NF)&" faces au total.");
 376   %Détermination des coordonnées des sommets
 377   message("Création des sommets.");
 378   s_:=readfrom nomsolide;
 379   if debut=0:
 380     for k=0 upto NbS-1:
 381       s_:=readfrom nomsolide;
 382       if s_<>EOF:
 383         ss1 := loptok s_;
 384         n_ := if ss1="%": 0 else: 1 fi;
 385         forever:
 386           ss[incr n_] := loptok s_;
 387           exitif ss[n_]="";
 388         endfor
 389       else: false
 390       fi;
 391       Sommet[k]:=(Mexp ((Mlog_str ss1) Mdiv (Mlog echelle)),Mexp ((Mlog_str ss3) Mdiv (Mlog echelle)),Mexp ((Mlog_str ss2) Mdiv (Mlog echelle)));
 392       %Sommet[k]:=(Mexp (Mlog_str ss1),Mexp (Mlog_str ss3),Mexp (Mlog_str ss2))/echelle;
 393       %Sommet[k]:=(Mexp (Mlog_str ss1)/echelle,Mexp (Mlog_str ss3)/echelle,Mexp (Mlog_str ss2)/echelle);
 394     endfor;
 395   else:
 396     for k=1 upto NbS:
 397       s_:=readfrom nomsolide;
 398       if s_<>EOF:
 399         ss1 := loptok s_;
 400         n_ := if ss1="%": 0 else: 1 fi;
 401         forever:
 402           ss[incr n_] := loptok s_;
 403           exitif ss[n_]="";
 404         endfor
 405       else: false
 406       fi;
 407       Sommet[k]:=(Mexp ((Mlog_str ss1) Mdiv (Mlog echelle)),Mexp ((Mlog_str ss3) Mdiv (Mlog echelle)),Mexp ((Mlog_str ss2) Mdiv (Mlog echelle)));
 408       %Sommet[k]:=(Mexp (Mlog_str ss1),Mexp (Mlog_str ss3),Mexp (Mlog_str ss2))/echelle;
 409       %Sommet[k]:=(Mexp (Mlog_str ss1)/echelle,Mexp (Mlog_str ss3)/echelle,Mexp (Mlog_str ss2)/echelle);
 410     endfor;
 411   fi;
 412   message("Création des faces.");
 413   %Détermination des faces
 414   apj:=0;color cc,dd;
 415   nbfvues:=0;
 416   for nf=-4000 upto (-4000+NF)-1:
 417     s_:=readfrom nomsolide;
 418      if s_<>EOF:
 419       ss1 := loptok s_;
 420       n_ := if ss1="%": 0 else: 1 fi;
 421       forever:
 422         ss[incr n_] := loptok s_;
 423         exitif ss[n_]="";
 424       endfor
 425     else: false
 426     fi;
 427     apj:=apj+1;
 428     ns[apj][0]:=Mexp Mlog_str ss1;%pour savoir le nb de sommets par face
 429     for nl=1 upto ns[apj][0]:
 430       Fc[apj][nl]:=Sommet[round(Mexp Mlog_str ss[nl+1])];
 431     endfor;
 432     dd:=Oeil-Fc[apj][1];
 433     cc:=invnormale*Normal(Fc[apj][1],Fc[apj][2],Fc[apj][3]);
 434     if (ProduitScalaire(dd,cc)>=0):
 435       Vue[apj]:=true;
 436       nbfvues:=nbfvues+1;
 437     else:
 438       if Creux=true:
 439         Vue[apj]:=false;
 440       else:
 441         apj:=apj-1;
 442       fi;
 443     fi;
 444   endfor;
 445   message("Faces vues déterminees : il y en a "&decimal(nbfvues)&".");
 446   closefrom nomsolide;
 447   DessineObjetNew;
 448 enddef;
 449
 450 vardef LectureOBJ(expr nomfichier)=
 451   OBJ:=true;
 452   string s_;
 453   string ss[];
 454   nbss:=1;
 455   apj:=0;
 456   forever:
 457     s_:=readfrom nomfichier;
 458     if s_<>EOF:
 459       ss0 := loptok s_;
 460       if ss0="v":
 461         n_:=0;
 462         forever:
 463           ss[incr n_] := loptok s_;
 464           exitif ss[n_]="";
 465         endfor
 466         Sommet[nbss]:=(Mexp((Mlog_str ss1) Mdiv (Mlog echelle)),Mexp ((Mlog_str ss3) Mdiv (Mlog echelle)),Mexp ((Mlog_str ss2) Mdiv (Mlog echelle)));
 467         nbss:=incr nbss;
 468       elseif ss0="f":
 469         n_:=0;
 470         forever:
 471           ss[incr n_] := loptok s_;
 472           exitif ss[n_]="";
 473         endfor;
 474         ns[apj][0]:=n_-1;
 475         for k=1 upto ns[apj][0]:
 476           if invnormale=1:
 477             Fc[apj][ns[apj][0]-k+1] := Sommet[round(Mexp(Mlog_str ss[k]))]
 478             %if unknown OTFc.@[apj][OTFc.@[apj].nb-k+1]:
 479           %  show OTFc.@[apj][OTFc.@[apj].nb-k+1];
 480           %fi;
 481           else:
 482             Fc[apj][k] := Sommet[round(Mexp(Mlog_str ss[k]))]
 483             %if unknown OTFc.@[apj][k]:
 484           %  show OTFc.@[apj][k];
 485           %fi;
 486           fi;
 487         endfor;
 488         if ProduitScalaire(Oeil-Fc[apj][1],Normal(Fc[apj][1],Fc[apj][2],Fc[apj][3]))>=0:
 489           Vue[apj]:=true;
 490           apj:=apj+1;
 491         else:
 492           if Creux=true:
 493             Vue[apj]:=false;
 494             apj:=apj+1;
 495           fi;
 496         fi;
 497       fi;
 498     fi;
 499     exitif s_=EOF;
 500   endfor;
 501   apj:=apj-1;
 502   closefrom nomfichier;
 503   DessineObjetNew
 504 enddef;
 505
 506 nb:=8;
 507
 508 %tube
 509 vardef Tube(expr Fn,dp,ds,rayon,tmin,nbp,pas)=%f,f',f'',rayon du tube,paramètre départ,nb points, pas
 510   save _tube;
 511   picture _tube;
 512   scantokens("vardef F(expr t)="&Fn&" enddef;");
 513   scantokens("vardef Fp(expr t)="&dp&" enddef;");
 514   scantokens("vardef Fd(expr t)="&ds&" enddef;");
 515   color G[][];
 516 %nb point sur le cercle
 517   NB:=nbp;%nb de pas
 518   for l=0 upto NB:
 519     for k=0 upto nb:
 520       G[l][k]=F(tmin+l*pas)+rayon*(cosd(k*(360/nb))*VN(tmin+l*pas)+sind(k*(360/nb))*VBN(tmin+l*pas));
 521     endfor;
 522   endfor;
 523   apj:=0;
 524   for l=0 upto (NB-1):
 525     for k=0 upto nb-1:
 526       apj:=apj+1;
 527       cpt[apj]:=apj;
 528       Fc[apj][1]:=G[l][k];
 529       Fc[apj][2]:=G[l][k+1];
 530       Fc[apj][3]:=G[l+1][k+1];
 531       Fc[apj][4]:=G[l+1][k];
 532       Fc[apj].iso:=(Fc[apj][1]+Fc[apj][2]+Fc[apj][3]+Fc[apj][4])/4;
 533       ALT[apj]:=-Zpart(GCoord(Fc[apj][1]));
 534     endfor;
 535   endfor;
 536   QS(1,apj);
 537   _tube=image(
 538     for k=1 upto apj:
 539       fill for l=1 upto 4:
 540         Projette(Fc[cpt[k]][l])--
 541       endfor
 542       cycle withcolor if arcenciel: lumin(cpt[k])*Hsvtorgb((floor((cpt[k]/apj)*360),satu,lum))
 543         else: lumin(cpt[k])*outcolor fi;
 544       draw for l=1 upto 4:
 545         Projette(Fc[cpt[k]][l])--
 546       endfor
 547       cycle withpen pencircle scaled0.25bp;
 548     endfor;
 549     );
 550   _tube
 551 enddef;
 552
 553 %Tube new :)
 554 %pour les tubes new :)
 555 vardef T(expr t)=Fp(t)
 556 enddef;
 557
 558 vardef VTn(expr t)=if Norm(T(t))=0:
 559     T(t)
 560   else:
 561     T(t)/Norm(T(t))
 562   fi
 563 enddef;
 564
 565 vardef VNn(expr t)=
 566   save _VNBis,__VN;
 567   color _VNBis,__VN;
 568   __VN=T(t) Vectprod ((T(t+nn)-T(t-nn))/2);
 569   if __VN=(0,0,0):
 570     _VNBis=(0,0,1);
 571   else:
 572     __VN:=__VN/Norm(__VN);
 573     if ProduitScalaire(VNbisprec[t-nn],__VN)>0:
 574       _VNBis=__VN/Norm(__VN)
 575     else:
 576       _VNBis=-(__VN/Norm(__VN))
 577     fi;
 578   fi;
 579   VNbisprec[t]=_VNBis;
 580   _VNBis
 581 enddef;
 582
 583 vardef VBNn(expr t)=VTn(t) Vectprod VNn(t)
 584 enddef;
 585
 586 vardef Tuben(expr Fn,dp,rayon,tmin,nbp,pas)=%f,f',f'',rayon du tube,paramètre départ,nb points, pas
 587   save _tuben;
 588   picture _tuben;
 589   scantokens("vardef F(expr t)="&Fn&" enddef;");
 590   scantokens("vardef Fp(expr t)="&dp&" enddef;");
 591   color G[][];
 592 %nb point sur le cercle
 593   NB:=nbp;%nb de pas
 594   nn:=pas;
 595   %pour gérer le "flip" aux points d'inflexion
 596   color VNbisprec[];
 597   VNbisprec[tmin-nn]=T(tmin-nn) Vectprod ((T(tmin)-T(tmin-2*nn))/2);
 598   %
 599   ang:=360/nb;
 600   for l=0 upto NB:
 601     for k=0 upto nb:
 602       G[l][k]=F(tmin+l*pas)+rayon*(cosd(k*ang)*VNn(tmin+l*pas)+sind(k*ang)*VBNn(tmin+l*pas));
 603     endfor;
 604   endfor;
 605   apj:=0;
 606   for l=0 upto (NB-1):
 607     for k=0 upto nb-1:
 608       apj:=apj+1;
 609       cpt[apj]:=apj;
 610       Fc[apj][1]:=G[l][k];
 611       Fc[apj][2]:=G[l][k+1];
 612       Fc[apj][3]:=G[l+1][k+1];
 613       Fc[apj][4]:=G[l+1][k];
 614       Fc[apj].iso:=(Fc[apj][1]+Fc[apj][2]+Fc[apj][3]+Fc[apj][4])/4;
 615       ALT[apj]:=-Zpart(GCoord(Fc[apj].iso));
 616     endfor;
 617   endfor;
 618   QS(1,apj);
 619   _tuben=image(
 620     for k=1 upto apj:
 621       fill for l=1 upto 4:
 622         Projette(Fc[cpt[k]][l])--
 623       endfor
 624       cycle withcolor if arcenciel: lumin(cpt[k])*Hsvtorgb((floor((cpt[k]/apj)*360),satu,lum))
 625         else: lumin(cpt[k])*outcolor fi;
 626       draw for l=1 upto 4:
 627         Projette(Fc[cpt[k]][l])--
 628       endfor
 629       cycle withpen pencircle scaled0.25bp;
 630     endfor;
 631     );
 632   _tuben
 633 enddef;
 634
 635
 636 vardef Fonction(expr fn,tmin,tmax,pas)=%fonction
 637   scantokens("vardef F(expr t)="&fn&" enddef;");
 638   save _f;
 639   path _f;
 640   _f=Projette(F(tmin))
 641   for k=tmin+pas step pas until tmax:
 642     --Projette(F(k))
 643   endfor;
 644   _f
 645 enddef;
 646
 647 color outcolor,incolor,outcolorbis;
 648
 649 boolean Spar;
 650 Spar:=false;
 651
 652 vardef Sparam(expr fn,umin,umax,upas,vmin,vmax,vpas)=%fonction
 653   Spar:=true;
 654   save _Sparam;
 655   scantokens("vardef Famille(expr u,v)="&fn&" enddef;");
 656   apj:=0;
 657   picture _Sparam;
 658   %On crée les facettes et on calcule la profondeur en Z.
 659   for k=umin step upas until umax:
 660     for l=vmin step vpas until vmax:
 661       apj:=apj+1;
 662       cpt[apj]:=apj;
 663       Fc[apj][1]:=Image(Famille(k+upas,l));
 664       Fc[apj][2]:=Image(Famille(k,l));
 665       Fc[apj][3]:=Image(Famille(k,l+vpas));
 666       Fc[apj][4]:=Image(Famille(k+upas,l+vpas));
 667       Fc[apj].iso:=(Fc[apj][1]+Fc[apj][3])/2;%(Fc[apj][1]+Fc[apj][2]+Fc[apj][3]+Fc[apj][4])/4;
 668       ALT[apj]:=-Zpart(GCoord(Fc[apj].iso));
 669       if ProduitScalaire(Oeil-Fc[apj].iso,invnormale*Normal(Fc[apj].iso,Fc[apj][1],Fc[apj][2]))>=0:
 670         Vue[apj]:=true
 671       else:
 672         Vue[apj]:=false
 673       fi;
 674     endfor;
 675   endfor;
 676   %On range les faces par un QS en fonction de leur profondeur
 677   QS(1,apj);
 678   %On affiche toutes les faces par ordre décroissant de profondeur.
 679   _Sparam=image(
 680     for k=1 upto apj:
 681       fill for l=1 upto 4:
 682         Projette(Fc[cpt[k]][l])--
 683       endfor
 684       cycle withcolor if Vue[cpt[k]]:
 685         if arcenciel: lumin(cpt[k])*Hsvtorgb((floor((cpt[k]/apj)*360),satu,lum))
 686         else: lumin(cpt[k])*outcolor fi
 687       else:lumin(cpt[k])*incolor fi;
 688       draw for l=1 upto 4:
 689         Projette(Fc[cpt[k]][l])--
 690       endfor
 691       cycle;
 692     endfor;
 693     );
 694   Spar:=false;
 695   _Sparam
 696 enddef;
 697
 698 vardef Revolution(expr fn,umin,umax,upas,vmin,vmax,vpas)=
 699   Sparam("(xpart(point(u) of "&fn&")*cos(v),xpart(point(u) of "&fn&")*sin(v),ypart(point(u) of "&fn&"))",umin,umax,upas,vmin,vmax,vpas)
 700 enddef;
 701
 702 boolean traits;%sur une idée d'Herbert Voss à propos de pst-solides :)
 703 %sert à désactiver le tracer des traits
 704 %15/07/08:pour l'instant implanter uniquement pour les surfaces en z
 705 traits=true;
 706
 707 boolean arcenciel;%pour essayer d'obtenir des dégradés tels que pst-solides :)
 708 arcenciel=false;
 709
 710 boolean surfz;
 711 surfz:=false;
 712
 713 vardef SurfZ(expr fn,xmin,xmax,ymin,ymax,nblignes,nbpoints)=
 714   surfz:=true;
 715   save _SurfZ;
 716   picture _SurfZ;
 717   %boolean Vue[];
 718   color alt[];
 719   scantokens("vardef Fz(expr X,Y)="&fn&" enddef;");
 720   apj:=0;sign:=0;
 721   IncX:=(xmax-xmin)/nbpoints;
 722   IncY:=(ymax-ymin)/nblignes;
 723   color Yc[][],Xc[][],Fc[][];
 724   for ligne=0 upto nblignes:
 725     y:=ymax-ligne*IncY;
 726     x:=xmin;
 727     Yc[ligne][0]=(x,y,Fz(x,y));
 728     for k=1 upto nbpoints:
 729       Yc[ligne][k]=((xmin+k*IncX,y,Fz(xmin+k*IncX,y)));
 730     endfor;
 731   endfor;
 732   %for ligne=0 step 3 until nbpoints:
 733   %  x:=xmax-ligne*IncX;
 734   %  for k=0 step 3 until nblignes:
 735   %    Xc[ligne div 3][k div 3]=(x,ymin+k*IncY,Fz(x,ymin+k*IncY));
 736   %  endfor;
 737   %endfor;
 738   for k=0 upto (nblignes-1):
 739     for l=0 step 3 until (nbpoints-3):
 740       apj:=apj+1;
 741       cpt[apj]:=apj;
 742       Fc[apj][1]:=Yc[k][l];
 743       Fc[apj][2]:=Yc[k][l+3];
 744       Fc[apj][3]:=Yc[k+1][l+3];
 745       Fc[apj][4]:=Yc[k+1][l];
 746       Fc[apj].iso:=(Fc[apj][1]+Fc[apj][2]+Fc[apj][3]+Fc[apj][4])/4;
 747       ALT[apj]:=-Zpart(GCoord(Fc[apj].iso));
 748       if ProduitScalaire(Oeil-Fc[apj].iso,Normal(Fc[apj].iso,Fc[apj][2],Fc[apj][1]))>=0:
 749         Vue[apj]:=true;
 750       else:
 751         Vue[apj]:=false
 752       fi;
 753     endfor;
 754   endfor;
 755   %On range les faces par un QS en fonction de leur profondeur
 756   QS(1,apj);
 757   %On affiche toutes les faces par ordre décroissant de profondeur.
 758   _SurfZ=image(
 759     pickup pencircle scaled 0.25bp;
 760     for k=1 upto apj:
 761       fill for l=1 upto 4:
 762         Projette(Fc[cpt[k]][l])--
 763       endfor
 764       cycle withcolor if Vue[cpt[k]]:
 765         if arcenciel: lumin(cpt[k])*Hsvtorgb((floor((cpt[k]/apj)*360),satu,lum))
 766         else:lumin(cpt[k])*outcolor fi
 767       else: lumin(cpt[k])*incolor fi;
 768       if traits=true:
 769         draw for l=1 upto 4:
 770           Projette(Fc[cpt[k]][l])--
 771         endfor
 772         cycle;
 773       fi;
 774     endfor;
 775     );
 776   surfz:=false;
 777   _SurfZ
 778 enddef;
 779
 780 %Objet prédéfinis
 781 vardef ObjetCube(expr ar)=
 782   NbS:=8;
 783   Sommet1:=(ar,0,0);
 784   Sommet2:=(ar,ar,0);
 785   Sommet3:=(0,ar,0);
 786   Sommet4:=(0,0,0);
 787   Sommet5:=(0,0,ar);
 788   Sommet6:=(ar,0,ar);
 789   Sommet7:=(ar,ar,ar);
 790   Sommet8:=(0,ar,ar);
 791 %%Faces
 792   NF:=6;
 793   ns[1][0]:=4;Fc[1][1]:=Sommet1;Fc[1][2]:=Sommet2;Fc[1][3]:=Sommet3;Fc[1][4]:=Sommet4;
 794   ns[2][0]:=4;Fc[2][1]:=Sommet4;Fc[2][2]:=Sommet3;Fc[2][3]:=Sommet8;Fc[2][4]:=Sommet5;
 795   ns[3][0]:=4;Fc[3][1]:=Sommet1;Fc[3][2]:=Sommet4;Fc[3][3]:=Sommet5;Fc[3][4]:=Sommet6;
 796   ns[4][0]:=4;Fc[4][1]:=Sommet5;Fc[4][2]:=Sommet8;Fc[4][3]:=Sommet7;Fc[4][4]:=Sommet6;
 797   ns[5][0]:=4;Fc[5][1]:=Sommet2;Fc[5][2]:=Sommet7;Fc[5][3]:=Sommet8;Fc[5][4]:=Sommet3;
 798   ns[6][0]:=4;Fc[6][1]:=Sommet1;Fc[6][2]:=Sommet6;Fc[6][3]:=Sommet7;Fc[6][4]:=Sommet2;
 799   %
 800   %Détermination des faces
 801   apj:=0;color cc,dd;
 802   for nf=1 upto NF:
 803     apj:=apj+1;
 804     dd:=Oeil-Fc[nf][1];
 805     cc:=invnormale*Normal(Fc[nf][1],Fc[nf][2],Fc[nf][3]);
 806     if (ProduitScalaire(dd,cc)>=0):
 807       Vue[apj]=true
 808     else:
 809       Vue[apj]=false;
 810     fi;
 811   endfor;
 812 enddef;
 813
 814 %coloriage et lumière
 815 vardef Hsvtorgb(expr CC)=%CC couleur donnée en hsv d'après http://en.wikipedia.org/wiki/HSL_color_space
 816   save $;
 817   color $;
 818   SSw:=floor(Xpart(CC)/60);
 819   SSh:=SSw mod 6;
 820   SSf:=(Xpart(CC)/60)-floor(SSw);
 821   SSs:=Ypart((CC));
 822   SSv:=Zpart((CC));
 823   SSp:=SSv*(1-SSs);
 824   SSq:=SSv*(1-SSf*SSs);
 825   SSt:=SSv*(1-(1-SSf)*SSs);
 826   if SSh=0: $=(SSv,SSt,SSp) elseif SSh=1:$=(SSq,SSv,SSp) elseif SSh=2:$=(SSp,SSv,SSt) elseif SSh=3:$=(SSp,SSq,SSv) elseif SSh=4:$=(SSt,SSp,SSv) elseif SSh=5:$=(SSv,SSp,SSq) fi;
 827   $
 828 enddef;
 829
 830 marron=Hsvtorgb((60,1,0.3));
 831
 832 intensite:=2;
 833 boolean eclairage;
 834 eclairage:=true;
 835
 836 color Lumiere;
 837
 838 invnormalelum:=1;
 839
 840 vardef lumin(expr nbt)=
 841   save $;
 842   numeric $;
 843   if eclairage=false:
 844     $=1;
 845   else:
 846     color uu,vv;
 847     uu=Lumiere-Fc[nbt].iso;
 848     uu:=uu/Norm(uu);
 849     vv=invnormalelum*Normal(Fc[nbt].iso,Fc[nbt][1],Fc[nbt][2]);
 850     if Norm(vv)<>0:
 851       vv:=vv/Norm(vv)
 852     fi;
 853     if (surfz) or (Spar) or (OFF) or (OBJ):
 854       $=intensite*abs(ProduitScalaire(vv,uu))
 855     else:
 856       $=intensite*(ProduitScalaire(vv,uu));
 857     fi;
 858     if $>1:
 859       $:=1;
 860     fi;
 861   fi;
 862   $
 863 enddef;
 864
 865
 866 %%sucre
 867
 868 %%Transparence fait par Anthony Phan
 869 picture alphapict_; alphapict_=nullpicture;
 870 color fillcolor; fillcolor=gris;
 871 fgalpha := 0.5; % usual alpha parameter
 872 bgalpha:= 1; % alpha parameter with respect to the background
 873
 874 vardef transparence expr c =
 875   alphapict_ := nullpicture;
 876   alphafill_(currentpicture, c);
 877   addto currentpicture also alphapict_;
 878 enddef;
 879
 880 def alphafill_(expr p, c) =
 881   begingroup
 882     save p_, xmax_, xmin_, ymax_, ymin_; picture p_;
 883     p_ = nullpicture;
 884     (xmin_, ymin_) = llcorner c; (xmax_, ymax_) = urcorner c;
 885     addto p_ contour c withcolor bgalpha[background, fillcolor];
 886     for p__ within p:
 887       numeric xmin__, xmax__, ymin__, ymax__;
 888       (xmin__, ymin__) = llcorner p__; (xmax__, ymax__) = urcorner p__;
 889       if (xmax__<= xmin_) or (xmin__ >= xmax_):
 890       else:
 891         if (ymax__<= ymin_) or (ymin__ >= ymax_):
 892         else:
 893           if (not clipped p__) and (not bounded p__):
 894             addto p_ also p__ withcolor
 895             fgalpha[(redpart p__, greenpart p__, bluepart p__),
 896             fillcolor];
 897           else:
 898             begingroup save alphapict_;
 899               picture alphapict_; alphapict_ = nullpicture;
 900               alphafill_(p__, pathpart p__);
 901               addto p_ also alphapict_;
 902             endgroup;
 903           fi
 904         fi
 905       fi
 906     endfor
 907     clip p_ to c;
 908     addto alphapict_ also p_;
 909   endgroup;
 910 enddef;
 911
 912 endinput;