// Last edited on 2019-04-13 01:11:54 by stolfilocal

#macro wave_all(L,R,N,sym,kind,tt)
  // A traveling wave illustration of symmetry {sym} and type {kind}
  // at time {tt}.
  
  // The medium is is illustrated as a solid cylinder if {sym=0}
  // or spherical cone if {sym=1}, that starts at the origin
  // and has the X axis as it symetry axis.  It extends {L} units
  // along the X-axis, and has maximum radius {R}.  If {sym=1} and
  // {R} is {L} or more, the medium is the whole ball of radius {L}.
  // The positive octant is cut away for better clarity.
  //
  // The medium is sliced into {N} parallel layers (plane or spherical)
  // filled with a translucent medium of varying color, to indicate the
  // field's strength.
  
  // TO DO: extend to vector-valued waves.
  
  // TO DO: option to show the unidimensional function graph {G(d,t)}.
  
  // Domain clipping parameter:
  #if (sym = 0)
    // Cylindrical domain: {Rcut} is cylinder radius.
    #local Rcut = R;
  #elseif (sym = 1)
    // Spherical domain: {Rcut} is slope of clipping cone, or 0 if none..
    #if (R < L)
      #local Rcut = R/sqrt(L*L - R*R);
    #else
      #local Rcut = 0.0; // No cone cutting.
    #end
  #else
    #debug concat("!! sym(1) = ", str(sym,8,8), "\n")
  #end
  
  union{
    #for(i,0,N-1)
      #local xx = L*(i+0.5)/N;  // Linear displacement from origin to mid-layer.
      
      #local vv = wave_scalar_value(L,sym,kind, xx,tt);

      #local xlo = L*(i+0.01)/N; // Min X/radius of layer.
      #local xhi = L*(i+0.99)/N; // Max X/radius of layer.
      object{ wave_layer(xlo,xhi,Rcut,sym, vv) }
    #end
  }
  
#end

#macro wave_layer(xlo,xhi,Rcut,sym, vv)
  #local thk = xhi-xlo; // Layer thickness.
  #local eps = 0.1*thk; // Fudge amount.
  // Define the layer object {layer} and the extent {xmin,xmax,rmax} 
  // of the box for cutting away the poitive octant:
  #if (sym = 0)
    // Cylindrical domain, plane wave:
    #local R = Rcut;
    #local cutter = box{ < xlo-eps, 0, 0 >, < xhi+eps, +1.02*R, +1.02*R > }
    #local layer = cylinder{ < xlo, 0, 0>, <xhi, 0, 0>, R }
    #local xmin = xlo - eps;
    #local xmax = xhi + eps;
    #local rmax = R + eps;
  #elseif (sym = 1)
    // Spherical domain, spherical wave:
    #local shell = difference{ sphere{ <0,0,0>, xhi } sphere{ <0,0,0>, xlo } }
    #if (Rcut = 0.0)
      #local R = xhi; // Cylinder radius.
      #local layer = obect{ shell }
      #local xmin = 0.0;
      #local xmax = xhi + eps;
      #local rmax = xhi + eps;
    #else
      #local xmin = xlo/sqrt(1.0 - Rcut*Rcut) - eps;
      #if (xmin < 0) #local xmin = 0; #end
      #local xmax = xhi + eps;
      #local rmax = Rcut*xhi + eps;
      #local layer = 
        intersection{
          object{ shell }
          cone{ < xmin,0,0>, xmin*Rcut, <xmax,0,0>, xmax*Rcut }
        }
    #end
  #else
    #debug concat("!! sym(2) = ", str(sym,8,8), "\n")
  #end
  #local cutter = box{ < xmin, 0, 0 >, < xmax, rmax, rmax > }
  difference{ 
    object{ layer }
    object{ cutter }
    wave_layer_texture(vv)
  }
#end

