diff --git a/Manifest.toml b/Manifest.toml
index f56c5db..9ac08f3 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.10.2"
 manifest_format = "2.0"
-project_hash = "d033fa29a644db60fc5e6aeb38e2ece2b7d0fc92"
+project_hash = "3ec2d41cb347d8597a69892af02440dbe8605a8b"
 
 [[deps.AbstractFFTs]]
 deps = ["LinearAlgebra"]
@@ -73,9 +73,9 @@ version = "7.9.0"
 
 [[deps.ArrayLayouts]]
 deps = ["FillArrays", "LinearAlgebra"]
-git-tree-sha1 = "6404a564c24a994814106c374bec893195e19bac"
+git-tree-sha1 = "0330bc3e828a05d1073553fb56f9695d73077370"
 uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
-version = "1.8.0"
+version = "1.9.1"
 weakdeps = ["SparseArrays"]
 
     [deps.ArrayLayouts.extensions]
@@ -177,9 +177,9 @@ version = "3.24.0"
 
 [[deps.ColorTypes]]
 deps = ["FixedPointNumbers", "Random"]
-git-tree-sha1 = "eb7f0f8307f71fac7c606984ea5fb2817275d6e4"
+git-tree-sha1 = "b10d0b65641d57b8b4d5e234446582de5047050d"
 uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
-version = "0.11.4"
+version = "0.11.5"
 
 [[deps.ColorVectorSpace]]
 deps = ["ColorTypes", "FixedPointNumbers", "LinearAlgebra", "Requires", "Statistics", "TensorCore"]
@@ -467,9 +467,9 @@ version = "2.13.1+0"
 
 [[deps.FreeTypeAbstraction]]
 deps = ["ColorVectorSpace", "Colors", "FreeType", "GeometryBasics"]
-git-tree-sha1 = "055626e1a35f6771fe99060e835b72ca61a52621"
+git-tree-sha1 = "2493cdfd0740015955a8e46de4ef28f49460d8bc"
 uuid = "663a7486-cb36-511b-a19d-713bb74d65c9"
