Surfaces éclairées

%@AUTEUR: Jean-Michel Sarlat

% u : unité des dimensions 
% xp, yp, zp : coordonnées de la source lumineuse  
% kp : facteur d'éclairement 
% lacouleur : couleur des dalles 

numeric u,xp,yp,zp,kp; 
color lacouleur; 

u = 2cm; 

% valeurs par défaut 

xp := 3; 
yp := 3; 
zp := 5; 
kp := 30; 

lacouleur := red+green; 

% O, I, J, K : points de base du repère 

pair O,I,J,K; 
O = (0,0); 
I = (-.7,-.7) scaled u; 
J = (1,0) scaled u; 
K = (0,1) scaled u; 

% Pour visualiser la base du repère  

vardef repereXYZ = 
 drawarrow O--I; 
 drawarrow O--J; 
 drawarrow O--K; 
enddef; 

% Passage de l'espace à la représentation plane en perspective  

vardef XYZtoXY(expr x,y,z) = 
 x*I+y*J+z*K 
enddef; 

% Différentations centrales à la louche 

vardef diffx(suffix f)(expr x,y) = 
 (f(x+0.01,y)-f(x-0.01,y))/0.02 
enddef; 
vardef diffy(suffix f)(expr x,y) = 
 (f(x,y+0.01)-f(x,y-0.01))/0.02 
enddef; 

% Calcul du facteur d'éclairement pour tenir compte de 
% l'orientation de la dalle et de l'incidence 

vardef facteur(suffix f)(expr x,y,z) = 
 numeric dfx,dfy,ca,cb,cc; 
 dfx := diffx(f,x,y); 
 dfy := diffy(f,x,y); 
 ca := (zp-z)-dfy*(yp-y)-dfx*(xp-x); 
 cb := sqrt(1+dfx*dfx+dfy*dfy); 
 cc := sqrt((z-zp)*(z-zp)+(y-yp)*(y-yp)+(x-xp)*(x-xp)); 
 kp*ca/(cb*cc*cc*cc) 
enddef; 

% La « procédure » . Elle est lourde c'est un premier jet, 
% elle nécessite de gonfler sérieusement MetaPost pour 
% digérer les deux tableaux qu'elle embarque ! 

vardef surfZ(suffix f)(expr xmin,xmax,ymin,ymax,nx,ny) = 
 numeric dx,dy,xt,yt,zt,couleur[][]; 
 pair Z[][]; 
 dx := (xmax-xmin)/nx; 
 dy := (ymax-ymin)/ny; 
 for i=0 upto nx: 
   xt := xmin+i*dx; 
   for j=0 upto ny: 
     yt := ymin+j*dy; 
     zt := f(xt,yt); 
     Z[i][j] = XYZtoXY(xt,yt,zt); 
     couleur[i][j] := facteur(f,xt,yt,zt); 
   endfor 
 endfor 
 for i = 0 upto nx-1: 
   for j= 0 upto ny-1: 
     fill Z[i][j]--Z[i][j+1]--Z[i+1][j+1]--Z[i+1][j]--cycle 
          withcolor couleur[i][j]*lacouleur; 
   endfor 
 endfor 
enddef;

% Quelques valeurs 

numeric Pi; 
numeric E; 
Pi := 3.14159; 
E  := 2.71828; 

% Le sinus en radians 

vardef sin(expr x) = 
 sind(x/Pi*180) 
enddef;

beginfig(1);
  xp := 3; 
  yp := 3; 
  zp := 10; 
  kp := 100; 
  lacouleur := red; 
 
 vardef f(expr x,y) = sin(x*Pi)+sin(y*Pi) enddef; 

 surfZ(f,-2,2,-2,2,100,100); 
 
 label(btex $z=\sin x+\sin y$ etex scaled 2,(0,-7cm));
endfig;

beginfig(2);
  xp := 3; 
  yp := 3; 
  zp := 10; 
  kp := 100; 
  lacouleur := (0.85,0.65,0.13);
       
  vardef f(expr x,y) = sin((x+y)*Pi) enddef; 

  surfZ(f,-2,2,-2,2,100,100); 

  label(btex $z = \sin (x+y)$ etex scaled 2,(0,-5.5cm));
endfig;

beginfig(3);
  xp := 3; 
  yp := 3; 
  zp := 10; 
  kp := 50; 
  lacouleur := blue; 
 
  vardef f(expr x,y) = sqrt(1+4*(x*x+y*y)) enddef; 

  surfZ(f,-2,2,-2,2,70,70); 

  label(btex $4x^2+4y^2-z^2+1=0$ etex scaled 2,(0,-1cm));
endfig;

beginfig(4);
  xp := 3; 
  yp := 3; 
  zp := 10; 
  kp := 50; 
  lacouleur := blue+green; 
 
  vardef norme(expr x,y) = 
   sqrt(x*x+y*y) 
  enddef; 
 
  vardef f(expr x,y) = 
   numeric no; 
   no := norme (x,y)*4+0.001; 
   4*sin(no)/no 
  enddef; 

  surfZ(f,-2,2,-2,2,100,100); 
  
  label(btex 
	$z={\displaystyle \sin(x^2+y^2)\over \displaystyle x^2+y^2}$ etex 
 	scaled 2,(4cm,6cm));
endfig;

beginfig(5);
  xp := 3; 
  yp := -3; 
  zp := 10; 
  kp := 50; 
  lacouleur := red+green; 
 
  vardef f(expr x,y) = 
   x*x-x*y+2*x+1 
  enddef; 
 
  surfZ(f,-1,1,-2,3,29,29); 

  label(btex $z= x^2-xy+2x+1$ etex scaled 2,(3cm,6cm));
endfig;

end