#macro wave_scalar_value(L,sym,kind, xx,tt)
  #local ampl = 1.0; // Max amplitude.
  // Compute the basic value:
  #local vv = 999.0; // Just in case.
  
  #if (kind = 0)
    // Sinusoidal plane wave
    #local lamb = L/5;
    #local speed = 2*L/3;
    #local phs = pi/3;
    #local tau = 0.0; // Means infinte (no decay).
    #local vv = wave_scalar_cosine_decay(lamb, phs, tau, xx - speed*tt);
  #end
  
  #if (kind = 1)
    // Generic traveling wave.
    #local gdev = L/10;      // Deviation of gaussian bell envelope.
    #local lamblo = L/5;     // Low freq lamb.
    #local lambhi = L/10;    // High freq lamb.
    #local speed = 2*L/3;    // Traveling speed.
    #local vv = wave_scalar_chirp(gdev, lamblo, lambhi, xx - speed*tt);
  #end
  
  #if (kind = 3)
    // Standing wave, generic time and space factors.
    #local gdev = L/10;      // Deviation of gaussian bell envelope.
    #local lamblo = L/5;     // Low freq lamb.
    #local lambhi = L/10;    // High freq lamb.
    #local sfactor = wave_scalar_chirp(gdev, lamblo, lambhi, xx - L/2); // Space factor.
    
    #local period = 4;       // Period of pulsation.
    #local rh3 = 0.6;        // Rel strength of harmonic 3.
    
    // Time factor, fundamental:
    #local ttau1 = 1.5;     // Characteristic decay time of fundamental.
    #local tpha1 = 0.0;     // Phase of fundamental.
    #local tamp1 = sqrt(1 - rh3*rh3); // Rel amp of fundamental.
    #local tfac1 = tamp1*wave_scalar_cosine_decay(period, tpha1, ttau1, zz); // Time factor (fund).
    
    // Time factor, harmonic 3:
    #local ttau3 = 0.7;     // Characteristic decay time of harmonic 3.
    #local tpha3 = 0.0;     // Phase of harmonic 3.
    #local tamp3 = rh3;     // Rel amp of harmonic 3.
    #local tfac3 = tamp3*wave_scalar_cosine_decay(period/3, tpha3, ttau3, zz); // Time factor (3rd).
    
    #local tfactor = tfac1+tfac3;
    
    #local vv = sfactor*tfactor;
  #end
    
  #if (kind = 4)
    // Generic planar wave.
    // Decaying wave in {+X} direction and decaying reflected at {X=L},
    // both with strong harmonic 3.

    #local lamb = L/5;  // Wavelength of fundamental.
    #local rh3 = 0.6;              // Rel strength of harmonic 3 at source.
    #local rh1 = sqrt(1-rh3*rh3);  // Rel strength of fundamental at source.
    #local rcoef31 = 0.5; // Reflection coef of harmonic 3.
   
    // Fundamental wave:
    #local rcof1 = 0.8;    // Reflection coef of fundamental.
    #local famp1 = rh1;    // Rel. amp. of fundamental at source.
    #local fpha1 = 0.0;    // Phase shift of fundamental.
    #local tau1 = 2.0*L;   // Charac decay length of fundamental.
    #local fwav1 = wave_scalar_cosine_decay(lamb, fpha1, tau, xx - speed*tt);  // Forward wave:
    #local rwav1 = rcof1*wave_scalar_cosine_decay(lamb, fpha1, tau1, 2*L - xx + speed*tt);  // Reflected wave:
    #local wav1 = famp1*(fwav1-rwav1); 
    
    // Harmonic 3 wave:
    #local rcof3 = 0.5;   // Reflection coef of harmonic 3.
    #local famp3 = rh3;   // Rel. amp. of harmonic 3 at source.
    #local fpha3 = 0.0;   // Phase shift of harmonic 3.
    #local tau3 = 3.0*L;  // Charac decay length of harmonic 3.
    #local fwav3 = wave_scalar_cosine_decay(lamb, fpha3, tau, xx - speed*tt);  // Forward wave:
    #local rwav3 = rcof3*wave_scalar_cosine_decay(lamb, fpha3, tau3, 2*L - xx + speed*tt);  // Reflected wave:
    #local wav3 = famp3*(fwav3-rwav3); 

    #local vv = wav1+wav3;
  #end
  ampl*vv
#end
    
#macro wave_scalar_cosine_decay(lamb, phs, tau, zz)
  // A cosinoid of {zz} with
  // period (wavelength) {lamb}, backwards angular 
  // shift {phs}, decay length {tau}.
  // If {tau} is zero, means infinite (no decay)
  #if (tau = 0.0)
    #local ramp = 1.0;
  #else 
    #local ramp = exp(-zz/tau);
  #end
  ramp*cos(2*pi*zz/lamb + phs)
