// Last edited on 2013-02-02 21:23:59 by stolfilocal
// Various devices to show stress in solids. 

#include "arrow.inc"

#declare stress_particle_wall_rel_thk = 0.001; // Relative thickness of particle wall.
#declare stress_rel_eps = 0.0001;     // Relative fudge to avoid coincident surfaces.

#macro stress_particle(RP,RF txp)
  // A spherical shell with texture {txp}, and a fog ball.
  // with radius {RF} (may be larger than {RP}).

  #local eps = RP*stress_rel_eps;
  #local thk = RP*stress_particle_wall_rel_thk;

  union{
    difference{ 
      sphere{ 0,RP }
      sphere{ 0,RP-thk }
      texture{ txp }
    }
    #if (RF > 0)
      #if (RF <= RP + eps) 
        // If {RF} is less than {RP}, make it less than {RP-thk}:
        #local RF = min(RF, RP - thk - eps);
      #end
      sphere{ 0,RF 
        hollow
        texture{ pigment{ color rgb <1,1,1> filter 1 } finish { diffuse 0 ambient 0 } } 
        interior{
          media{
            // scattering{ 1, 2*< 0.100, 0.100, 0.150 > }
            absorption 2*< 0.100, 0.100, 0.150 > 
            emission 2*< 0.100, 0.100, 0.150 > 
          }
        }
      }
    #end
  }
#end

#macro stress_particle_cutter(RP,un,uct,vn,vct, txp)
  // An object used to get cutaway views of a particle with radius {PR} and
  // center at the origin. The cut faces are perpendicular to {un,vn} and 
  // removes anything that is beyond {uct} in direction {un}
  // and beyong {vct} in direction {vn}.
  intersection{
    sphere{ 0,1.1*RP }
    union{
      plane{ un,0 translate uct*un }
      plane{ vn,0 translate vct*vn }
    }
    texture{ txp }
  }
#end

#macro stress_surfel(un, RP,srad, txd)
  // A surface element on the boundary of the particle
  // with normal direction along {un}.
  // {un}    Patch normal.
  // {RP}    Radius of particle.
  // {srad}  Angular radius of patch (radians).
  // {txd}   Texture for surface patch.
  
  #local eps = RP*stress_rel_eps;
  
  // Disk parameters:
  #local dthk = RP*stress_particle_wall_rel_thk/cos(srad);
  #local dctr = RP*cos(srad)*un;  // Center of outer nominal surface of surfel.
  #local drad = RP*sin(srad);
  
  // Shrink surfel by {eps} to avoid coincident surfaces:
  #local dp0 = dctr - (dthk - eps)*un;
  #local dp1 = dctr - eps*un;
  #local dre = drad - eps;
  
  cylinder{ dp0, dp1, dre texture{ txd } }
#end

#macro stress_surfel_puncher(un, RP,srad, txp)
  // A puncher to open up space on the particle's boundary 
  // in order to accomodate the surface element along direction {un}.
  // {un}    Patch normal.
  // {RP}    Radius of particle.
  // {srad}  Angular radius of patch (radians).
  // {txp}   Texture of particle hole's walls.
  #local eps = RP*stress_rel_eps;
  #local un = vnormalize(un);
  #local dthk = RP*stress_particle_wall_rel_thk/cos(srad);
  #local dctr = RP*cos(srad)*un;  // Center of outer surface of surfel.
  #local drad = RP*sin(srad);
  
  // Enlarge hole by {eps} to avoid coincident surfaces:
  #local dp0 = dctr - (dthk + eps)*un;
  #local dp1 = (RP + eps)*un;
  #local dre = drad + eps;
  cylinder{ dp0, dp1, dre texture{ txp } }
#end

