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 vardef TraceGrille(expr NB)=
 156   color ppt[];
 157   for k=0 upto NB:
 158     ppt[k]:=(k,0,0);
 159     draw Projette(ppt[k]+(0,0,NB))--Projette(ppt[k])--Projette(ppt[k]+(0,NB,0));
 160     ppt[k]:=(0,k,0);
 161     draw Projette(ppt[k]+(NB,0,0))--Projette(ppt[k])--Projette(ppt[k]+(0,0,NB));
 162     ppt[k]:=(0,0,k);
 163     draw Projette(ppt[k]+(NB,0,0))--Projette(ppt[k])--Projette(ppt[k]+(0,NB,0));
 164   endfor;
 165 enddef;
 166
 167 primarydef u Vectprod v =
 168   (Ypart u * Zpart v - Zpart u * Ypart v,
 169     Zpart u * Xpart v - Xpart u * Zpart v,
 170     Xpart u * Ypart v - Ypart u * Xpart v)
 171 enddef;
 172
 173 vardef Normal(expr vecun,vecde,vectr)=
 174   save aa;
 175   color aa;
 176   P1:=redpart(vecde-vecun);
 177   P2:=greenpart(vecde-vecun);
 178   P3:=bluepart(vecde-vecun);
 179   Q1:=redpart(vectr-vecun);
 180   Q2:=greenpart(vectr-vecun);
 181   Q3:=bluepart(vectr-vecun);
 182   aa=(P2*Q3-Q2*P3,P3*Q1-Q3*P1,P1*Q2-Q1*P2);
 183   aa
 184 enddef;
 185
 186 vardef ProduitScalaire(expr wec,mor)=
 187   %Mexp(Mlog redpart(wec) Mmul Mlog redpart(mor))+Mexp(Mlog greenpart(wec) Mmul Mlog greenpart(mor))+Mexp(Mlog bluepart(wec) Mmul Mlog bluepart(mor))
 188   Xpart(wec)*Xpart(mor)+Ypart(wec)*Ypart(mor)+Zpart(wec)*Zpart(mor)
 189 enddef;
 190 %pour les rotations et translations
 191
 192 vardef RotX(expr ptx)=
 193   (Xpart(ptx),Ypart(ptx)*cosd(angx)-Zpart(ptx)*sind(angx),Ypart(ptx)*sind(angx)+Zpart(ptx)*cosd(angx))
 194 enddef;
 195
 196 vardef RotY(expr ptx)=
 197   (Xpart(ptx)*cosd(angy)+Zpart(ptx)*sind(angy),Ypart(ptx),-Xpart(ptx)*sind(angy)+Zpart(ptx)*cosd(angy))
 198 enddef;
 199
 200 vardef RotZ(expr ptx)=
 201   (Xpart(ptx)*cosd(angz)-Ypart(ptx)*sind(angz),Xpart(ptx)*sind(angz)+Ypart(ptx)*cosd(angz),Zpart(ptx))
 202 enddef;
 203
 204 vardef RotXYZ(expr ptx)=
 205   RotZ(RotY(RotX(ptx)))
 206 enddef;
 207
 208 angx:=0;
 209 angy:=0;
 210 angz:=0;
 211 color TR;
 212 TR:=(0,0,0);
 213 %pour les tubes
 214 vardef VT(expr t)=Fp(t)/Norm(Fp(t))
 215 enddef;
 216
 217 vardef VN(expr t)=Fd(t)/Norm(Fd(t))
 218 enddef;
 219
 220 vardef VBN(expr t)=VT(t) Vectprod VN(t)
 221 enddef;
 222
 223 path feuillet;
 224 numeric _tfig,_nfig;
 225 pair coinbg,coinbd,coinhd,coinhg;
 226 _nfig:=0;
 227
 228 def feuille(expr xa,ya,xb,yb) =
 229   feuillet := (xa,ya)--(xa,yb)--(xb,yb)--(xb,ya)--cycle;
 230   coinbg := (xa,ya);
 231   coinbd := (xb,ya);
 232   coinhd := (xb,yb);
 233   coinhg := (xa,yb);
 234   %modifié le 29.09.04
 235   z.so=(xpart(coinbg/1cm),ypart(coinbg/1cm));
 236   z.ne=(xpart(coinhd/1cm),ypart(coinhd/1cm));
 237   %fin modification
 238   extra_endfig := "clip currentpicture to feuillet;" & extra_endfig;
 239 enddef;
 240
 241 def figureespace(expr xa,ya,xb,yb) =
 242   _nfig:=_nfig+1;
 243   beginfig(_nfig);
 244     feuille(xa,ya,xb,yb);
 245     _tfig:= if (xb-xa)>(yb-ya): xb-xa else: yb-ya fi;
 246         _tfig:=2*_tfig;
 247 enddef;
 248
 249 def finespace=
 250 endfig;
 251 enddef;
 252
 253 def QS(expr ndeb,nfin)=
 254   begingroup
 255     save v,m,x;
 256     if ndeb<nfin:
 257       v:=ALT[cpt[ndeb]];
 258       m:=ndeb;
 259       for i=(ndeb+1) upto nfin:
 260         if ALT[cpt[i]]<v:
 261           m:=m+1;
 262           x:=cpt[m];cpt[m]:=cpt[i];cpt[i]:=x;
 263         fi
 264       endfor;
 265       x:=cpt[m];cpt[m]:=cpt[ndeb];cpt[ndeb]:=x;
 266       QS(ndeb,m-1);
 267       QS(m+1,nfin);
 268     fi
 269   endgroup
 270 enddef;
 271
 272 boolean Pointilles;
 273 Pointilles:=false;
 274
 275 boolean Vue[];
 276 color Fc[].iso;
 277 boolean Creux;
 278 Creux:=false;
 279
 280 vardef DessineObjetNew=
 281   numeric ALT[];
 282   %on détermine les zmax dans le repère spatial de l'écran
 283   for k=1 upto apj:
 284     Fc[k].iso:=(0,0,0) for l=1 upto ns[k][0]:
 285       +Fc[k][l]
 286     endfor;
 287     Fc[k].iso:=Fc[k].iso/ns[k][0];
 288   endfor;
 289   for k=1 upto apj:
 290     %ALT[k]=-Zpart(GCoord(Fc[k].iso));
 291     zmax:=infinity;
 292     for l=1 upto ns[k][0]:
 293       t:=-Zpart(GCoord(Fc[k][l]));
 294       if t<zmax:
 295         zmax:=t;
 296       fi;
 297     endfor;
 298     ALT[k]=zmax;
 299     cpt[k]:=k;
 300   endfor;
 301   %On trie suivant les zmax
 302   QS(1,apj);
 303   %On dessine
 304   if Pointilles:
 305     for k=1 upto apj:
 306       draw for l=1 upto ns[cpt[k]][0]:
 307         Projette(Fc[cpt[k]][l])--
 308       endfor
 309       cycle if Vue[cpt[k]]=false: dashed evenly fi;
 310     endfor;
 311   else:
 312     for k=1 upto apj:
 313       if Creux=false:
 314         if Vue[cpt[k]]:
 315           fill for l=1 upto ns[cpt[k]][0]:
 316             Projette(Fc[cpt[k]][l])--
 317           endfor
 318           cycle withcolor
 319           if arcenciel:
 320             lumin(cpt[k])*Hsvtorgb(((cpt[k]/apj)*360,satu,lum))
 321           else:
 322             lumin(cpt[k])*outcolor
 323           fi;
 324           if traits:
 325             draw for l=1 upto ns[cpt[k]][0]:
 326               Projette(Fc[cpt[k]][l])--
 327             endfor
 328             cycle;
 329           fi;
 330         fi;
 331       else:
 332         fill for l=1 upto ns[cpt[k]][0]:
 333           Projette(Fc[cpt[k]][l])--
 334         endfor
 335         cycle withcolor
 336         if Vue[cpt[k]]:
 337           if arcenciel:
 338             lumin(cpt[k])*Hsvtorgb(((cpt[k]/apj)*360,satu,lum))
 339           else:
 340             lumin(cpt[k])*outcolor
 341           fi
 342         else:
 343           lumin(cpt[k])*incolor
 344         fi;
 345         if traits:
 346           draw for l=1 upto ns[cpt[k]][0]:
 347             Projette(Fc[cpt[k]][l])--
 348           endfor
 349           cycle;
 350         fi;
 351       fi;
 352     endfor;
 353   fi;
 354 enddef;
 355
 356 color Fc[][];%représente les sommets des différentes faces;
 357
 358 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
 359 boolean OFF,OBJ;%pour utiliser la lumière "correctement"
 360 OFF:=false;
 361 OBJ:=false;
 362
 363 vardef LectureOFF(expr nomsolide)=
 364   OFF:=true;
 365   %Détermination du nombre de sommets et de faces.
 366   string s_;
 367   s_=readfrom nomsolide;
 368   string ss[];
 369   if s_<>EOF:
 370     ss1 := loptok s_;
 371     t_ := if ss1="%": 0 else: 1 fi;
 372     forever:
 373       ss[incr t_] := loptok s_;
 374       exitif ss[t_]="";
 375     endfor
 376   else: false
 377   fi;
 378   NbS:=round(Mexp Mlog_str ss1);
 379   NF:=round(Mexp Mlog_str ss2);
 380   message("Il y a "&decimal(NbS)&" sommets.");
 381   message("Il y a "&decimal(NF)&" faces au total.");
 382   %Détermination des coordonnées des sommets
 383   message("Création des sommets.");
 384   s_:=readfrom nomsolide;
 385   if debut=0:
 386     for k=0 upto NbS-1:
 387       s_:=readfrom nomsolide;
 388       if s_<>EOF:
 389         ss1 := loptok s_;
 390         n_ := if ss1="%": 0 else: 1 fi;
 391         forever:
 392           ss[incr n_] := loptok s_;
 393           exitif ss[n_]="";
 394         endfor
 395       else: false
 396       fi;
 397       Sommet[k]:=(Mexp ((Mlog_str ss1) Mdiv (Mlog echelle)),Mexp ((Mlog_str ss3) Mdiv (Mlog echelle)),Mexp ((Mlog_str ss2) Mdiv (Mlog echelle)));
 398       %Sommet[k]:=(Mexp (Mlog_str ss1),Mexp (Mlog_str ss3),Mexp (Mlog_str ss2))/echelle;
 399       %Sommet[k]:=(Mexp (Mlog_str ss1)/echelle,Mexp (Mlog_str ss3)/echelle,Mexp (Mlog_str ss2)/echelle);
 400     endfor;
 401   else:
 402     for k=1 upto NbS:
 403       s_:=readfrom nomsolide;
 404       if s_<>EOF:
 405         ss1 := loptok s_;
 406         n_ := if ss1="%": 0 else: 1 fi;
 407         forever:
 408           ss[incr n_] := loptok s_;
 409           exitif ss[n_]="";
 410         endfor
 411       else: false
 412       fi;
 413       Sommet[k]:=(Mexp ((Mlog_str ss1) Mdiv (Mlog echelle)),Mexp ((Mlog_str ss3) Mdiv (Mlog echelle)),Mexp ((Mlog_str ss2) Mdiv (Mlog echelle)));
 414       %Sommet[k]:=(Mexp (Mlog_str ss1),Mexp (Mlog_str ss3),Mexp (Mlog_str ss2))/echelle;
 415       %Sommet[k]:=(Mexp (Mlog_str ss1)/echelle,Mexp (Mlog_str ss3)/echelle,Mexp (Mlog_str ss2)/echelle);
 416     endfor;
 417   fi;
 418   message("Création des faces.");
 419   %Détermination des faces
 420   apj:=0;color cc,dd;
 421   nbfvues:=0;
 422   for nf=-4000 upto (-4000+NF)-1:
 423     s_:=readfrom nomsolide;
 424      if s_<>EOF:
 425       ss1 := loptok s_;
 426       n_ := if ss1="%": 0 else: 1 fi;
 427       forever:
 428         ss[incr n_] := loptok s_;
 429         exitif ss[n_]="";
 430       endfor
 431     else: false
 432     fi;
 433     apj:=apj+1;
 434     ns[apj][0]:=Mexp Mlog_str ss1;%pour savoir le nb de sommets par face
 435     for nl=1 upto ns[apj][0]:
 436       Fc[apj][nl]:=Sommet[round(Mexp Mlog_str ss[nl+1])];
 437     endfor;
 438     dd:=Oeil-Fc[apj][1];
 439     cc:=invnormale*Normal(Fc[apj][1],Fc[apj][2],Fc[apj][3]);
 440     if (ProduitScalaire(dd,cc)>=0):
 441       Vue[apj]:=true;
 442       nbfvues:=nbfvues+1;
 443     else:
 444       if Creux=true:
 445         Vue[apj]:=false;
 446       else:
 447         apj:=apj-1;
 448       fi;
 449     fi;
 450   endfor;
 451   message("Faces vues déterminees : il y en a "&decimal(nbfvues)&".");
 452   closefrom nomsolide;
 453   DessineObjetNew;
 454 enddef;
 455
 456 vardef LectureOBJ(expr nomfichier)=
 457   OBJ:=true;
 458   string s_;
 459   string ss[];
 460   nbss:=1;
 461   apj:=0;
 462   forever:
 463     s_:=readfrom nomfichier;
 464     if s_<>EOF:
 465       ss0 := loptok s_;
 466       if ss0="v":
 467         n_:=0;
 468         forever:
 469           ss[incr n_] := loptok s_;
 470           exitif ss[n_]="";
 471         endfor
 472         Sommet[nbss]:=(Mexp((Mlog_str ss1) Mdiv (Mlog echelle)),Mexp ((Mlog_str ss3) Mdiv (Mlog echelle)),Mexp ((Mlog_str ss2) Mdiv (Mlog echelle)));
 473         nbss:=incr nbss;
 474       elseif ss0="f":
 475         n_:=0;
 476         forever:
 477           ss[incr n_] := loptok s_;
 478           exitif ss[n_]="";
 479         endfor;
 480         ns[apj][0]:=n_-1;
 481         for k=1 upto ns[apj][0]:
 482           if invnormale=1:
 483             Fc[apj][ns[apj][0]-k+1] := Sommet[round(Mexp(Mlog_str ss[k]))]
 484             %if unknown OTFc.@[apj][OTFc.@[apj].nb-k+1]:
 485           %  show OTFc.@[apj][OTFc.@[apj].nb-k+1];
 486           %fi;
 487           else:
 488             Fc[apj][k] := Sommet[round(Mexp(Mlog_str ss[k]))]
 489             %if unknown OTFc.@[apj][k]:
 490           %  show OTFc.@[apj][k];
 491           %fi;
 492           fi;
 493         endfor;
 494         if ProduitScalaire(Oeil-Fc[apj][1],Normal(Fc[apj][1],Fc[apj][2],Fc[apj][3]))>=0:
 495           Vue[apj]:=true;
 496           apj:=apj+1;
 497         else:
 498           if Creux=true:
 499             Vue[apj]:=false;
 500             apj:=apj+1;
 501           fi;
 502         fi;
 503       fi;
 504     fi;
 505     exitif s_=EOF;
 506   endfor;
 507   apj:=apj-1;
 508   closefrom nomfichier;
 509   DessineObjetNew
 510 enddef;
 511
 512 nb:=8;
 513
 514 %tube
 515 vardef Tube(expr Fn,dp,ds,rayon,tmin,nbp,pas)=%f,f',f'',rayon du tube,paramètre départ,nb points, pas
 516   save _tube;
 517   picture _tube;
 518   scantokens("vardef F(expr t)="&Fn&" enddef;");
 519   scantokens("vardef Fp(expr t)="&dp&" enddef;");
 520   scantokens("vardef Fd(expr t)="&ds&" enddef;");
 521   color G[][];
 522 %nb point sur le cercle
 523   NB:=nbp;%nb de pas
 524   for l=0 upto NB:
 525     for k=0 upto nb:
 526       G[l][k]=F(tmin+l*pas)+rayon*(cosd(k*(360/nb))*VN(tmin+l*pas)+sind(k*(360/nb))*VBN(tmin+l*pas));
 527     endfor;
 528   endfor;
 529   apj:=0;
 530   for l=0 upto (NB-1):
 531     for k=0 upto nb-1:
 532       apj:=apj+1;
 533       cpt[apj]:=apj;
 534       Fc[apj][1]:=G[l][k];
 535       Fc[apj][2]:=G[l][k+1];
 536       Fc[apj][3]:=G[l+1][k+1];
 537       Fc[apj][4]:=G[l+1][k];
 538       Fc[apj].iso:=(Fc[apj][1]+Fc[apj][2]+Fc[apj][3]+Fc[apj][4])/4;
 539       ALT[apj]:=-Zpart(GCoord(Fc[apj][1]));
 540     endfor;
 541   endfor;
 542   QS(1,apj);
 543   _tube=image(
 544     for k=1 upto apj:
 545       fill for l=1 upto 4:
 546         Projette(Fc[cpt[k]][l])--
 547       endfor
 548       cycle withcolor if arcenciel: lumin(cpt[k])*Hsvtorgb((floor((cpt[k]/apj)*360),satu,lum))
 549         else: lumin(cpt[k])*outcolor fi;
 550       draw for l=1 upto 4:
 551         Projette(Fc[cpt[k]][l])--
 552       endfor
 553       cycle withpen pencircle scaled0.25bp;
 554     endfor;
 555     );
 556   _tube
 557 enddef;
 558
 559 %Tube new :)
 560 %pour les tubes new :)
 561 vardef T(expr t)=Fp(t)
 562 enddef;
 563
 564 vardef VTn(expr t)=if Norm(T(t))=0:
 565     T(t)
 566   else:
 567     T(t)/Norm(T(t))
 568   fi
 569 enddef;
 570
 571 vardef VNn(expr t)=
 572   save _VNBis,__VN;
 573   color _VNBis,__VN;
 574   __VN=T(t) Vectprod ((T(t+nn)-T(t-nn))/2);
 575   if __VN=(0,0,0):
 576     _VNBis=(0,0,1);
 577   else:
 578     __VN:=__VN/Norm(__VN);
 579     if ProduitScalaire(VNbisprec[t-nn],__VN)>0:
 580       _VNBis=__VN/Norm(__VN)
 581     else:
 582       _VNBis=-(__VN/Norm(__VN))
 583     fi;
 584   fi;
 585   VNbisprec[t]=_VNBis;
 586   _VNBis
 587 enddef;
 588
 589 vardef VBNn(expr t)=VTn(t) Vectprod VNn(t)
 590 enddef;
 591
 592 vardef Tuben(expr Fn,dp,rayon,tmin,nbp,pas)=%f,f',f'',rayon du tube,paramètre départ,nb points, pas
 593   save _tuben;
 594   picture _tuben;
 595   scantokens("vardef F(expr t)="&Fn&" enddef;");
 596   scantokens("vardef Fp(expr t)="&dp&" enddef;");
 597   color G[][];
 598 %nb point sur le cercle
 599   NB:=nbp;%nb de pas
 600   nn:=pas;
 601   %pour gérer le "flip" aux points d'inflexion
 602   color VNbisprec[];
 603   VNbisprec[tmin-nn]=T(tmin-nn) Vectprod ((T(tmin)-T(tmin-2*nn))/2);
 604   %
 605   ang:=360/nb;
 606   for l=0 upto NB:
 607     for k=0 upto nb:
 608       G[l][k]=F(tmin+l*pas)+rayon*(cosd(k*ang)*VNn(tmin+l*pas)+sind(k*ang)*VBNn(tmin+l*pas));
 609     endfor;
 610   endfor;
 611   apj:=0;
 612   for l=0 upto (NB-1):
 613     for k=0 upto nb-1:
 614       apj:=apj+1;
 615       cpt[apj]:=apj;
 616       Fc[apj][1]:=G[l][k];
 617       Fc[apj][2]:=G[l][k+1];
 618       Fc[apj][3]:=G[l+1][k+1];
 619       Fc[apj][4]:=G[l+1][k];
 620       Fc[apj].iso:=(Fc[apj][1]+Fc[apj][2]+Fc[apj][3]+Fc[apj][4])/4;
 621       ALT[apj]:=-Zpart(GCoord(Fc[apj].iso));
 622     endfor;
 623   endfor;
 624   QS(1,apj);
 625   _tuben=image(
 626     for k=1 upto apj:
 627       fill for l=1 upto 4:
 628         Projette(Fc[cpt[k]][l])--
 629       endfor
 630       cycle withcolor if arcenciel: lumin(cpt[k])*Hsvtorgb((floor((cpt[k]/apj)*360),satu,lum))
 631         else: lumin(cpt[k])*outcolor fi;
 632       draw for l=1 upto 4:
 633         Projette(Fc[cpt[k]][l])--
 634       endfor
 635       cycle withpen pencircle scaled0.25bp;
 636     endfor;
 637     );
 638   _tuben
 639 enddef;
 640
 641
 642 vardef Fonction(expr fn,tmin,tmax,pas)=%fonction
 643   scantokens("vardef F(expr t)="&fn&" enddef;");
 644   save _f;
 645   path _f;
 646   _f=Projette(F(tmin))
 647   for k=tmin+pas step pas until tmax:
 648     --Projette(F(k))
 649   endfor;
 650   _f
 651 enddef;
 652
 653 color outcolor,incolor,outcolorbis;
 654
 655 boolean Spar;
 656 Spar:=false;
 657
 658 vardef Sparam(expr fn,umin,umax,upas,vmin,vmax,vpas)=%fonction
 659   Spar:=true;
 660   save _Sparam;
 661   scantokens("vardef Famille(expr u,v)="&fn&" enddef;");
 662   apj:=0;
 663   picture _Sparam;
 664   %On crée les facettes et on calcule la profondeur en Z.
 665   for k=umin step upas until umax:
 666     for l=vmin step vpas until vmax:
 667       apj:=apj+1;
 668       cpt[apj]:=apj;
 669       Fc[apj][1]:=Image(Famille(k+upas,l));
 670       Fc[apj][2]:=Image(Famille(k,l));
 671       Fc[apj][3]:=Image(Famille(k,l+vpas));
 672       Fc[apj][4]:=Image(Famille(k+upas,l+vpas));
 673       Fc[apj].iso:=(Fc[apj][1]+Fc[apj][3])/2;%(Fc[apj][1]+Fc[apj][2]+Fc[apj][3]+Fc[apj][4])/4;
 674       ALT[apj]:=-Zpart(GCoord(Fc[apj].iso));
 675       if ProduitScalaire(Oeil-Fc[apj].iso,invnormale*Normal(Fc[apj].iso,Fc[apj][1],Fc[apj][2]))>=0:
 676         Vue[apj]:=true
 677       else:
 678         Vue[apj]:=false
 679       fi;
 680     endfor;
 681   endfor;
 682   %On range les faces par un QS en fonction de leur profondeur
 683   QS(1,apj);
 684   %On affiche toutes les faces par ordre décroissant de profondeur.
 685   _Sparam=image(
 686     for k=1 upto apj:
 687       fill for l=1 upto 4:
 688         Projette(Fc[cpt[k]][l])--
 689       endfor
 690       cycle withcolor if Vue[cpt[k]]:
 691         if arcenciel: lumin(cpt[k])*Hsvtorgb((floor((cpt[k]/apj)*360),satu,lum))
 692         else: lumin(cpt[k])*outcolor fi
 693       else:lumin(cpt[k])*incolor fi;
 694       if traits=true:
 695         draw for l=1 upto 4:
 696             Projette(Fc[cpt[k]][l])--
 697         endfor
 698         cycle;
 699       else:
 700         draw for l=1 upto 4:
 701             Projette(Fc[cpt[k]][l])--
 702         endfor
 703         cycle withcolor if Vue[cpt[k]]:
 704             if arcenciel: lumin(cpt[k])*Hsvtorgb((floor((cpt[k]/apj)*360),satu,lum))
 705           else: lumin(cpt[k])*outcolor fi
 706         else:lumin(cpt[k])*incolor fi;
 707       fi;
 708     endfor;
 709     );
 710   Spar:=false;
 711   _Sparam
 712 enddef;
 713
 714 vardef Revolution(expr fn,umin,umax,upas,vmin,vmax,vpas)=
 715   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)
 716 enddef;
 717
 718 boolean traits;%sur une idée d'Herbert Voss à propos de pst-solides :)
 719 %sert à désactiver le tracer des traits
 720 %15/07/08:pour l'instant implanter uniquement pour les surfaces en z
 721 traits=true;
 722
 723 boolean arcenciel;%pour essayer d'obtenir des dégradés tels que pst-solides :)
 724 arcenciel=false;
 725
 726 boolean surfz;
 727 surfz:=false;
 728 boolean couleurz;
 729 couleurz=false;
 730 color cz[];
 731 boolean Mcir;
 732 Mcir:=false;
 733 rayd:=0;%sauf si division par 0 :(
 734 %ajout de angtotal et angd pour les cas dans les différents cadrans
 735 angtotal:=360;
 736 angd:=0;
 737
 738 vardef SurfZ(text t_)=
 739   surfz:=true;
 740   save _SurfZ;
 741   picture _SurfZ;
 742   color alt[];
 743   nbargument:=0;
 744   for p_=t_ :
 745     if string p_:
 746       scantokens("vardef Fz(expr X,Y)="&p_&" enddef;");
 747     else:
 748       nbargument:=nbargument+1;
 749       NN[nbargument]:=p_;
 750     fi;
 751   endfor;
 752   apj:=0;sign:=0;
 753   Zmax:=-infinity;
 754   Zmin:=infinity;
 755   if nbargument=6:
 756     xmin:=NN1;
 757     xmax:=NN2;
 758     ymin:=NN3;
 759     ymax:=NN4;
 760     nblignes:=NN5;
 761     nbpoints:=NN6;
 762     IncX:=(xmax-xmin)/nbpoints;
 763     IncY:=(ymax-ymin)/nblignes;
 764     color Yc[][],Xc[][],Fc[][];
 765     for ligne=0 upto nblignes:
 766       y:=ymax-ligne*IncY;
 767       x:=xmin;
 768       Yc[ligne][0]=(x,y,Fz(x,y));
 769       for k=1 upto nbpoints:
 770         Yc[ligne][k]=((xmin+k*IncX,y,Fz(xmin+k*IncX,y)));
 771       endfor;
 772     endfor;
 773     for k=0 upto (nblignes-1):
 774       for l=0 step 3 until (nbpoints-3):
 775         apj:=apj+1;
 776         cpt[apj]:=apj;
 777         Fc[apj][1]:=Yc[k][l];
 778         Fc[apj][2]:=Yc[k][l+3];
 779         Fc[apj][3]:=Yc[k+1][l+3];
 780         Fc[apj][4]:=Yc[k+1][l];
 781         Fc[apj].iso:=(Fc[apj][1]+Fc[apj][2]+Fc[apj][3]+Fc[apj][4])/4;
 782         ALT[apj]:=-Zpart(GCoord(Fc[apj].iso));
 783         if Zpart(Fc[apj].iso)>Zmax:
 784           Zmax:=Zpart(Fc[apj].iso);
 785         fi;
 786         if Zpart(Fc[apj].iso)<Zmin:
 787           Zmin:=Zpart(Fc[apj].iso);
 788         fi;
 789         if ProduitScalaire(Oeil-Fc[apj].iso,Normal(Fc[apj].iso,Fc[apj][2],Fc[apj][1]))>=0:
 790           Vue[apj]:=true;
 791         else:
 792           Vue[apj]:=false
 793         fi;
 794       endfor;
 795     endfor;
 796   else:
 797     raymax:=NN1;
 798     pray:=NN2;
 799     nblignes:=NN3;
 800       for r=rayd step pray until (raymax-pray):
 801       for y=0 step 1 until (nblignes-1):
 802         apj:=apj+1;
 803         cpt[apj]:=apj;
 804         Fc[apj][1]:=(r*cosd(angd+y*(angtotal/nblignes)),r*sind(angd+y*(angtotal/nblignes)),Fz(r*cosd(angd+y*(angtotal/nblignes)),r*sind(angd+y*(angtotal/nblignes))));
 805         Fc[apj][2]:=(r*cosd(angd+(y+1)*(angtotal/nblignes)),r*sind(angd+(y+1)*(angtotal/nblignes)),Fz(r*cosd(angd+(y+1)*(angtotal/nblignes)),r*sind(angd+(y+1)*(angtotal/nblignes))));
 806         Fc[apj][3]:=((r+pray)*cosd(angd+(y+1)*(angtotal/nblignes)),(r+pray)*sind(angd+(y+1)*(angtotal/nblignes)),Fz((r+pray)*cosd(angd+(y+1)*(angtotal/nblignes)),(r+pray)*sind(angd+(y+1)*(angtotal/nblignes))));
 807         Fc[apj][4]:=((r+pray)*cosd(angd+y*(angtotal/nblignes)),(r+pray)*sind(angd+y*(angtotal/nblignes)),Fz((r+pray)*cosd(angd+y*(angtotal/nblignes)),(r+pray)*sind(angd+y*(angtotal/nblignes))));
 808         Fc[apj].iso:=(Fc[apj][1]+Fc[apj][2]+Fc[apj][3]+Fc[apj][4])/4;
 809         if Zpart(Fc[apj].iso)>Zmax:
 810           Zmax:=Zpart(Fc[apj].iso);
 811         fi;
 812         if Zpart(Fc[apj].iso)<Zmin:
 813           Zmin:=Zpart(Fc[apj].iso);
 814         fi;
 815         ALT[apj]:=-Zpart(GCoord(Fc[apj].iso));
 816         if ProduitScalaire(Oeil-Fc[apj].iso,Normal(Fc[apj].iso,Fc[apj][2],Fc[apj][1]))>=0:
 817           Vue[apj]:=true;
 818         else:
 819           Vue[apj]:=false
 820         fi;
 821       endfor;
 822     endfor;
 823   fi;
 824   %On range les faces par un QS en fonction de leur profondeur
 825   QS(1,apj);
 826   %On affiche toutes les faces par ordre décroissant de profondeur.
 827   _SurfZ=image(
 828     pickup pencircle scaled 0.25bp;
 829     for k=1 upto apj:
 830       fill for l=1 upto 4:
 831         Projette(Fc[cpt[k]][l])--
 832       endfor
 833       cycle withcolor if couleurz:
 834           if unknown cz1:Hsvtorgb((floor(((Zpart(Fc[cpt[k]].iso)-Zmin)/(Zmax-Zmin))*360),satu,lum))
 835         else:
 836           ((Zpart(Fc[cpt[k]].iso)-Zmin)/(Zmax-Zmin))[cz2,cz1]
 837         fi;
 838       else:
 839         if Vue[cpt[k]]:
 840           if arcenciel: lumin(cpt[k])*Hsvtorgb((floor((cpt[k]/apj)*360),satu,lum))
 841           else:lumin(cpt[k])*outcolor fi
 842         else: lumin(cpt[k])*incolor fi;
 843         if traits=true:
 844           draw for l=1 upto 4:
 845               Projette(Fc[cpt[k]][l])--
 846           endfor
 847           cycle;
 848         fi;
 849       fi;
 850       if traits=true:
 851         draw for l=1 upto 4:
 852             Projette(Fc[cpt[k]][l])--
 853         endfor
 854         cycle;
 855       fi;
 856     endfor;
 857     );
 858   surfz:=false;
 859   _SurfZ
 860 enddef;
 861
 862 %Objet prédéfinis
 863 vardef ObjetCube(expr ar)=
 864   NbS:=8;
 865   Sommet1:=(ar,0,0);
 866   Sommet2:=(ar,ar,0);
 867   Sommet3:=(0,ar,0);
 868   Sommet4:=(0,0,0);
 869   Sommet5:=(0,0,ar);
 870   Sommet6:=(ar,0,ar);
 871   Sommet7:=(ar,ar,ar);
 872   Sommet8:=(0,ar,ar);
 873 %%Faces
 874   NF:=6;
 875   ns[1][0]:=4;Fc[1][1]:=Sommet1;Fc[1][2]:=Sommet2;Fc[1][3]:=Sommet3;Fc[1][4]:=Sommet4;
 876   ns[2][0]:=4;Fc[2][1]:=Sommet4;Fc[2][2]:=Sommet3;Fc[2][3]:=Sommet8;Fc[2][4]:=Sommet5;
 877   ns[3][0]:=4;Fc[3][1]:=Sommet1;Fc[3][2]:=Sommet4;Fc[3][3]:=Sommet5;Fc[3][4]:=Sommet6;
 878   ns[4][0]:=4;Fc[4][1]:=Sommet5;Fc[4][2]:=Sommet8;Fc[4][3]:=Sommet7;Fc[4][4]:=Sommet6;
 879   ns[5][0]:=4;Fc[5][1]:=Sommet2;Fc[5][2]:=Sommet7;Fc[5][3]:=Sommet8;Fc[5][4]:=Sommet3;
 880   ns[6][0]:=4;Fc[6][1]:=Sommet1;Fc[6][2]:=Sommet6;Fc[6][3]:=Sommet7;Fc[6][4]:=Sommet2;
 881   %
 882   %Détermination des faces
 883   apj:=0;color cc,dd;
 884   for nf=1 upto NF:
 885     apj:=apj+1;
 886     dd:=Oeil-Fc[nf][1];
 887     cc:=invnormale*Normal(Fc[nf][1],Fc[nf][2],Fc[nf][3]);
 888     if (ProduitScalaire(dd,cc)>=0):
 889       Vue[apj]=true
 890     else:
 891       Vue[apj]=false;
 892     fi;
 893   endfor;
 894 enddef;
 895
 896 %coloriage et lumière
 897 vardef Hsvtorgb(expr CC)=%CC couleur donnée en hsv d'après http://en.wikipedia.org/wiki/HSL_color_space
 898   save $;
 899   color $;
 900   SSw:=floor(Xpart(CC)/60);
 901   SSh:=SSw mod 6;
 902   SSf:=(Xpart(CC)/60)-floor(SSw);
 903   SSs:=Ypart((CC));
 904   SSv:=Zpart((CC));
 905   SSp:=SSv*(1-SSs);
 906   SSq:=SSv*(1-SSf*SSs);
 907   SSt:=SSv*(1-(1-SSf)*SSs);
 908   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;
 909   $
 910 enddef;
 911
 912 marron=Hsvtorgb((60,1,0.3));
 913
 914 intensite:=2;
 915 boolean eclairage;
 916 eclairage:=true;
 917
 918 color Lumiere;
 919
 920 invnormalelum:=1;
 921
 922 vardef lumin(expr nbt)=
 923   save $;
 924   numeric $;
 925   if eclairage=false:
 926     $=1;
 927   else:
 928     color uu,vv;
 929     uu=Lumiere-Fc[nbt].iso;
 930     uu:=uu/Norm(uu);
 931     vv=invnormalelum*Normal(Fc[nbt].iso,Fc[nbt][1],Fc[nbt][2]);
 932     if Norm(vv)<>0:
 933       vv:=vv/Norm(vv)
 934     fi;
 935     if (surfz) or (Spar) or (OFF) or (OBJ):
 936       $=intensite*abs(ProduitScalaire(vv,uu))
 937     else:
 938       $=intensite*(ProduitScalaire(vv,uu));
 939     fi;
 940     if $>1:
 941       $:=1;
 942     fi;
 943   fi;
 944   $
 945 enddef;
 946
 947
 948 %%sucre
 949
 950 %%Transparence fait par Anthony Phan
 951 picture alphapict_; alphapict_=nullpicture;
 952 color fillcolor; fillcolor=gris;
 953 fgalpha := 0.5; % usual alpha parameter
 954 bgalpha:= 1; % alpha parameter with respect to the background
 955
 956 vardef transparence expr c =
 957   alphapict_ := nullpicture;
 958   alphafill_(currentpicture, c);
 959   addto currentpicture also alphapict_;
 960 enddef;
 961
 962 def alphafill_(expr p, c) =
 963   begingroup
 964     save p_, xmax_, xmin_, ymax_, ymin_; picture p_;
 965     p_ = nullpicture;
 966     (xmin_, ymin_) = llcorner c; (xmax_, ymax_) = urcorner c;
 967     addto p_ contour c withcolor bgalpha[background, fillcolor];
 968     for p__ within p:
 969       numeric xmin__, xmax__, ymin__, ymax__;
 970       (xmin__, ymin__) = llcorner p__; (xmax__, ymax__) = urcorner p__;
 971       if (xmax__<= xmin_) or (xmin__ >= xmax_):
 972       else:
 973         if (ymax__<= ymin_) or (ymin__ >= ymax_):
 974         else:
 975           if (not clipped p__) and (not bounded p__):
 976             addto p_ also p__ withcolor
 977             fgalpha[(redpart p__, greenpart p__, bluepart p__),
 978             fillcolor];
 979           else:
 980             begingroup save alphapict_;
 981               picture alphapict_; alphapict_ = nullpicture;
 982               alphafill_(p__, pathpart p__);
 983               addto p_ also alphapict_;
 984             endgroup;
 985           fi
 986         fi
 987       fi
 988     endfor
 989     clip p_ to c;
 990     addto alphapict_ also p_;
 991   endgroup;
 992 enddef;
 993
 994 endinput;