%% 1er essai ou presque, brut de brut
%% transcription de la macro surz proposee par JM Sarlat pour MetaPost

-4 4 setxrange
-4 4 setyrange

0 setlinewidth 

%% la source d'eclairage
/xp 3 def
/yp 3 def
/zp 7 def
/kp 8 def

/nx 100 def
/ny 100 def

%% f (x, y) = sin (x^2 + y^2)
/f {
  2 copy 
  0 eq
   {
      0 eq
         {1}
         {dup mul exch dup mul add dup Sin exch div}
      ifelse
   }
   {
      pop
      dup mul exch dup mul add dup Sin exch div
   }
  ifelse
  
} def

/vect_I {-.7 dup} def
/vect_J {1 0} def
/vect_K {0 1} def

%% syntaxe : x y z XYZtoXY --> X Y
/XYZtoXY {
3 dict begin
   /z exch def
   /y exch def
   /x exch def
   vect_I x mulv
   vect_J y mulv addv
   vect_K z mulv addv
end
} def

%% syntaxe : x y (f) diffx 
/diffx {
3 dict begin
   /la_fonction exch cvx def
   /y exch def
   /x exch def
   x .01 add y la_fonction
   x .01 sub y la_fonction sub
   .02 div
end
} def

%% syntaxe : x y (f) diffy
/diffy {
3 dict begin
   /la_fonction exch cvx def
   /y exch def
   /x exch def
   x y .01 add la_fonction
   x y .01 sub la_fonction sub
   .02 div
end
} def

%% syntaxe : x y z (f) facteur
/facteur {
4 dict begin
   /fonction exch cvx def
   /z exch def
   /y exch def
   /x exch def
   /dfx x y (fonction) diffx def
   /dfy x y (fonction) diffy def
   /ca 
      zp z sub
      yp y sub dfy mul sub
      xp x sub dfx mul sub
   def
   /cb
      dfx dup mul
      dfy dup mul add
      1 add
      sqrt
   def
   /cc
      z zp sub dup mul
      y yp sub dup mul add
      x xp sub dup mul add
      sqrt
   def
   kp ca mul 
   cc dup dup mul mul cb mul div
end
} def

/dx 
   xmax xmin sub nx div 
def
/dy
   ymax ymin sub ny div 
def

/fillstyle {
   [ 1 .65 0 ] {kp coeff_couleur abs mul mul} apply
   aload pop setrgbcolor
%   setgray 
   fill
} def

/i 0 def
nx 1 sub {
   /xt xmin i dx mul add def
   /xt+1 xt dx add def
   /j 0 def
   ny 1 sub {
      /yt ymin j dy mul add def
      /yt+1 yt dy add def
      /coeff_couleur xt yt 2 copy f (f) facteur def
      
      xt yt 2 copy f XYZtoXY /A0 defpoint
      xt yt+1 2 copy f XYZtoXY /A1 defpoint
      xt+1 yt+1 2 copy f XYZtoXY /A2 defpoint
      xt+1 yt 2 copy f XYZtoXY /A3 defpoint
      
      gsave
   	 [A0 A1 A2 A3 A0] polygone*
      grestore
      /j j 1 add store
   } repeat
   /i i 1 add store
} repeat