#macro stress_vectors(un,vn, vf, RP,srad,LA,vmax,rar,nrng, txf,txg)
  // One or more vectors showing the stress on a surface element
  // with normal direction along {un} on the boundary of a particle.
  // If {nrng} is zero draws a single arrow at the center fo the surfel
  // If {nrng} is positive draws {6*nrng} arrows in {nrng} rings, omitting the centeral one.
  // {un}    Patch normal.
  // {vn}    A vector not collinear with {un} (only used if {nrng > 0}).
  // {vf}    The vector to show.
  // {RP}    Radius of particle.
  // {srad}  Angular radius of patch (radians).
  // {LA}    Max arrow length for scaling.
  // {vmax}  Max field magnitude to assume for scaling.
  // {rar}   Radius of arrow shaft (controls head size too).
  // {nrng}  Number of arrow rings. If positive omits central arrow.
  // {txf}   Texture for outward stress arrows.
  // {txg}   Texture for inward stress arrows.
  
  #debug "    ENTER stress_vectors\n"

  #local eps = RP*stress_rel_eps;
  
  #local un = vnormalize(un);
  #debug "      normalized {un}\n"
  #if (nrng > 0)
    // Get vectors {vn,wn} orthogonal to {un}:
    #debug "      normalizing {vn,wn}\n"
    #local vn = vn - vdot(vn,un)*un;
    #debug concat("      orthized {vn} rel {un}, norm = ", str(vlength(vn),10,8), "\n")
    #local vn = vnormalize(vn);
    #local wn = vnormalize(vcross(un,vn));
  #end
  
  // Disk parameters:
  #debug "      computing disk dims\n"
  #local dthk = RP*stress_particle_wall_rel_thk/cos(srad);
  #local dctr = (RP*cos(srad) - dthk/2)*un;  // Center of surfel, halfway through its thickness.
  #local drad = RP*sin(srad);
  
  // Arrow parameters:
  #local va_f = (LA/vmax)*vf; // Field arrow vector.
  #local sh_f = 8*rar;         // Length of head.
  #local rh_f = 3*rar;         // Radius of head.
  #local rt_f = 0.5*rar;       // Fattening radius for head.
  #local rb_f = 2.0*rar;       // Radius of little ball at base of arrow. 
  
  #debug "      starting loop\n"

  union{
    #local dt = vdot(un,va_f);
    #if (nrng > 0) #local ir = 1; #else #local ir = 0; #end
    #while (ir <= nrng)
      #debug concat("    begin ring ", str(ir,0,0), "\n")
      #if (nrng = 0)
        #local nt = 1;     // Number of arrows in ring.
        #local rrad = 0;   // Radius of ring.
        #local ph0 = 0;    // Starting phase.
      #end
      #if (nrng = 1)
        #local nt = 6;          // Number of arrows in ring.
        #local rrad = 2*drad/3; // Radius of ring.
        #local ph0 = 0;       // Starting ph0.
      #end
      #if (nrng = 2)
        #local aa = 1/sqrt(3); // Rel radius of first band.
        #if (ir = 1)
          #local nt = 6;               // Number of arrows in ring.
          #local rrad = 2*aa/3*drad;   // Radius of ring.
          #local ph0 = 0;              // Starting phase.
        #else
          #local nt = 12;               // Number of arrows in ring.
          #local rrad = (1-aa/3)*drad;  // Radius of ring.
          #local ph0 = pi/12;           // Starting phase.
        #end
      #end
      #local nib = intersection{ sphere{ 0,rb_f } plane{ -un,0 } }
      #local it = 0;
      #while (it < nt)
        #debug concat("      begin arrow ", str(it,0,0), "\n")
        #if (ir > 0)
          #local arg = (2*pi*it + ph0)/nt; // Argument angle (radians).
          #local bp = dctr + rrad*(cos(arg)*vn + sin(arg)*wn);
        #else
          #local bp = dctr;
        #end
        #if (dt >= 0)
          object{ arrow(va_f,rar,sh_f,rh_f,rt_f) translate bp texture{ txf } } 
          object{ nib translate bp texture{ txf } }
        #else
          object{ arrow(va_f,rar,sh_f,rh_f,rt_f) translate bp - va_f texture{ txg } } 
          object{ nib translate bp texture{ txg } }
        #end
        
        #local it = it + 1;
      #end
      #local ir = ir + 1;
    #end
  }
  #debug "    EXIT stress_vectors\n"
#end