mp-solid.mp

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