%{
    @author: Fujun Luan
    @email : fl356@cornell.edu
%}

function yarn = RollPly(params, plies)
    plyNum = double(params.ply_num);
    
    yarn = [];
    new_fiber_flag = [Inf Inf Inf];
    
    e_l = params.ellipse_long; e_s = params.ellipse_short; b_r = sqrt(e_l * e_s);
    
    for plyId = 1 : plyNum
        ply = plies{plyId};
        fibers = ply.fibers;
        
        fiberNum = length(fibers);
        
        for fiberId = 1 : fiberNum            
            fiber = fibers{fiberId};
            vertices = fiber.vertices;
            vertexNum = length(vertices);
            
            for vertexId = 1 : vertexNum
                vertex = vertices{vertexId};
                z = vertex(3);
                yarn_theta = iif(params.yarn_clock_wise, -z * 2.0 * pi / params.yarn_alpha, z * 2.0 * pi / params.yarn_alpha);
                
                local_x = vertex(1) / b_r;
                local_y = vertex(2) / b_r;
                
                short_axis = normalize(ply.base_center);    
                long_axis = [-short_axis(2) short_axis(1) 0];
                
                local_p = [local_x local_y 0];
                local_x_tmp = dot(local_p, short_axis) * e_s;  
                local_y_tmp = dot(local_p, long_axis) * e_l;
                local_p = local_x_tmp * short_axis + local_y_tmp * long_axis;
                local_x = local_p(1);    local_y = local_p(2);
                
                world_x_before_ply_rotation = local_x + ply.base_center(1);
                world_y_before_ply_rotation = local_y + ply.base_center(2);
                world_x = world_x_before_ply_rotation * cos(yarn_theta) - world_y_before_ply_rotation * sin(yarn_theta);
				world_y = world_y_before_ply_rotation * cos(yarn_theta) + world_x_before_ply_rotation * sin(yarn_theta);
                    
                new_vertex = [world_x world_y z];
                
                yarn = [yarn; new_vertex];
            end
            
            yarn = [yarn; new_fiber_flag];
        end
    end

end


function dot_v = dot(v1, v2)
    dot_v = v1 * v2';
end

function norm_v = normalize(v) 
    norm_v = v ./ norm(v);
end

function out = iif(cond, a, b)
    if cond == 1
        out = a; 
    else 
        out = b;
    end
end

 