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..8e856e01 100644
--- a/src/filament.jl
+++ b/src/filament.jl
@@ -59,27 +59,25 @@ 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)"
@@ -87,7 +85,7 @@ function velocity_3D_bound_vortex!(
         # Project onto core radius
         cross3!(r2Xr0, r2, r0)
         nr0sq = nr0 * nr0
-        nr2Xr0 = norm3(r2Xr0)
+        nr2Xr0 = smooth_norm3(r2Xr0)
         d_r1_r0 = dot3(r1, r0)
         d_r2_r0 = dot3(r2, r0)
         @inbounds for k in 1:3
@@ -98,9 +96,9 @@ function velocity_3D_bound_vortex!(
         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 +154,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 +162,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 +188,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 +198,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,20 +270,18 @@ 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]
@@ -293,14 +289,14 @@ function velocity_3D_trailing_vortex_semiinfinite!(
         @inbounds for k in 1:3
             cross_tmp[k] = r1[k]/nr1 - Vf[k]
         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] +
                          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 +320,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..25d05727 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!(
diff --git a/test/filament/test_semi_infinite_filament.jl b/test/filament/test_semi_infinite_filament.jl
index ab870fb5..fcf87c62 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!(