use <scad-utils/transformations.scad>
use <scad-utils/trajectory_path.scad>
use <scad-utils/trajectory.scad>
use <scad-utils/shapes.scad>
use <scad-utils/spline.scad>
use <Round-Anything/polyround.scad>


use <list-comprehension-demos/extrusion.scad>
use <list-comprehension-demos/skin.scad>
use <list-comprehension-demos/sweep.scad>

function construct_torsion_minimizing_rotations(tangents) = [
        for (i = [0:len(tangents)-2])
                rotate_from_to(tangents[i],tangents[i+1])
];

function accumulate_rotations(rotations,acc_=[]) = let(i = len(acc_))
        i ==  len(rotations) ? acc_ :
        accumulate_rotations(rotations,
                i == 0 ? [rotations[0]] : concat(acc_, [ rotations[i] * acc_[i-1] ])
        );

function construct_transform_path3(path) = let(
        l = len(path),
        tangents = [ for (i=[0:l-1]) tangent_path(path, i)],
        local_rotations = construct_torsion_minimizing_rotations(concat([[0,0,1]],tangents)),
        rotations = accumulate_rotations(concat(local_rotations))
) [ for (i = [0:l-1]) construct_rt(rotations[i], path[i]) ]; 
    
function perpendicular(v, left) = let(s=left?-1:1)
    (v[0] == 0) ? [s*v[1], 0] : (v[1] == 0) ? [0, -s*v[0]] : [s, -s*v[0] / v[1]];


function thickify(path, thickness=0.05, left=true) = let(
    n = len(path),
    tangents = concat(
        [path[1]-path[0]],
        [for (i=[1:1:n-2]) (path[i+1]-path[i-1])/2 ],
        [path[n-1]-path[n-2]]
    ),
    orthogonal = [for (i=[0:1:n-1]) perpendicular(tangents[i],left)],
    shifted = [for (i=[0:1:n-1]) path[i]+orthogonal[i]/norm(orthogonal[i])*thickness]
) concat(path, [for (i=[n-1:-1:0]) shifted[i]]);
    
function dot_product(v1,v2) = v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];
    
function angle(v1,v2) = 
    sign(dot_product([0, 0, 1], cross(v1,v2))) *
    acos(dot_product(v1,v2)/(norm(v1)*norm(v2)));

function ellipsis_transformation(v1,v2) = 
    translation(v1) * 
    rotation(axis=[0,0,angle([0, 1, 0], v2-v1)]) *
    rotation(axis=[270,0,0]) 
    ;

function ellipsis_trajectory(a, b, n=360) = let(
    path = [for (i=[0:360/n:359])  [a * cos(i), b * sin(i), 0]],
    l = len(path)
) concat([
        for (i=[0:1:len(path)-2])
        ellipsis_transformation(path[i], path[i+1])
    ], [
        ellipsis_transformation(path[l-1], path[0])
    ]
);


module shell () {
    PI = 3.141592653589;
    PATH_STEPS = 400;
    UNIT = 10;
    epsilon = 0.001;
    
    function f(a, b, t) = [a * cos(t), b * sin(t), 0];

    NSTEPS=50;
    //path_shape = [for (i=[0:360/NSTEPS:360]) f(UNIT*1.5,UNIT, i)];
    //path = construct_transform_path3(path_shape);
    // path2 is equivalent to path, except that it seems more smooth
    path2 = ellipsis_trajectory(UNIT*1.5, UNIT, n=50);

    
    function distance(alpha) = (alpha < 0.5) ? alpha/0.5 : (1-alpha)/0.5;
    function top_distance(alpha) = (alpha < 0.5)?(alpha/0.5)^6:1;
    function bottom_distance(alpha) = (alpha > 0.5)?sqrt(((1-alpha)/0.5)):1;

    
    function spline_profile(shape) = 
        let(
            spl1=spline_args(shape),
            points=[for(t=[0:0.2:len(shape)-1]) spline(spl1, t)],
            points2d=[for (i=[0:1:len(points)-1]) [points[i][0], points[i][1]]]
        ) points2d;

    function shell_profile(alpha) = let(
        d = distance(alpha),
        td = top_distance(alpha),
        bd = bottom_distance(alpha),
        width = max(epsilon, 30*((1-td)*d^4 + bd*td*d)),
        height = max(epsilon, 50*d^5),
        xalpha=0.8,
        yalpha=0.3,
        shape=[
            [0,0,0],
            [width*xalpha, xalpha*height*yalpha, 0],
            [width, height, 0]
        ],
        path = spline_profile(shape)
    ) thickify(path, thickness=1, left=false);
            
    function ring_profile(alpha) = let(
        d = distance(alpha),
        td = top_distance(alpha),
        bd = bottom_distance(alpha),
        width = max(epsilon, 30*((1-td)*d^4+bd*td*d)),
        height = max(epsilon, 50*d^5),
        path=[
            [width, height, 0],
            [width+2.5, height, 0],
            [width+5, height, 0]
        ]
    ) thickify(path, thickness=1, left=false);
    
    function ring_position(alpha) = let(
        d = distance(alpha),
        td = top_distance(alpha),
        bd = bottom_distance(alpha),
        width = max(epsilon, 30*((1-td)*d^4+bd*td*d)),
        height = max(epsilon, 50*d^5)
    ) [width, height, 0];



    let() {
        ring_positions = [for (i=[0:len(path2)-1]) transform(path2[i], [ring_position(i/len(path2))])[0]];
        echo(ring_positions);
        ring_transforms = construct_transform_path3(ring_positions);
        ring_transformed_profiles = [for (i=[0:1:len(ring_transforms)-1]) transform(ring_transforms[i], thickify([[0,0], [5, 0]], thickness=1, left=false))];
        skin(ring_transformed_profiles);

        trans = [for (i=[0:len(path2)-1]) transform(path2[i], shell_profile(i/len(path2)))];
        skin(trans);

        *let() {
            trans2 = [for (i=[0:len(path2)-1]) transform(path2[i], ring_profile(i/len(path2)))];
            echo(trans2[0]);
            for (slice = trans2) {
                for (n = slice)
                    translate(n)sphere(0.5);
            }
            *skin(trans2);
        }
    }
}

shell();
//polygon(spline_profile(0.2, 0.2, 0.9, 0.1));
*translate([-7,-20, 36])scale([0.15, 0.15, 0.15])rotate([90,0,0])shell();
*#import("./face.stl");