diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl
index 8ca7c003..53d968ef 100644
--- a/src/body_aerodynamics.jl
+++ b/src/body_aerodynamics.jl
@@ -625,7 +625,7 @@ function find_center_of_pressure(
         normal[1] = v1y*v2z - v1z*v2y
         normal[2] = v1z*v2x - v1x*v2z
         normal[3] = v1x*v2y - v1y*v2x
-        normal_norm = norm3(normal)
+        normal_norm = smooth_norm3(normal)
         if normal_norm != 0
             normal[1] /= normal_norm
             normal[2] /= normal_norm
@@ -771,8 +771,8 @@ function calculate_results(
             panel, alpha_dist[i])
         panel_width_array[i] = panel.width
         va_norm = va_norm_array[i]
-        x_norm = norm3(panel.x_airf)
-        z_norm = norm3(panel.z_airf)
+        x_norm = smooth_norm3(panel.x_airf)
+        z_norm = smooth_norm3(panel.z_airf)
         if va_norm == 0.0 || x_norm == 0.0 || z_norm == 0.0
             alpha_geometric[i] = NaN
         else
@@ -824,7 +824,7 @@ function calculate_results(
     total_area > 0.0 || throw(ArgumentError(
         "Total panel area must be positive."))
     reference_speed = sqrt(weighted_speed_sq / total_area)
-    direction_norm = norm3(va_ref_vector)
+    direction_norm = smooth_norm3(va_ref_vector)
     if direction_norm <= 0.0
         va_ref_vector .= (1.0, 0.0, 0.0)
         direction_norm = 1.0
@@ -833,7 +833,7 @@ function calculate_results(
         va_ref_vector[k] = va_ref_vector[k] / direction_norm *
                            reference_speed
     end
-    va_ref_mag = norm3(va_ref_vector)
+    va_ref_mag = smooth_norm3(va_ref_vector)
     va_ref_mag > 0.0 || throw(ArgumentError(
         "Reference freestream magnitude must be positive."))
     va_ref_unit = body_aero.work_vectors[1]
@@ -843,7 +843,7 @@ function calculate_results(
     end
     dir_lift_ref = body_aero.work_vectors[2]
     cross3!(dir_lift_ref, va_ref_vector, spanwise_direction)
-    dir_lift_ref_norm = norm3(dir_lift_ref)
+    dir_lift_ref_norm = smooth_norm3(dir_lift_ref)
     dir_lift_ref_norm > 0.0 || throw(ArgumentError(
         "Reference lift direction is undefined because " *
         "reference flow is parallel to spanwise direction."))
@@ -874,14 +874,14 @@ function calculate_results(
             induced_va_airfoil[k] = c_alpha * panel.x_airf[k] +
                                     s_alpha * panel.z_airf[k]
         end
-        normalize3!(induced_va_airfoil)
+        smooth_normalize3!(induced_va_airfoil)
 
         cross3!(dir_lift_induced_va,
                 induced_va_airfoil, panel.y_airf)
-        normalize3!(dir_lift_induced_va)
+        smooth_normalize3!(dir_lift_induced_va)
         cross3!(dir_drag_induced_va,
                 spanwise_direction, dir_lift_induced_va)
-        normalize3!(dir_drag_induced_va)
+        smooth_normalize3!(dir_drag_induced_va)
 
         q_lift = 0.5 * density * v_a_dist[i]^2
         lift_i = cl_array[i] * q_lift * chord_array[i]
@@ -900,7 +900,7 @@ function calculate_results(
         q_panel = 0.5 * density * va_panel_mag^2
         cross3!(dir_lift_prescribed_va,
                 panel.va, spanwise_direction)
-        normalize3!(dir_lift_prescribed_va)
+        smooth_normalize3!(dir_lift_prescribed_va)
 
         lift_prescribed_va =
             dot3(lift_induced_va, dir_lift_prescribed_va) +
diff --git a/src/filament.jl b/src/filament.jl
index e2c5a772..fce05ff9 100644
--- a/src/filament.jl
+++ b/src/filament.jl
@@ -59,48 +59,51 @@ function velocity_3D_bound_vortex!(
     r2 .= XVP .- filament.x2
 
     # Cut-off radius
-    nr0 = norm3(r0)
+    nr0 = smooth_norm3(r0)
     epsilon = core_radius_fraction * nr0
 
     cross3!(r1Xr2, r1, r2)
     cross3!(r1Xr0, r1, r0)
-    nr1 = norm3(r1)
-    nr2 = norm3(r2)
+    nr1 = smooth_norm3(r1)
+    nr2 = smooth_norm3(r2)
     @inbounds for k in 1:3
         r1r2norm[k] = r1[k]/nr1 - r2[k]/nr2
     end
 
     # Check point location relative to filament
-    nr1Xr0 = norm3(r1Xr0)
+    nr1Xr0 = smooth_norm3(r1Xr0)
     if nr1Xr0 / nr0 > epsilon
-        nr1Xr2 = norm3(r1Xr2)
+        nr1Xr2 = smooth_norm3(r1Xr2)
         coeff = (gamma / (4π)) / (nr1Xr2^2) * dot3(r0, r1r2norm)
         @inbounds for k in 1:3
             vel[k] = coeff * r1Xr2[k]
         end
-    elseif nr1Xr0 / nr0 < 1e-12 * epsilon
-        vel .= 0.0
     else
         @debug "inside core radius"
         @debug "distance from control point to filament: $(nr1Xr0 / nr0)"
 
-        # Project onto core radius
-        cross3!(r2Xr0, r2, r0)
+        # Project onto core cylinder along the radial direction
+        # (perpendicular component of r1 to r0). r1 and r2 share the
+        # same radial component, so one radial unit vector serves both.
         nr0sq = nr0 * nr0
-        nr2Xr0 = norm3(r2Xr0)
         d_r1_r0 = dot3(r1, r0)
         d_r2_r0 = dot3(r2, r0)
+        r_rad = r1Xr0  # reuse buffer; no longer needed below
+        @inbounds for k in 1:3
+            r_rad[k] = r1[k] - d_r1_r0 * r0[k] / nr0sq
+        end
+        nr_rad = smooth_norm3(r_rad)
         @inbounds for k in 1:3
             r1_proj[k] = d_r1_r0 * r0[k] / nr0sq +
-                         epsilon * r1Xr0[k] / nr1Xr0
+                         epsilon * r_rad[k] / nr_rad
             r2_proj[k] = d_r2_r0 * r0[k] / nr0sq +
-                         epsilon * r2Xr0[k] / nr2Xr0
+                         epsilon * r_rad[k] / nr_rad
         end
         cross3!(r1_projXr2_proj, r1_proj, r2_proj)
 
-        nr1pXr2p = norm3(r1_projXr2_proj)
-        nr1_proj = norm3(r1_proj)
-        nr2_proj = norm3(r2_proj)
+        nr1pXr2p = smooth_norm3(r1_projXr2_proj)
+        nr1_proj = smooth_norm3(r1_proj)
+        nr2_proj = smooth_norm3(r2_proj)
         d_sum = 0.0
         @inbounds for k in 1:3
             d_sum += r0[k] * (r1_proj[k]/nr1_proj -
@@ -156,7 +159,7 @@ as implemented in KiteAeroDyn".
     r2 .= XVP .- filament.x2
 
     # Vector perpendicular to core radius
-    nr0 = norm3(r0)
+    nr0 = smooth_norm3(r0)
     nr0sq = nr0 * nr0
     d_r1_r0 = dot3(r1, r0)
     @inbounds for k in 1:3
@@ -164,22 +167,22 @@ as implemented in KiteAeroDyn".
     end
 
     # Cut-off radius
-    epsilon = sqrt(4 * ALPHA0 * NU * norm3(r_perp) / v_a)
+    epsilon = sqrt(4 * ALPHA0 * NU * smooth_norm3(r_perp) / v_a)
 
     cross3!(r1Xr2, r1, r2)
     cross3!(r1Xr0, r1, r0)
     cross3!(r2Xr0, r2, r0)
 
-    nr1 = norm3(r1)
-    nr2 = norm3(r2)
+    nr1 = smooth_norm3(r1)
+    nr2 = smooth_norm3(r2)
     @inbounds for k in 1:3
         normr1r2[k] = r1[k]/nr1 - r2[k]/nr2
     end
 
     # Check point location relative to filament
-    nr1Xr0 = norm3(r1Xr0)
+    nr1Xr0 = smooth_norm3(r1Xr0)
     if nr1Xr0 / nr0 > epsilon
-        nr1Xr2 = norm3(r1Xr2)
+        nr1Xr2 = smooth_norm3(r1Xr2)
         coeff = (gamma / (4π)) / (nr1Xr2^2) * dot3(r0, normr1r2)
         @inbounds for k in 1:3
             vel[k] = coeff * r1Xr2[k]
@@ -190,7 +193,7 @@ as implemented in KiteAeroDyn".
         # Project onto core radius — reuse r_perp, normr1r2
         r1_proj = r_perp
         r2_proj = normr1r2
-        nr2Xr0 = norm3(r2Xr0)
+        nr2Xr0 = smooth_norm3(r2Xr0)
         d_r2_r0 = dot3(r2, r0)
         @inbounds for k in 1:3
             r1_proj[k] = d_r1_r0 * r0[k] / nr0sq +
@@ -200,9 +203,9 @@ as implemented in KiteAeroDyn".
         end
 
         cross3!(r1Xr2, r1_proj, r2_proj)
-        nr1Xr2_val = norm3(r1Xr2)
-        nr1_proj = norm3(r1_proj)
-        nr2_proj = norm3(r2_proj)
+        nr1Xr2_val = smooth_norm3(r1Xr2)
+        nr1_proj = smooth_norm3(r1_proj)
+        nr2_proj = smooth_norm3(r2_proj)
         d_sum = 0.0
         @inbounds for k in 1:3
             d_sum += r0[k] * (r1_proj[k]/nr1_proj -
@@ -272,35 +275,34 @@ function velocity_3D_trailing_vortex_semiinfinite!(
     @inbounds for k in 1:3
         r_perp[k] = d_r1_Vf * Vf[k]
     end
-    epsilon = sqrt(4 * ALPHA0 * NU * norm3(r_perp) / v_a)
+    epsilon = sqrt(4 * ALPHA0 * NU * smooth_norm3(r_perp) / v_a)
 
     cross3!(r1XVf, r1, Vf)
 
-    nr1XVf = norm3(r1XVf)
-    nVf = norm3(Vf)
-    nr1 = norm3(r1)
+    nr1XVf = smooth_norm3(r1XVf)
+    nVf = smooth_norm3(Vf)
+    nr1 = smooth_norm3(r1)
     if nr1XVf / nVf > epsilon
         K = GAMMA / (4π) / (nr1XVf^2) * (1 + d_r1_Vf / nr1)
         @inbounds for k in 1:3
             vel[k] = K * r1XVf[k]
         end
-    elseif nr1XVf / nVf < 1e-12 * epsilon
-        vel .= 0.0
     else
         r1_proj = work_vectors[4]
         cross_tmp = work_vectors[5]
-        # temp = r1/nr1 - Vf
+        # Radial direction: perpendicular component of r1 to Vf.
+        nVfsq = nVf * nVf
         @inbounds for k in 1:3
-            cross_tmp[k] = r1[k]/nr1 - Vf[k]
+            cross_tmp[k] = r1[k] - d_r1_Vf * Vf[k] / nVfsq
         end
-        n_tmp = norm3(cross_tmp)
+        n_tmp = smooth_norm3(cross_tmp)
         @inbounds for k in 1:3
-            r1_proj[k] = d_r1_Vf * Vf[k] +
+            r1_proj[k] = d_r1_Vf * Vf[k] / nVfsq +
                          epsilon * cross_tmp[k] / n_tmp
         end
         cross3!(cross_tmp, r1_proj, Vf)
-        K = GAMMA / (4π) / (norm3(cross_tmp)^2) *
-            (1 + dot3(r1_proj, Vf) / norm3(r1_proj))
+        K = GAMMA / (4π) / (smooth_norm3(cross_tmp)^2) *
+            (1 + dot3(r1_proj, Vf) / smooth_norm3(r1_proj))
         @inbounds for k in 1:3
             vel[k] = K * cross_tmp[k]
         end
@@ -324,10 +326,10 @@ Compute cross product of 3D vectors in-place.
     nothing
 end
 
-@inline norm3(a) = sqrt(a[1]*a[1] + a[2]*a[2] + a[3]*a[3])
+@inline smooth_norm3(a) = sqrt(a[1]*a[1] + a[2]*a[2] + a[3]*a[3] + 1e-18)
 @inline dot3(a, b) = a[1]*b[1] + a[2]*b[2] + a[3]*b[3]
-@inline function normalize3!(v)
-    n = norm3(v)
-    n > 0 && (v[1] /= n; v[2] /= n; v[3] /= n)
+@inline function smooth_normalize3!(v)
+    n = smooth_norm3(v)
+    v[1] /= n; v[2] /= n; v[3] /= n
     nothing
 end
diff --git a/src/panel.jl b/src/panel.jl
index 966b448f..f196a257 100644
--- a/src/panel.jl
+++ b/src/panel.jl
@@ -596,7 +596,7 @@ function calculate_velocity_induced_bound_2D!(
     cross3!(cross_, r0, r3)
 
     # Calculate induced velocity
-    coeff = norm3(r0) / (2π * norm3(cross_)^2)
+    coeff = smooth_norm3(r0) / (2π * smooth_norm3(cross_)^2)
     @inbounds for k in 1:3
         U_2D[k] = cross_[k] * coeff
     end
diff --git a/src/solver.jl b/src/solver.jl
index 4cd4bc79..c9cde9ee 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -354,13 +354,13 @@ function solve!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics, gamma_dist
             dir_iva[k] = c_alpha * panel.x_airf[k] +
                          s_alpha * panel.z_airf[k]
         end
-        normalize3!(dir_iva)
+        smooth_normalize3!(dir_iva)
 
         # Calculate lift and drag directions
         cross3!(dir_lift, dir_iva, panel.y_airf)
-        normalize3!(dir_lift)
+        smooth_normalize3!(dir_lift)
         cross3!(dir_drag, spanwise_direction, dir_lift)
-        normalize3!(dir_drag)
+        smooth_normalize3!(dir_drag)
 
         # Calculate force vectors
         li = lift[i]
@@ -656,7 +656,7 @@ function solve_base!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics, gamma
     nothing
 end
 
-@inline smooth_sqrt(x) = sqrt(x + 1e-30)
+@inline smooth_sqrt(x) = sqrt(x + 1e-18)
 
 @inline function update_gamma_candidate!(
     gamma_out,
diff --git a/test/filament/test_bound_filament.jl b/test/filament/test_bound_filament.jl
index 511af4c3..c5022df4 100644
--- a/test/filament/test_bound_filament.jl
+++ b/test/filament/test_bound_filament.jl
@@ -133,35 +133,35 @@ end
     @testset "Around Core Radius" begin
         filament = create_test_filament()
         delta = 1e-5
-        
+
         points = [
             [0.5, core_radius_fraction - delta, 0.0],
             [0.5, core_radius_fraction, 0.0],
             [0.5, core_radius_fraction + delta, 0.0]
         ]
-        
+
         velocities = [zeros(3) for p in points]
         [
             velocity_3D_bound_vortex!(velocities[i], filament, p, gamma, core_radius_fraction, work_vectors)
             for (i, p) in enumerate(points)
         ]
-        
+
         # Check for NaN and finite values
         @test all(!any(isnan.(v)) for v in velocities)
         @test all(all(isfinite.(v)) for v in velocities)
-        
+
         # Check magnitude is maximum at core radius
         @test norm(velocities[2]) > norm(velocities[1])
         @test norm(velocities[2]) > norm(velocities[3])
-        
+
         # Check continuity around core radius
         @test isapprox(velocities[1], velocities[2], rtol=1e-2)
-        
+
         # Check non-zero velocities
         @test !all(isapprox.(velocities[1], zeros(3), atol=1e-10))
         @test !all(isapprox.(velocities[2], zeros(3), atol=1e-10))
         @test !all(isapprox.(velocities[3], zeros(3), atol=1e-10))
-        
+
         # Check symmetry
         v_neg = zeros(3)
         velocity_3D_bound_vortex!(
@@ -174,4 +174,82 @@ end
         )
         @test isapprox(velocities[2], -v_neg)
     end
+
+    @testset "Velocity is azimuthal (perpendicular to axis and radius)" begin
+        # A real vortex's induced velocity is purely azimuthal:
+        # perpendicular to BOTH the filament axis r0 AND the radial
+        # direction from the axis to the field point. The old Branch 3
+        # projected the field point in the azimuthal direction instead
+        # of the radial direction, producing a velocity with a nonzero
+        # radial component. This test catches that.
+        filament = create_test_filament()
+        r0 = [1.0, 0.0, 0.0]
+
+        # Probe at multiple distances spanning both branches (in-core
+        # and outside-core) and azimuthal positions.
+        for d in (0.25, 0.5, 0.99, 1.0, 1.01, 2.0) .* core_radius_fraction
+            for phi in (0.0, π/4, π/2, π, -π/3)
+                p = [0.5, d * cos(phi), d * sin(phi)]
+                v = zeros(3)
+                velocity_3D_bound_vortex!(
+                    v, filament, p, gamma,
+                    core_radius_fraction, work_vectors)
+
+                # Radial direction at p (perpendicular to axis)
+                r_radial = [0.0, p[2], p[3]]
+
+                # v must be perpendicular to both axis and radial
+                @test isapprox(dot(v, r0), 0.0; atol=1e-10)
+                @test isapprox(dot(v, r_radial), 0.0; atol=1e-8)
+            end
+        end
+    end
+
+    @testset "Solid-body rotation inside core" begin
+        # Inside the core, velocity magnitude scales linearly with
+        # d_perp (Branch 3 = linear ramp from 0 on axis to peak at
+        # boundary). Direction stays constant along a radial line.
+        filament = create_test_filament()
+
+        v_half = zeros(3)
+        v_quarter = zeros(3)
+        velocity_3D_bound_vortex!(
+            v_half, filament,
+            [0.5, 0.5 * core_radius_fraction, 0.0],
+            gamma, core_radius_fraction, work_vectors)
+        velocity_3D_bound_vortex!(
+            v_quarter, filament,
+            [0.5, 0.25 * core_radius_fraction, 0.0],
+            gamma, core_radius_fraction, work_vectors)
+
+        # Magnitudes scale linearly with d_perp
+        @test isapprox(norm(v_quarter), 0.5 * norm(v_half); rtol=1e-3)
+        # Directions are parallel
+        @test isapprox(normalize(v_quarter), normalize(v_half); atol=1e-8)
+    end
+
+    @testset "Branch 1 / Branch 3 agree at core boundary" begin
+        # The two branches are designed to be continuous at d_perp =
+        # epsilon. Probe just inside and just outside the boundary with
+        # enough delta to clearly land in different branches; results
+        # must agree in both magnitude AND direction. Old code returned
+        # orthogonal directions here.
+        filament = create_test_filament()
+        delta = 0.05 * core_radius_fraction  # 5% in/out
+        v_inside = zeros(3)
+        v_outside = zeros(3)
+        velocity_3D_bound_vortex!(
+            v_inside, filament,
+            [0.5, core_radius_fraction - delta, 0.0],
+            gamma, core_radius_fraction, work_vectors)
+        velocity_3D_bound_vortex!(
+            v_outside, filament,
+            [0.5, core_radius_fraction + delta, 0.0],
+            gamma, core_radius_fraction, work_vectors)
+
+        # Magnitudes within a few % of each other (peak is at boundary)
+        @test isapprox(norm(v_inside), norm(v_outside); rtol=0.1)
+        # Directions match — would fail with old azimuthal-projection bug
+        @test isapprox(normalize(v_inside), normalize(v_outside); atol=1e-2)
+    end
 end
\ No newline at end of file
diff --git a/test/filament/test_semi_infinite_filament.jl b/test/filament/test_semi_infinite_filament.jl
index ab870fb5..f041a224 100644
--- a/test/filament/test_semi_infinite_filament.jl
+++ b/test/filament/test_semi_infinite_filament.jl
@@ -73,7 +73,10 @@ end
         ]
         induced_velocity = zeros(3)
 
-        # Filament start point is singular in the current implementation.
+        # Singular geometry: smooth_norm3 floors the cross-product
+        # magnitude, so no division-by-zero produces NaN. Result must
+        # be finite; structural zero cross-products on the axis give
+        # zero velocity.
         velocity_3D_trailing_vortex_semiinfinite!(
             induced_velocity,
             filament,
@@ -83,7 +86,7 @@ end
             filament.vel_mag,
             work_vectors
         )
-        @test all(isnan.(induced_velocity))
+        @test all(isfinite.(induced_velocity))
 
         for point in test_points
             velocity_3D_trailing_vortex_semiinfinite!(
@@ -118,7 +121,7 @@ end
         filament = create_test_filament2()
         vel_pos = zeros(3)
         vel_neg = zeros(3)
-        
+
         velocity_3D_trailing_vortex_semiinfinite!(
             vel_pos,
             filament,
@@ -137,7 +140,55 @@ end
             filament.vel_mag,
             work_vectors
         )
-        
+
         @test isapprox(vel_pos, -vel_neg)
     end
+
+    @testset "Velocity is azimuthal (perpendicular to axis and radius)" begin
+        # Vortex induced velocity is purely azimuthal: perpendicular
+        # to BOTH the filament axis and the radial direction. Old
+        # Branch 3 projected along the azimuthal direction instead of
+        # the radial direction, giving a non-azimuthal velocity.
+        filament = create_test_filament2()
+        direction = filament.direction
+
+        # Lamb-Oseen epsilon at x=0.5 is ~sqrt(4·α₀·ν·0.5) ≈ 6e-3.
+        # Probe across that range.
+        for d in (1e-4, 1e-3, 5e-3, 1e-2, 1e-1)
+            for phi in (0.0, π/4, π/2, π, -π/3)
+                p = [0.5, d * cos(phi), d * sin(phi)]
+                v = zeros(3)
+                velocity_3D_trailing_vortex_semiinfinite!(
+                    v, filament, direction, p, gamma,
+                    filament.vel_mag, work_vectors)
+
+                r_radial = [0.0, p[2], p[3]]
+                @test isapprox(dot(v, direction), 0.0; atol=1e-10)
+                # Looser tol — semi-infinite velocity has an axial
+                # contribution from the (1 + r1·Vf/|r1|) factor.
+                @test isapprox(dot(v, r_radial) / norm(v), 0.0; atol=1e-6)
+            end
+        end
+    end
+
+    @testset "Solid-body rotation inside core" begin
+        # At small d_perp (inside the Lamb-Oseen core), Branch 3 gives
+        # a linear ramp in magnitude with constant direction.
+        filament = create_test_filament2()
+        v_a = filament.vel_mag
+
+        # epsilon at x=0.5 ≈ sqrt(4·α₀·ν·0.5/v_a) ≈ 6.1e-3
+        # Probe well inside the core.
+        d_inside = 1e-4   # ~60× smaller than epsilon — definitely inside
+        v1 = zeros(3); v2 = zeros(3)
+        velocity_3D_trailing_vortex_semiinfinite!(
+            v1, filament, filament.direction,
+            [0.5, d_inside, 0.0], gamma, v_a, work_vectors)
+        velocity_3D_trailing_vortex_semiinfinite!(
+            v2, filament, filament.direction,
+            [0.5, 2 * d_inside, 0.0], gamma, v_a, work_vectors)
+
+        @test isapprox(norm(v2), 2 * norm(v1); rtol=1e-3)
+        @test isapprox(normalize(v2), normalize(v1); atol=1e-8)
+    end
 end