use <Naca_sweep.scad>
// Demo code for affine intersection path of two linear extruded shapes 
// 

// view switches
showproduct = true; 
showintersection = true; 
showshapes=true; 


// definition of the two shapes. !! 
// In case of instabilities (parallel faces rotate shape by a small angle
shape1 = Tx(0, circle(10, 300)); 
shape2 = Tx(-3.75, Rz(-0.001, circle(13.9, 3))); // triangle
//shape2 = Tx(0, Rz(-0.1, circle(10.01, 4)));      // square
//shape2 = Tx(-1.5, Rz(-0.1, circle(11.51, 5)));   // pentagon
//shape2 = Tx(0, Rz(-0.1, circle(10.01, 6)));      // hexagon
//shape2 = Tx(0, Rz(-0.1, circle(10.01, 8)));      // octagon

w1 = 45;  // angle shape1
w2 = -w1;   // angle shape2

// display shapes
if(showshapes)
{
  #rotate([w1, 0, 0])linear_extrude(100, center= true) polygon(vec2(shape1)); 
  #rotate([w2, 0, 0])linear_extrude(100, center= true) polygon(vec2(shape2)); 
}  
   
// calc intersection paths   
res1 = intersect( shape1, shape2, w1, w2); 
//res2 = intersect( shape2, shape1, w2, w1);  // vice versa

ma= Max(res1); 
path3=select(res1, ma[0]); 

function select(L, thres) = let(L_= Rx(45,L)) 
let(L=
 [for(x=L_)  if(
 x[0]<=0 && x[1]<=0 && x[2]<0  ||
 x[0]<=0 && x[1]>=0 && x[2]>=0 || 
 x[0]>=0 && x[1]>=0 && x[2]>=0 ||
 x[0]>=0 && x[1]<=0 && x[2]<0 ) x]) 
 Rx(-45, L); 

//// define intersection path "visually" with
//// showintersection = true
//path3 = Rz(0,concat(mid(res1, 0, 30), 
//mid(res1, 61, 90),
//mid(res1, 152, 182),
//mid(res1, 215, 243)
//)); 
//color("green")   Points(path3, r=.5); 

 
 
 // show intersection path
if(showintersection)
{
  color("red")   Points(res1, r=.5);   // all points
  color("blue")   Points(path3, r=.5); // selected points
}  

h = 20;    // height
th=2;      // thickness

// compute scaling
max = Max(path3);   
sx=(max[0]-th)/max[0]; 
sy=(max[1]-th)/max[1]; 

// prepare paths for circular sweep
resA = Tz(h/2, path3);      // outer up
resB = S(sx, sy,1, resA);  // inner up
resC = Tz(-h, resB);      // inner down
resD = Tz(-h, resA);      // outer down

// prepare circular sweep
path = [for(i=[0:len(resA)-1])[resA[i],resB[i], resC[i], resD[i]]]; 
dx = 22;  // distance of the two views

if (showproduct)
{
  translate([dx, 0, 0])sweep(path, close=true); 
  rotate([0,0,180])translate([dx, 0, 0])sweep(path, close=true); 
}
$vpr = showproduct?[w1, 0, 0]:$vpr; 

//res1_ = mid(res1, 32, 63); 


/////// library code part /////////////////////////////

// intersection point of line AB and plane P1P2P3
function line_plane_intersection(A, B, P1, P2, P3) =
    let(
        N = cross(P2 - P1, P3 - P1), // normal of plane
        d = N*P1,                    // plane equation N * P = d
        AB = B - A,                  // line direction vector
        t_numer = d - N*A,
        t_denom = N*AB,
        t = t_numer/ t_denom
    )
    assert(t_denom != 0, "N,AB are parallel") 
    (A + t * AB); 
    

// distance between point P and line L1,L2
function distance_point_to_line(L1, L2, P) =
  let(
      v = L2 - L1, 
      w = P - L1, 
      proj_length = (w*v) / (v*v), // Projektion von w auf v
      proj_point = L1 + proj_length * v, // Lotfußpunkt auf der Gerade
      distance = norm(P - proj_point) // Lotrechter Abstand
  ) distance;

// construct a circle shape
function circle(r, $fn=$fn?$fn:360/$fa)=  [for(i=[0:360/$fn:359.99]) r*[cos(i), sin(i), 0]]; 

// rotate vector or list (of list of...) vectors around x, y, or z
function Rx(x, A) = A[0][0]!=undef?[for(i=A) Rx(x, i)]:
    A*[[1, 0, 0], [0, cos(x), sin(x)], [0, -sin(x), cos(x)]]; 
function Ry(y, A) = A[0][0]!=undef?[for(i=A) Ry(y, i)]:
    A*[[cos(y), 0, sin(y)], [0, 1, 0], [-sin(y), 0, cos(y)]]; 
function Rz(z, A) = A[0][0]!=undef? [for(i=A) Rz(z, i)]: 
    A*[[cos(z), sin(z), 0], [-sin(z), cos(z), 0], [0, 0, 1]]; 
    
// translate vector or list (of list of...) vectors with x
function Tx(x=0, v=undef) = v[0][0]!=undef?[for (i=v) Tx(x,i)]:v+[x,0,0]; 
function Ty(y=0, v=undef) = v[0][0]!=undef?[for (i=v) Ty(y,i)]:v+[0,y,0]; 

// reduce 3D point list to 2D
function vec2(v) = v[0][0]==undef?[v[0], v[1]]:[for(a=v) vec2(a)]; 

// show list of points as point cloud 
module Points(L, r)  assert(is_list(L) && is_list(L[0])) 
   for(p=L) translate(p) sphere(r); 

// intersect two 2D shapes at an angle w
function intersect(Sh1, Sh2, w1=w1, w2=w2) = 
let(Sh1 = Rx(w1, Sh1), Sh2 = Rx(w2, Sh2))
let(N1 =  Rx(w1, [0,0,1]), N2 = Rx(w2, [0,0,1]))
let(L1 = len(Sh1), L2=len(Sh2))
let(sign = L1<L2)
let(result =
  [
    for(i = [0:L2-1], j = [0:L1-1])
    // define the current face of Sh2
    let(P1 = Sh2[i], P2=Sh2[(i+1)%L2]) 
    let(P3 = P1+N2, P4 = P2+N2)
    let(d = norm(P1-P2))
    
    // calc the intersection of all lines of Sh1 with current face of Sh2
    let(A = Sh1[j], B = A+N1)   
      let(I = line_plane_intersection(A, B, P1, P2, P3))  
      let(dist1 = distance_point_to_line(P1, P3, I))
      let(dist2 = distance_point_to_line(P2, P4, I))
    // filter points that closely touch the face
      if(dist1+dist2<=d*1.0001) 
      I
  ])
  result; 
; 

// max() vector of a multidim list
 function Max(L) = 
   let(dim = is_list(L[0])?len(L[0]):1) 
   [
    for(d=[0:dim-1]) let(x = [for(i=L)i[d]]) max(x)];