// Last edited on 2023-06-12 20:49:07 by stolfi

// #include "curvsamp_dir6.inc"
#include "curvsamp_dir8.inc"
// #include "curvsamp_dir12.inc"
#if (bis)
  #include "curvsamp_dir24bis.inc"
  #include "curvsamp_dir48bis.inc"
#else
  #include "curvsamp_dir24.inc"
  #include "curvsamp_dir48.inc"
#end

#macro curvsamp_show_dirs(Rs,Rt,rt, bis, show_axes)
  // Displays a set of directions defined in "curvsamp_dir{48,6,12}.inc"
  
  //   {Rs} radius of shell of directions
  //   {Rt} major radius of tori.
  //   {rt} minor radius of tori.
  
  union{
    // object{ curvsamp_show_some_dirs(6, dir6,   Rs,Rt,rt, show_axes, tx_dir6) }
    object{ curvsamp_show_some_dirs(8, dir8,   Rs,Rt,rt, show_axes, tx_dir8) }
    // object{ curvsamp_show_some_dirs(12, dir12,   Rs,Rt,rt, show_axes, tx_dir12) }
    #if (bis)
      object{ curvsamp_show_some_dirs(24, dir24bis, Rs,Rt,rt, show_axes, tx_dir24) }
      object{ curvsamp_show_some_dirs(48, dir48bis, Rs,Rt,rt, show_axes, tx_dir48) }
    #else
      object{ curvsamp_show_some_dirs(24, dir24, Rs,Rt,rt, show_axes, tx_dir24) }
      object{ curvsamp_show_some_dirs(48, dir48, Rs,Rt,rt, show_axes, tx_dir48) }
    #end
  }
#end

#macro curvsamp_show_some_dirs(n,dirs, Rs,Rt,rt, show_axes, tx)
  // Display set of directions {dirs[0..n-1]}

  //   {n}             number of directions
  //   {dirs[0..n-1]}  array of unit directions
  //   {Rs} radius of  shell of directions
  //   {Rt}            major radius of torus.
  //   {rt}            minor radius of torus.
  
  #local show_skel = 0;
  
  union{
    #local k = 0;
    #while (k < n)
      #local ctr = Rs*dirs[k];
      object{ curvsamp_torus(Rt,rt, ctr, dirs[k], show_axes, tx) }
      #local k = k + 1;
    #end
    
    #if (show_skel)
      curvsamp_show_skeleton(n, dirs, Rs,rt/5, show_axes, tx)
    #end 
  }

#end

#macro curvsamp_show_skeleton(n, dirs, Rs,rb, show_axes, tx)

  union{

    #local k = 0;
    #while (k < n)
      #local ctr = Rs*dirs[k];
      sphere{ ctr, rb texture{ tx } }
      #local j = 0;
      #local d1 = 99999; #local j1 = -1;
      #local d2 = 99999; #local j2 = -1;
      #local d3 = 99999; #local j3 = -1;
      #local d4 = 99999; #local j4 = -1;
      #while (j < n)
        #if (j != k)
          #local dj = vlength(dirs[k] - dirs[j]);
          #if (dj < d1)
            #local d4 = d3; #local j4 = j3;
            #local d3 = d2; #local j3 = j2;
            #local d2 = d1; #local j2 = j1;
            #local d1 = dj; #local j1 = j;
          #elseif (dj < d2)
            #local d4 = d3; #local j4 = j3;
            #local d3 = d2; #local j3 = j2;
            #local d2 = dj; #local j2 = j;
          #elseif (dj < d3)
            #local d4 = d3; #local j4 = j3;
            #local d3 = dj; #local j3 = j;
          #elseif (dj < d4)
            #local d4 = dj; #local j4 = j;
          #end
        #end
        #local j = j + 1;
      #end
      #debug concat("  jnear = ", str(j1,0,0), " ", str(j2,0,0), " ", str(j3,0,0), "\n")
      #if (j1 < k) 
        cylinder{ Rs*dirs[k], Rs*dirs[j1], rb texture{ tx } }
      #end
      #if (j2 < k) 
        cylinder{ Rs*dirs[k], Rs*dirs[j2], rb texture{ tx } }
      #end
      #if (j3 < k) 
        cylinder{ Rs*dirs[k], Rs*dirs[j3], rb texture{ tx } }
      #end
       #if ((j4 < k) & (d4 < 1.01*d1))
        cylinder{ Rs*dirs[k], Rs*dirs[j4], rb texture{ tx } }
      #end
       
      #local k = k + 1;
    #end
  }

#end