#end

#macro wave_scalar_chirp(gdev, lamblo, lambhi, zz)
  // A chirp wave with Gaussian bell amplitude profile,
  // with maximum at {zz=0}, deviation {gdev}, 
  // differential wavelength {lamblo} at {zz=-gdev} and
  // {lambhi} at {zz=+gdev}.
  
  #local beta = log(lambhi/lamblo)/2;                  // Eponent coefficient.
  #local gamma = sqrt(lamblo*lambhi)*beta/(2*pi*gdev); // Angular speed at {zz=0}.
  #local tt = zz/gdev;              // Relative position.
  #local ang = exp(beta*tt)/gamma;   // Angular parameter.
  #local wav = cos(ang);             // Chirp with amplitude 1.
  #local env = exp(-zz*zz/2); // Envelope.
  env*wav
#end

#macro wave_layer_texture(vv)
  #local cor = wave_layer_color(vv);
  texture{
    // pigment{ color rgb cor filter 0.80 }
    // finish{ diffuse 0.03 reflection 0.25 ambient 0.02 specular 0.00 roughness 0.005 }
    pigment{ color rgb cor }
    finish{ diffuse 0.95 reflection 0.00 ambient 0.05 specular 0.00 roughness 0.005 }
  }
#end

#macro wave_layer_color(vv)
  // Color for a layer with value {vv} (from -1 to +1).
  #if (vv > 0)
    #local cor = < 1.000, 1.000 - 0.2*vv, 1.000 - 0.5*vv >;
  #elseif (vv < 0.0)
    #local cor = < 1.000 + 0.5*vv, 1.000 + 0.3*vv, 1.000 >;
  #else
    #local cor = < 1.000, 1.000, 1.000 >;
  #end
  // #local cor0 = < 0.000, 0.880, 0.500 >; // 0.5830
  // #local cor1 = < 0.200, 0.715, 1.000 >; // 0.5970
  // #local cor = rgb cor0*(1-sp) + cor1*sp;
  cor
#end

// #macro vortex_layer_texture(sp)
//   // Tetxure for tube of layer {sp} (from 0 to 1).
//   #debug "ENTER vortex_layer_texture\n"
//   #local col0 = vortex_layer_color(sp);
//   #local wht = <1,1,1>;
//   
//   #local mix_dust = 0.65;
//   #local col_dust = (1-mix_dust)*wht + mix_dust*col0;
//   
//   #local mix_glas = 0.12;
//   #local col_glas = (1-mix_glas)*wht + mix_glas*col0;
//   #local tx = 
//     texture{
//       pigment{ color col_glas filter 0.98 }
//       finish{ diffuse 0.00 reflection 0 ambient 0.00 specular 0.00 roughness 0.005 }
//     }
//     texture{
//       pigment{ color col_dust transmit 0.50 }
//       finish{ diffuse 0.80 reflection 0 ambient 0.20 specular 0.00 roughness 0.005 }
//     }
//   #local bogus = 0;
//   tx
//   #debug "EXIT vortex_layer_texture\n"
// #end
// 
// #macro vortex_particle_texture(sp)
//   // Tetxure for particles of layer {sp} (from 0 to 1).
//   #local col0 = vortex_layer_color(sp);
//   #local col1 = < 0.950, 0.800, 0.200 >; // 0.7775
//   #local mix = 0.50;
//   #local col = (1-mix)*col1 + mix*col0;
//   #local tx = 
//     texture{
//       pigment{ color col }
//       finish{ diffuse 0.60 reflection 0 ambient 0.40 specular 0.00 roughness 0.005 }
//     }
//   #local bogus = 0;
//     
//   tx
// #end
// 
// #macro vortex_particle_orbit_texture(sp)
//   #local col = vortex_layer_color(sp);
//   #local tx = 
//     texture{
//       pigment{ color 0.50*col }
//       finish{ diffuse 0.90 reflection 0 ambient 0.10 specular 0.00 roughness 0.005 }
//     }
//   #local bogus = 0;
//   tx
// #end