/*
Pov3 Parametric equation - Ennepers minimal surface. 
With smooth triangles. Colored by surface curvature, which somewhat
destroys the smooth triangle effect - but what the hey. Its easy 
to try a different parametric eqn. Just do a search and replace on
the 3 eqn elements.
Parse time exceeds trace time !
*/
#version 3.0
global_settings { assumed_gamma 1.0 }

#declare umin  = -2            //endpoints in the
#declare umax  =  2            //parametric rectangle
#declare vmin  = -2  
#declare vmax  =  2
#declare uiter = 40            //number of steps
#declare viter = 40
#declare iu    = (umax-umin)/uiter    //increment per step
#declare iv    = (vmax-vmin)/viter 

union{
  #declare uu = umin    
  #while (uu<=umax)                        //outer loop
    #declare vv = vmin
    #while (vv<=vmax)                      //inner loop
//POINTS      
      #declare x1=uu-(uu*uu*uu/3)+uu*vv*vv //Calculate the
      #declare y1=vv-(vv*vv*vv/3)+uu*uu*vv //four points in
      #declare z1=uu*uu-vv*vv              //the rectangle defined
      #declare uu=uu+iu                    //by uu,uu+iu,vv,vv+iv
      #declare x2=uu-(uu*uu*uu/3)+uu*vv*vv //Uses uu and vv because
      #declare y2=vv-(vv*vv*vv/3)+uu*uu*vv //u and v are built in
      #declare z2=uu*uu-vv*vv              //constants.
      #declare vv=vv+iv  
      #declare x3=uu-(uu*uu*uu/3)+uu*vv*vv
      #declare y3=vv-(vv*vv*vv/3)+uu*uu*vv
      #declare z3=uu*uu-vv*vv
      #declare uu=uu-iu
      #declare x4=uu-(uu*uu*uu/3)+uu*vv*vv
      #declare y4=vv-(vv*vv*vv/3)+uu*uu*vv
      #declare z4=uu*uu-vv*vv
      #declare vv=vv-iv
//NORMALS   
      #declare uu=uu+2*iu                   //Calculate surface normals
      #declare nx1=uu-(uu*uu*uu/3)+uu*vv*vv //to the nine rectangles in
      #declare ny1=vv-(vv*vv*vv/3)+uu*uu*vv //the immediate neighbourhood 
      #declare nz1=uu*uu-vv*vv
      #declare uu=uu-2*iu 
      #declare vv=vv+2*iv
      #declare nx2=uu-(uu*uu*uu/3)+uu*vv*vv
      #declare ny2=vv-(vv*vv*vv/3)+uu*uu*vv
      #declare nz2=uu*uu-vv*vv
      #declare vv=vv-2*iv 
      #declare uu=uu-iu
      #declare nx3=uu-(uu*uu*uu/3)+uu*vv*vv
      #declare ny3=vv-(vv*vv*vv/3)+uu*uu*vv
      #declare nz3=uu*uu-vv*vv
      #declare uu=uu+iu 
      #declare vv=vv-iv
      #declare nx4=uu-(uu*uu*uu/3)+uu*vv*vv
      #declare ny4=vv-(vv*vv*vv/3)+uu*uu*vv
      #declare nz4=uu*uu-vv*vv
      #declare vv=vv+iv 
      #declare n1=vcross(<nx2-x4,ny2-y4,nz2-z4>,<nx1-x2,ny1-y2,nz1-z2>)
      #declare n2=vcross(<x4-x3,y4-y3,z4-z3>,<nx2-x4,ny2-y4,nz2-z4>)   
      #declare n3=vcross(<nx3-x1,ny3-y1,nz3-z1>,<nx2-x4,ny2-y4,nz2-z4>)   
      #declare n4=vcross(<x3-x2,y3-y2,z3-z2>,<nx1-x2,ny1-y2,nz1-z2>) 
      #declare n5=vcross(<x1-x2,y1-y2,z1-z2>,<x3-x2,y3-y2,z3-z2>) 
      #declare n6=vcross(<nx3-x1,ny3-y1,nz3-z1>,<x4-x1,y4-y1,z4-z1>) 
      #declare n7=vcross(<nx1-x2,ny1-y2,nz1-z2>,<nx4-x1,ny4-y1,nz4-z1>) 
      #declare n8=vcross(<nx4-x1,ny4-y1,nz4-z1>,<x1-x2,y1-y2,z1-z2>)   
      #declare n9=vcross(<nx4-x1,ny4-y1,nz4-z1>,<nx3-x1,ny3-y1,nz3-z1>)  
//CURVATURE  as difference between normals - a kind of directional curvature    
      #declare ucur = vlength(vnormalize(n6)-vnormalize(n5))/2
      #declare vcur = vlength(vnormalize(n2)-vnormalize(n5))/2 
//TRIANGLES      
      object{
        smooth_triangle{
          <x1,y1,z1>,n5+n6+n8+n9,
          <x2,y2,z2>,n4+n5+n7+n8,
          <x3,y3,z3>,n1+n2+n4+n5
         }
          texture{ 
          pigment {color rgb 
//Fiddle colors for better result - multiplying curvatures 
//gives pseudo Gaussian curvature             
             <0.5+ucur*vcur*200,0.5-ucur*vcur*200,ucur/vcur-0.4>
                  }
          finish { ambient 0.2 specular 0.7 roughness 0.05 }
            }
          }
      object{
        smooth_triangle{
          <x1,y1,z1>, n5+n6+n8+n9,  
          <x3,y3,z3>, n1+n2+n4+n5   
          <x4,y4,z4>, n2+n3+n5+n6  
         }
          texture{
        pigment{ color rgb 
           <0.5+ucur*vcur*200,0.5-ucur*vcur*200 , ucur/vcur-0.4>
               }
        finish { ambient 0.2 specular 0.7 roughness 0.05 }
         }
        }
    
      #declare vv = vv+iv
   #end                 //inner loop
   #declare uu = uu+iu
 #end                   //outer loop
rotate <90,0,0>     //This is just changing the
rotate <0,90,0>     //coordinate system.
}

//Lights, camera, action
camera
{ location  < 15, 7,-15>
  direction  1*z
  look_at   < 0, 0, 0>
} 

light_source { <0, 0,-19> color rgb 1 }
light_source { <19, 2, 0> color rgb 1 }
object {box{-20,20 inverse}pigment{color rgb 0.8}}