Wasn't it Hans-Werner who wrote:
>Can I use POV-Ray also for doing area and volume calculations ?
>Here is an example:

>Suppose the mathematical function y=exp(-x^2). Greetings to Gauss !
>Now this function rotates around the z-Axis.
>Now I cut off a cylinder from this volume.
>Maybe along the x-Axis.
>Is it possible to calculate the volume that was cutted off ?

You can find a numerical approximation to the volume by firing lots of
trace() calls at it and adding up the volumes of the columns. The
smaller the Step size between the trace() positions, the more accurate
will be the calculated volume.

//-----------------
// Create the shape

#declare E = 2.71828;
#declare THING = intersection {
  isosurface {
    function { y-pow(E, -(x*x+z*z)) }
        max_gradient 1.4
        contained_by{box{<-3,0,-3><3,1,3>}}
    }
  cylinder {-x*3,x*3,1}  
}

// ----------------------------------------
// This bit is so that you can check that the
// THING is the shape that you expected it to be.
// It's not necessary for the volume calculation

camera { location  <3, 1, -4> look_at <0, 0, 0> angle 50}
light_source {<-100,200,-100> colour rgb 1}
object {THING
  pigment {rgb .9}
  finish {phong 0.5 phong_size 10}
}

// ----------------------------------------
// Set the range over which to iterate

#declare Xmin = min_extent(THING).x;
#declare Xmax = max_extent(THING).x;
#declare Zmin = min_extent(THING).z;
#declare Zmax = max_extent(THING).z;

// ----------------------------------------
// Trace the rays

#declare Norm = <0,0,0>;
#declare Total=0;
#declare Step = 0.05; // Set smaller for more accuracy
#declare X=Xmin;
#while (X<Xmax)
  #declare Z=Zmin;
  #while (Z<Zmax)
    #declare P1=trace(THING, <X,100,Z>, -y);
    #declare P2=trace(THING, <X,-100,Z>, y, Norm);
    // Test to see if the trace hit anything
    #if (vlength(Norm)>0)
      #declare Total = Total + P1.y - P2.y;
    #end
        #declare Z=Z+Step;
  #end
  #declare X=X+Step;
#end

//------------------------------------------
// Calculate the volume

#declare Volume = Total*Step*Step;
#debug concat ("Volume = ", str(Volume,-1,-1), " cubic POV units\n")

Post a reply to this message