-version = "0.10.1"
+version = "0.10.3"
 
 [[deps.FriBidi_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -495,9 +495,9 @@ version = "3.3.9+0"
 
 [[deps.GLMakie]]
 deps = ["ColorTypes", "Colors", "FileIO", "FixedPointNumbers", "FreeTypeAbstraction", "GLFW", "GeometryBasics", "LinearAlgebra", "Makie", "Markdown", "MeshIO", "ModernGL", "Observables", "PrecompileTools", "Printf", "ShaderAbstractions", "StaticArrays"]
-git-tree-sha1 = "7a411adf08375e01d864386fb1eaf384de5ac9e9"
+git-tree-sha1 = "b99a999cdcc7467769c7ea8cdc0978a5f059aa7b"
 uuid = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
-version = "0.9.9"
+version = "0.9.10"
 
 [[deps.GeoInterface]]
 deps = ["Extents"]
@@ -859,10 +859,10 @@ uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 version = "0.5.13"
 
 [[deps.Makie]]
-deps = ["Animations", "Base64", "CRC32c", "ColorBrewer", "ColorSchemes", "ColorTypes", "Colors", "Contour", "DelaunayTriangulation", "Distributions", "DocStringExtensions", "Downloads", "FFMPEG_jll", "FileIO", "FilePaths", "FixedPointNumbers", "Format", "FreeType", "FreeTypeAbstraction", "GeometryBasics", "GridLayoutBase", "ImageIO", "InteractiveUtils", "IntervalSets", "Isoband", "KernelDensity", "LaTeXStrings", "LinearAlgebra", "MacroTools", "MakieCore", "Markdown", "MathTeXEngine", "Observables", "OffsetArrays", "Packing", "PlotUtils", "PolygonOps", "PrecompileTools", "Printf", "REPL", "Random", "RelocatableFolders", "Scratch", "ShaderAbstractions", "Showoff", "SignedDistanceFields", "SparseArrays", "StableHashTraits", "Statistics", "StatsBase", "StatsFuns", "StructArrays", "TriplotBase", "UnicodeFun"]
-git-tree-sha1 = "27af6be179c711fb916a597b6644fbb5b80becc0"
+deps = ["Animations", "Base64", "CRC32c", "ColorBrewer", "ColorSchemes", "ColorTypes", "Colors", "Contour", "DelaunayTriangulation", "Distributions", "DocStringExtensions", "Downloads", "FFMPEG_jll", "FileIO", "FilePaths", "FixedPointNumbers", "Format", "FreeType", "FreeTypeAbstraction", "GeometryBasics", "GridLayoutBase", "ImageIO", "InteractiveUtils", "IntervalSets", "Isoband", "KernelDensity", "LaTeXStrings", "LinearAlgebra", "MacroTools", "MakieCore", "Markdown", "MathTeXEngine", "Observables", "OffsetArrays", "Packing", "PlotUtils", "PolygonOps", "PrecompileTools", "Printf", "REPL", "Random", "RelocatableFolders", "Scratch", "ShaderAbstractions", "Showoff", "SignedDistanceFields", "SparseArrays", "Statistics", "StatsBase", "StatsFuns", "StructArrays", "TriplotBase", "UnicodeFun"]
+git-tree-sha1 = "46ca613be7a1358fb93529726ea2fc28050d3ae0"
 uuid = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
-version = "0.20.8"
+version = "0.20.9"
 
 [[deps.MakieCore]]
 deps = ["Observables", "REPL"]
@@ -898,9 +898,9 @@ version = "0.4.11"
 
 [[deps.Missings]]
 deps = ["DataAPI"]
-git-tree-sha1 = "f66bdc5de519e8f8ae43bdc598782d35a25b1272"
+git-tree-sha1 = "ec4f7fbeab05d7747bdf98eb74d130a2a2ed298d"
 uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
-version = "1.1.0"
+version = "1.2.0"
 
 [[deps.Mmap]]
 uuid = "a63ad114-7e13-5084-954f-fe012c677804"
@@ -1077,12 +1077,6 @@ git-tree-sha1 = "eb3f9df2457819bf0a9019bd93cc451697a0751e"
 uuid = "2ae35dd2-176d-5d53-8349-f30d82d94d4f"
 version = "0.4.20"
 
-[[deps.PikaParser]]
-deps = ["DocStringExtensions"]
-git-tree-sha1 = "d6ff87de27ff3082131f31a714d25ab6d0a88abf"
-uuid = "3bbf5609-3e7b-44cd-8549-7c69f321e792"
-version = "0.6.1"
-
 [[deps.Pixman_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LLVMOpenMP_jll", "Libdl"]
 git-tree-sha1 = "64779bc4c9784fee475689a1752ef4d5747c5e87"
@@ -1367,12 +1361,6 @@ weakdeps = ["ChainRulesCore"]
     [deps.SpecialFunctions.extensions]
     SpecialFunctionsChainRulesCoreExt = "ChainRulesCore"
 
-[[deps.StableHashTraits]]
-deps = ["Compat", "PikaParser", "SHA", "Tables", "TupleTools"]
-git-tree-sha1 = "10dc702932fe05a0e09b8e5955f00794ea1e8b12"
-uuid = "c5dd0088-6c3f-4803-b00e-f31a60c170fa"
-version = "1.1.8"
-
 [[deps.StackViews]]
 deps = ["OffsetArrays"]
 git-tree-sha1 = "46e589465204cd0c08b4bd97385e4fa79a0c770c"
@@ -1414,9 +1402,9 @@ version = "1.7.0"
 
 [[deps.StatsBase]]
 deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
-git-tree-sha1 = "1d77abd07f617c4868c33d4f5b9e1dbb2643c9cf"
+git-tree-sha1 = "5cf7606d6cef84b543b483848d4ae08ad9832b21"
 uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
-version = "0.34.2"
+version = "0.34.3"
 
 [[deps.StatsFuns]]
 deps = ["HypergeometricFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
@@ -1511,11 +1499,6 @@ git-tree-sha1 = "4d4ed7f294cda19382ff7de4c137d24d16adc89b"
 uuid = "981d1d27-644d-49a2-9326-4793e63143c3"
 version = "0.1.0"
 
-[[deps.TupleTools]]
-git-tree-sha1 = "41d61b1c545b06279871ef1a4b5fcb2cac2191cd"
-uuid = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6"
-version = "1.5.0"
-
 [[deps.UUIDs]]
 deps = ["Random", "SHA"]
 uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
diff --git a/examples/demo_kabsch_rot.jl b/examples/demo_kabsch_rot.jl
new file mode 100644
index 0000000..1599181
--- /dev/null
+++ b/examples/demo_kabsch_rot.jl
@@ -0,0 +1,41 @@
+using Comodo
+using GLMakie
+using GeometryBasics
+using FileIO
+using Rotations
+
+#=
+In this demo a mesh is loaded, in this case a triangulated surface from an STL 
+file. Next the coordinates are rotated, and the unrotated and rotated meshes 
+are visualized. 
+=#
+
+# Loading a mesh
+fileName_mesh = joinpath(comododir(),"assets","stl","stanford_bunny_low.stl")
+M = load(fileName_mesh)
+
+# Obtain mesh faces and vertices
+F = faces(M)
+V1 = togeometrybasics_points(coordinates(M))
+F,V1 = mergevertices(F,V1)
+
+# Define a rotation tensor using Euler angles
+Q = RotXYZ(0.0,0.25*π,0.25*π)
+V2 = [GeometryBasics.Point{3, Float64}(Q*v) for v ∈ V1] 
+
+R = kabsch_rot(V2,V1)
+V3 = [GeometryBasics.Point{3, Float64}(R*v) for v ∈ V2] 
+
+# Rotate the coordinates
+fig = Figure(size = (800,800))
+ax = Axis3(fig[1, 1], aspect = :data)
+
+hp1 = poly!(ax, GeometryBasics.Mesh(V1,F), color=:green,transparency=false,shading = FastShading)
+hp2 = poly!(ax, GeometryBasics.Mesh(V2,F), strokewidth=2,color=:red,shading = FastShading)
+hp3 = wireframe!(ax, GeometryBasics.Mesh(V3,F), color=:red,linewidth=2)
+
+Legend(fig[1, 2],[hp1,hp2,hp3],["Initial","Rotated","Back rotated"])
+
+fig
+
+
diff --git a/examples/demo_sweeploft.jl b/examples/demo_sweeploft.jl
new file mode 100644
index 0000000..d10b630
--- /dev/null
+++ b/examples/demo_sweeploft.jl
@@ -0,0 +1,150 @@
+using Comodo
+using GeometryBasics
+using GLMakie
+using Rotations
+using Statistics
+using LinearAlgebra
+
+testCase = 1
+
+if testCase == 1    
+    # Define guide curve
+    nc = 51 # Number of points on guide curve
+    P = Vector{GeometryBasics.Point{3, Float64}}(undef,4)
+    P[1 ] = GeometryBasics.Point{3, Float64}( 0.0, 0.0, 0.0)
+    P[2 ] = GeometryBasics.Point{3, Float64}( 1.0, 0.0, 0.0)
+    P[3 ] = GeometryBasics.Point{3, Float64}( 1.0, 1.0, 0.0)
+    P[4 ] = GeometryBasics.Point{3, Float64}( 1.0, 1.0, 1.0)
+    Vc = nbezier(P,nc) # Get Bezier fit points
+    Vc = [vc.*10 for vc in Vc]
+    Vc,Sc = evenly_sample(Vc, nc)
+
+    # Define section curves
+    np = 25 # Number of section points
+    f(t) = 2.0 + 0.5.*sin(3*t)
+    V1 = circlepoints(f,np; dir=:acw)
+    V1,_ = evenly_sample(V1, np)
+    Q = RotXYZ(0.0,0.5*π,0.0) # Define a rotation tensor using Euler angles
+    V1 = [Q*v for v ∈ V1] # Rotate the coordinates
+
+    f(t) = 2.0 + 0.5*sin(3*t)
+    V2 = circlepoints(f,np; dir=:acw)
+    V2,_ = evenly_sample(V2, np)
+    V2 = [v2 .+ Vc[end] for v2 ∈ V2] 
+elseif testCase == 2
+        # Define guide curve
+        nc = 75 # Number of points on guide curve
+        P = Vector{GeometryBasics.Point{3, Float64}}(undef,4)
+        P[1 ] = GeometryBasics.Point{3, Float64}( 0.0, 0.0, 0.0)
+        P[2 ] = GeometryBasics.Point{3, Float64}( 1.0, 0.0, 0.0)
+        P[3 ] = GeometryBasics.Point{3, Float64}( 2.0, 0.0, 0.0)
+        P[4 ] = GeometryBasics.Point{3, Float64}( 3.0, 0.0, 0.0)
+        Vc = nbezier(P,nc) # Get Bezier fit points
+        Vc = [vc.*10 for vc in Vc]
+        Vc,Sc = evenly_sample(Vc, nc)
+    
+        # Define section curves
+        np = 35 # Number of section points
+        f(t) = 2.0 + 0.5.*sin(3*t)
+        V1 = circlepoints(f,np; dir=:acw)
+        V1,_ = evenly_sample(V1, np)
+        Q = RotXYZ(0.0,0.5*π,0.0) # Define a rotation tensor using Euler angles
+        V1 = [(Q*v) .+ Vc[1] for v ∈ V1] # Rotate the coordinates
+    
+        f(t) = 2.0 + 0.5*sin(3*t)
+        V2 = circlepoints(f,np; dir=:acw)
+        V2,_ = evenly_sample(V2, np)
+        Q = RotXYZ(0.0,0.5*π,0.0) # Define a rotation tensor using Euler angles
+        V2 = [(Q*v) .+ Vc[end] for v ∈ V2] 
+elseif testCase == 3
+    nc = 201 # Number of points on guide curve
+    r = 10.0
+    a = 4*π
+    Vc = [GeometryBasics.Point{3, Float64}(r*cos(t),r*sin(t),10.0*(t/(a/2))) for t ∈ range(0,a,nc)]
+
+    # Define section curves
+    np = 50 # Number of section points
+
+    # Section 1
+    f(t) = 1.5 + 0.5.*sin(3*t)
+    V1 = circlepoints(f,np; dir=:cw)
+    V1,_ = evenly_sample(V1, np)
+    Q = RotXYZ(0.5*π,0.0,0.0) # Define a rotation tensor using Euler angles
+    V1 = [Q*v for v ∈ V1] # Rotate the coordinates
+
+    # Ensure section is orthogonal to guide curve
+    n3_1 = Q*Vec3{Float64}(0.0,0.0,-1.0)                       
+    n2_1 = Q*Vec3{Float64}(1.0,0.0,0.0)   
+    n1_1 = normalizevector(cross(n3_1,n2_1))
+    S11 = mapreduce(permutedims,vcat,[n1_1,n2_1,n3_1])
+
+    n3_1 = normalizevector(normalizevector(Vc[2]-Vc[1]))                            
+    n2_1 = normalizevector(cross(normalizevector(V1[1]),n3_1))
+    n1_1 = normalizevector(cross(n3_1,n2_1))
+    S12 = mapreduce(permutedims,vcat,[n1_1,n2_1,n3_1])
+
+    R = RotMatrix3{Float64}(S12\S11)
+    V1 = [R*v for v ∈ V1]
+    V1= [v .+ Vc[1] for v ∈ V1] 
+
+    # Section 2
+    f(t) = 4 + 1.5*sin(5*t)
+    V2 = circlepoints(f,np; dir=:cw)
+    V2,_ = evenly_sample(V2, np)
+    Q1 = RotXYZ(0.5*π,0.0,0.0) # Define a rotation tensor using Euler angles
+    Q2 = RotXYZ(0.0,-0.25*π,0.0) # Define a rotation tensor using Euler angles
+    Q = Q2*Q1
+    V2 = [Q*v for v ∈ V2] # Rotate the coordinates
+
+    # Ensure section is orthogonal to guide curve
+    n3_2 = Q*Vec3{Float64}(0.0,0.0,-1.0)                       
+    n2_2 = Q*Vec3{Float64}(1.0,0.0,0.0)   
+    n1_2 = normalizevector(cross(n3_2,n2_2))
+    S21 = RotMatrix3{Float64}(mapreduce(permutedims,vcat,[n1_2,n2_2,n3_2]))
+
+    n3_2 = normalizevector(normalizevector(Vc[2]-Vc[1]))                            
+    n2_2 = normalizevector(cross(normalizevector(V1[1]),n3_1))
+    n1_2 = normalizevector(cross(n3_1,n2_1))
+    S22 = mapreduce(permutedims,vcat,[n1_1,n2_1,n3_1])
+
+    R = RotMatrix3{Float64}(S22\S21)    
+    V2 = [R*v for v ∈ V2] # Rotate the coordinates
+    V2 = [v .+ Vc[end] for v ∈ V2] 
+end
+
+
+
+#########
+
+# face_type=:quad
+F,V = sweeploft(Vc,V1,V2; face_type=:quad, num_twist=0)
+F = invert_faces(F)
+
+
+# Visualization
+fig = Figure(size = (800,800))
+ax = Axis3(fig[1, 1],aspect = :data,title="Swept lofting")
+
+stepRange1 = -4:1:4
+hSlider1 = Slider(fig[2, 1], range = stepRange1, startvalue = 0,linewidth=30)
+
+# scatter!(ax, Vc,markersize=8,color=:black)
+hp1 = lines!(ax, Vc,linewidth=4,color=:black)
+
+scatter!(ax, V1,markersize=8,color=:blue)
+hp2 = lines!(ax, V1,linewidth=4,color=:blue)
+
+scatter!(ax, V2,markersize=8,color=:red)
+hp3 = lines!(ax, V2,linewidth=4,color=:red)
+
+hp1 = poly!(ax, GeometryBasics.Mesh(V,F), strokecolor=:black, strokewidth=1,color=:white,transparency=false,shading = FastShading)
+# hp1 = mesh!(ax, GeometryBasics.Mesh(V,F), color=:white,transparency=false,shading = FastShading)
+
+on(hSlider1.value) do stepIndex1
+    F,V = sweeploft(Vc,V1,V2; face_type=:quad, num_twist=stepIndex1)
+    F = invert_faces(F)   
+    hp1[1] = GeometryBasics.Mesh(V,F)
+end
+
+fig 
+
diff --git a/src/Comodo.jl b/src/Comodo.jl
index cf3db50..f70840c 100644
--- a/src/Comodo.jl
+++ b/src/Comodo.jl
@@ -31,7 +31,7 @@ export circlepoints, loftlinear, trisurfslice
 export wrapindex, edgeangles, count_edge_face, boundaryedges, edges2curve
 export pointspacingmean, extrudecurve, meshgroup, ray_triangle_intersect
 export distmarch, mesh_curvature_polynomial, curve_length, evenly_sample
-export invert_faces
+export invert_faces, kabsch_rot, sweeploft
 
 # Export functions: Visualization related
 export slidercontrol
diff --git a/src/functions.jl b/src/functions.jl
index fd7141c..b36d153 100644
--- a/src/functions.jl
+++ b/src/functions.jl
@@ -1151,7 +1151,9 @@ function togeometrybasics_mesh(VM,FM)
 end
 
 
-function edgecrossproduct(F,V)    
+
+# function edgecrossproduct(F,V::Union{Array{Vec{m, T}, 1},Array{Point{m, T}, 1}})  where T<:Real where m 
+function edgecrossproduct(F,V) 
     C = Vector{GeometryBasics.Vec{3, Float64}}(undef,length(F)) # Allocate array cross-product vectors
     n =  length(F[1]) # Number of nodes per face    
     for q ∈ eachindex(F) # Loop over all faces
@@ -2772,4 +2774,129 @@ flipped, by reversing the node order for each face.
 """
 function invert_faces(F::Array{NgonFace{N, Int64}, 1}) where N    
     return map(f-> reverse(f),F) # [NgonFace{N, Int64}(reverse(f)) for f ∈ F]    
+end
+
+
+"""
+    R = kabsch_rot(V1::Array{Point{N, Float64}, 1},V2::Array{Point{N, Float64}, 1}) where N    
+
+# Description  
+Computes the rotation tensor `R` to rotate the points in `V1` to best match the 
+points in `V2`. 
+
+# Reference
+[Wolfgang Kabsch, A solution for the best rotation to relate two sets of vectors, Acta Crystallographica Section A, vol. 32, no. 5, pp. 922-923, 1976, doi: 10.1107/S0567739476001873](https://doi.org/10.1107/S0567739476001873)
+"""
+function kabsch_rot(V1::Array{Point{N, Float64}, 1},V2::Array{Point{N, Float64}, 1}) where N    
+    # Centre on means 
+    V1 = V1.-mean(V1)
+    V2 = V2.-mean(V2)
+        
+    # Compute A matrix
+    A = zeros(Float64,3,3)
+    for i = 1:3    
+        for j= 1:3
+            for q ∈ eachindex(V1)
+                @inbounds A[i,j] += V1[q][i]*V2[q][j] 
+            end
+        end
+    end
+
+    # Kabsch algorithm
+    U, _, V = svd(A)    
+    d = det(U)*det(V) #sign(det(U'*V))
+    D = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
+    D[3,3] = d 
+    return Rotations.RotMatrix3{Float64}(V*D*U')
+end
+
+
+"""
+    F,V = sweeploft(Vc,V1,V2; face_type=:quad, num_twist = 0, close_loop=true)   
+
+# Description
+This function implements swept lofting. The start curve `V1` is pulled allong the 
+guide curve `Vc` while also gradually (linearly) morphing into the end curve 
+`V2`. 
+The optional parameter `face_type` (default :quad) defines the type of mesh 
+faces uses. The same face types as `loftlinear` and `extrudecurve` are supported, 
+i.e. `:quad`, `:tri_slash`, or `tri`. 
+The optional parameter `num_twist` (default is 0) can be used to add an integer 
+number (negative or positive) of full twists to the loft. 
+Finally the optional parameter `close_loop` (default is `true`) determines if the
+section curves are deemed closed or open ended. 
+"""
+function sweeploft(Vc,V1,V2; face_type=:quad, num_twist = 0, close_loop=true)    
+    np = length(V1) # Number of section points
+    nc = length(Vc) # Number of curve steps
+
+    # Centre start/end sections arond means (should be curve start/end points)
+    V1b = [v.-Vc[1] for v ∈ V1] 
+    V2b = [v.-Vc[end] for v ∈ V2] 
+
+    # Determine rotation between sections 
+    n3_1 = facenormal([collect(1:length(V1))],V1)[1] # normalizevector(Vc[2]-Vc[1])
+    n1_1 = normalizevector(V1b[1]) 
+    n2_1 = normalizevector(cross(n1_1,n3_1))
+    n1_1 = normalizevector(cross(n3_1,n2_1))
+    S1p = mapreduce(permutedims,vcat,[n1_1,n2_1,n3_1])
+
+    n3_2 = facenormal([collect(1:length(V2))],V2)[1] # normalizevector(Vc[end]-Vc[end-1])
+    n1_2 = normalizevector(V2b[1]) 
+    n2_2 = normalizevector(cross(n1_2,n3_2))
+    n1_2 = normalizevector(cross(n3_2,n2_2))        
+    S2p = mapreduce(permutedims,vcat,[n1_2,n2_2,n3_2])
+
+    Q12 = RotMatrix3{Float64}(S1p\S2p) # nearest_rotation(S1p\S2p)
+    
+    # Rotate V2b to start orientation
+    V2b = [Q12*v for v ∈ V2b] 
+
+    # Linearly loft "alligned" sections to create draft intermediate sections
+    F,V = loftlinear(V1b,V2b;num_steps=nc,close_loop=close_loop,face_type=face_type)
+    
+    # Rotating and positioning all sections 
+    n1p = n1_1
+    for q = 1:nc # For all curve points 
+        
+        ind = (1:np) .+ (q-1)*np # Indices of the points to map
+        if q==1 # Just take the start section when we are at the start
+            V[ind] = V1
+        elseif q == nc # Just take the end section when we are at the start
+            V[ind] = V2            
+        else # Rotate and position intermediate section 
+            n3 = normalizevector(normalizevector(Vc[q]-Vc[q-1]) .+ normalizevector(Vc[q+1]-Vc[q]))                        
+            n2 = normalizevector(cross(n1p,n3))
+            n1 = normalizevector(cross(n3,n2))
+            S2 = mapreduce(permutedims,vcat,[n1,n2,n3])
+
+            Q12 = RotMatrix3{Float64}(S2\S1p)
+            V[ind] = [(Q12*v).+Vc[q] for v ∈ V[ind]]
+            n1p = n1
+        end   
+
+        if q == nc-1 # Once here, second to last, a potential rotational mismatch needs to be resolve
+            Q_fix = Rotations.RotMatrix3{Float64}(S2\S2p) # Rotation between last and second to last
+            t_a = Rotations.params(AngleAxis(Q_fix)) # Angle/axis representation
+            if dot(t_a[2:end],n3)<0 
+                β_fix = t_a[1]
+            else # Flip angle sign in this case 
+                β_fix = -t_a[1]
+            end
+
+            if β_fix>pi # Use shorter negative direction instead
+                β_fix = -(2.0*pi-β_fix)
+            end
+        
+            # Distribute the fix as around curve rotations for each intermediate step            
+            β_range = range(0,β_fix+(2.0*π*num_twist),nc) 
+            for q = 2:nc-1
+                ind = (1:np) .+ (q-1)*np
+                ns = normalizevector(normalizevector(Vc[q]-Vc[q-1]) .+ normalizevector(Vc[q+1]-Vc[q]))
+                Q = AngleAxis(β_range[q], ns[1], ns[2], ns[3])                  
+                V[ind] = [(Q*(v-Vc[q])).+Vc[q] for v ∈ V[ind]]
+            end
+        end
+    end      
+    return F,V
 end
\ No newline at end of file
diff --git a/test/Manifest.toml b/test/Manifest.toml
index a2406f6..590b772 100644
--- a/test/Manifest.toml
+++ b/test/Manifest.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.10.2"
 manifest_format = "2.0"
-project_hash = "df760090b2514e3cc5a3f02dbfe17b8fc168dacd"
+project_hash = "e9bc825cf8ca45dba6b20812af30e2dfca889490"
 
 [[deps.AbstractFFTs]]
 deps = ["LinearAlgebra"]
@@ -1127,6 +1127,12 @@ git-tree-sha1 = "9b23c31e76e333e6fb4c1595ae6afa74966a729e"
 uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 version = "2.9.4"
 
+[[deps.Quaternions]]
+deps = ["LinearAlgebra", "Random", "RealDot"]
+git-tree-sha1 = "994cc27cdacca10e68feb291673ec3a76aa2fae9"
+uuid = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0"
+version = "0.7.6"
+
 [[deps.REPL]]
 deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
 uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
@@ -1150,6 +1156,12 @@ weakdeps = ["FixedPointNumbers"]
     [deps.Ratios.extensions]
     RatiosFixedPointNumbersExt = "FixedPointNumbers"
 
+[[deps.RealDot]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "9f0a1b71baaf7650f4fa8a1d168c7fb6ee41f0c9"
+uuid = "c1ae055f-0cd5-4b69-90a6-9a35b1a98df9"
+version = "0.1.0"
+
 [[deps.RecipesBase]]
 deps = ["PrecompileTools"]
 git-tree-sha1 = "5c3d09cc4f31f5fc6af001c250bf1278733100ff"
@@ -1191,6 +1203,16 @@ git-tree-sha1 = "6ed52fdd3382cf21947b15e8870ac0ddbff736da"
 uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
 version = "0.4.0+0"
 
+[[deps.Rotations]]
+deps = ["LinearAlgebra", "Quaternions", "Random", "StaticArrays"]
+git-tree-sha1 = "2a0a5d8569f481ff8840e3b7c84bbf188db6a3fe"
+uuid = "6038ab10-8711-5258-84ad-4b1120ba62dc"
+version = "1.7.0"
+weakdeps = ["RecipesBase"]
+
+    [deps.Rotations.extensions]
+    RotationsRecipesBaseExt = "RecipesBase"
+
 [[deps.RoundingEmulator]]
 git-tree-sha1 = "40b9edad2e5287e05bd413a38f61a8ff55b9557b"
 uuid = "5eaf0fd0-dfba-4ccb-bf02-d820a40db705"
diff --git a/test/Project.toml b/test/Project.toml
index c6d1691..e47c957 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -2,5 +2,6 @@
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/runtests.jl b/test/runtests.jl
index e289b18..be1fc17 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,4 +1,4 @@
-using Test, FileIO, Comodo, Comodo.GeometryBasics, Statistics, LinearAlgebra, GLMakie
+using Test, FileIO, Comodo, Comodo.GeometryBasics, Statistics, LinearAlgebra, GLMakie, Rotations
 
 # ConnectivitySet
 
@@ -3527,7 +3527,6 @@ end
 
     # Even sampling should be nearly perfect for downsampling
     @testset "Evenly downsampling circular curve" begin
-
         # Example circle curve raw
         r = 2.25
         nc = 100    
@@ -3557,3 +3556,87 @@ end
     F_inv = [TriangleFace{Int64}(3, 2, 1),TriangleFace{Int64}(6, 5, 4)]   
     @test invert_faces(F)==F_inv    
 end
+
+@testset "kabsch_rot" begin
+    tol_level = 1e-6
+
+    M = cube(sqrt(3))
+    V1 = coordinates(M)
+    R_true = RotXYZ(0.25*π,0.25*π,0.25*π)
+    V2 = [R_true*v for v ∈ V1]
+    R_kabsch_forward = kabsch_rot(V1,V2)
+    R_kabsch_backward = kabsch_rot(V2,V1)
+    V1r = [R_kabsch_backward*v for v ∈ V2]
+    @test isapprox(R_kabsch_forward,R_true,atol=tol_level) # Able to retrieve forward rotation
+    @test isapprox(V1r,V1,atol=tol_level) # Check if backward rotation is succesful   
+end
+
+@testset "sweeploft" verbose = true begin
+    tol_level = 1e-6
+
+    # Define guide curve
+    nc = 25 # Number of points on guide curve
+    P = Vector{GeometryBasics.Point{3, Float64}}(undef,4)
+    P[1 ] = GeometryBasics.Point{3, Float64}( 0.0, 0.0, 0.0)
+    P[2 ] = GeometryBasics.Point{3, Float64}( 1.0, 0.0, 0.0)
+    P[3 ] = GeometryBasics.Point{3, Float64}( 1.0, 1.0, 0.0)
+    P[4 ] = GeometryBasics.Point{3, Float64}( 1.0, 1.0, 1.0)
+    Vc = nbezier(P,nc) # Get Bezier fit points
+    Vc = [vc.*10 for vc in Vc]
+    Vc,Sc = evenly_sample(Vc, nc)
+
+    # Define section curves
+    np = 20 # Number of section points
+    ff(t) = 2.0 + 0.5.*sin(3*t)
+    V1 = circlepoints(ff,np; dir=:acw)
+    V1,_ = evenly_sample(V1, np)
+    Q = RotXYZ(0.0,0.5*π,0.0) # Define a rotation tensor using Euler angles
+    V1 = [Q*v for v ∈ V1] # Rotate the coordinates
+
+    ff(t) = 2.0 + 0.5*sin(3*t)
+    V2 = circlepoints(ff,np; dir=:acw)
+    V2,_ = evenly_sample(V2, np)
+    V2 = [v2 .+ Vc[end] for v2 ∈ V2] 
+    
+    @testset "quad" begin
+        F,V = sweeploft(Vc,V1,V2; face_type=:quad, num_twist=0, close_loop=true)
+        @test length(F) == (nc-1)*np
+        @test isa(F,Vector{QuadFace{Int64}})
+        @test length(V) == nc*np
+        @test isa(V,Vector{Point3{Float64}})
+
+        F,V = sweeploft(Vc,V1,V2; face_type=:quad, num_twist=1, close_loop=false)
+        @test length(F) == (nc-1)*(np-1)
+        @test isa(F,Vector{QuadFace{Int64}})
+        @test length(V) == nc*np
+        @test isa(V,Vector{Point3{Float64}})        
+    end
+
+    @testset "tri_slash" begin
+        F,V = sweeploft(Vc,V1,V2; face_type=:tri_slash, num_twist=0, close_loop=true)
+        @test length(F) == (nc-1)*np*2
+        @test isa(F,Vector{TriangleFace{Int64}})
+        @test length(V) == nc*np
+        @test isa(V,Vector{Point3{Float64}})
+
+        F,V = sweeploft(Vc,V1,V2; face_type=:tri_slash, num_twist=1, close_loop=false)
+        @test length(F) == (nc-1)*(np-1)*2
+        @test isa(F,Vector{TriangleFace{Int64}})
+        @test length(V) == nc*np
+        @test isa(V,Vector{Point3{Float64}})        
+    end
+
+    @testset "tri" begin
+        F,V = sweeploft(Vc,V1,V2; face_type=:tri, num_twist=0, close_loop=true)
+        @test length(F) == (nc-1)*np*2
+        @test isa(F,Vector{TriangleFace{Int64}})
+        @test length(V) == nc*np
+        @test isa(V,Vector{Point3{Float64}})
+
+        F,V = sweeploft(Vc,V1,V2; face_type=:tri, num_twist=1, close_loop=false)
+        @test length(F) == (nc-1)*(np-1)*2
+        @test isa(F,Vector{TriangleFace{Int64}})
+        @test length(V) == nc*np
+        @test isa(V,Vector{Point3{Float64}})        
+    end
+end