diff --git a/src/equations/compressible_euler_1d.jl b/src/equations/compressible_euler_1d.jl
index a71750ff98..542a5b7039 100644
--- a/src/equations/compressible_euler_1d.jl
+++ b/src/equations/compressible_euler_1d.jl
@@ -412,10 +412,10 @@ See also
     return SVector(f1, f2, f3)
 end
 
-# While `normal_direction` isn't strictly necessary in 1D, certain solvers assume that 
-# the normal component is incorporated into the numerical flux. 
-# 
-# See `flux(u, normal_direction::AbstractVector, equations::AbstractEquations{1})` for a 
+# While `normal_direction` isn't strictly necessary in 1D, certain solvers assume that
+# the normal component is incorporated into the numerical flux.
+#
+# See `flux(u, normal_direction::AbstractVector, equations::AbstractEquations{1})` for a
 # similar implementation.
 @inline function flux_ranocha(u_ll, u_rr, normal_direction::AbstractVector,
                               equations::CompressibleEulerEquations1D)
@@ -943,6 +943,12 @@ end
     return rho
 end
 
+@inline function velocity(u, equations::CompressibleEulerEquations1D)
+    rho = u[1]
+    v1 = u[2] / rho
+    return v1
+end
+
 @inline function pressure(u, equations::CompressibleEulerEquations1D)
     rho, rho_v1, rho_e = u
     p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2) / rho)
diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl
index e964d9c2b2..7884e65d6b 100644
--- a/src/equations/compressible_euler_2d.jl
+++ b/src/equations/compressible_euler_2d.jl
@@ -1957,6 +1957,28 @@ end
     return rho
 end
 
+@inline function velocity(u, equations::CompressibleEulerEquations2D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    return SVector(v1, v2)
+end
+
+@inline function velocity(u, orientation::Int, equations::CompressibleEulerEquations2D)
+    rho = u[1]
+    v = u[orientation + 1] / rho
+    return v
+end
+
+@inline function velocity(u, normal_direction::AbstractVector,
+                          equations::CompressibleEulerEquations2D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    v = v1 * normal_direction[1] + v2 * normal_direction[2]
+    return v
+end
+
 @inline function pressure(u, equations::CompressibleEulerEquations2D)
     rho, rho_v1, rho_v2, rho_e = u
     p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2) / rho)
diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl
index 4f4d4553a8..79ade4e7ae 100644
--- a/src/equations/compressible_euler_3d.jl
+++ b/src/equations/compressible_euler_3d.jl
@@ -1692,6 +1692,30 @@ end
     return rho
 end
 
+@inline function velocity(u, equations::CompressibleEulerEquations3D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    v3 = u[4] / rho
+    return SVector(v1, v2, v3)
+end
+
+@inline function velocity(u, orientation::Int, equations::CompressibleEulerEquations3D)
+    rho = u[1]
+    v = u[orientation + 1] / rho
+    return v
+end
+
+@inline function velocity(u, normal_direction::AbstractVector,
+                          equations::CompressibleEulerEquations3D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    v3 = u[4] / rho
+    v = v1 * normal_direction[1] + v2 * normal_direction[2] + v3 * normal_direction[3]
+    return v
+end
+
 @inline function pressure(u, equations::CompressibleEulerEquations3D)
     rho, rho_v1, rho_v2, rho_v3, rho_e = u
     p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho)
diff --git a/src/equations/compressible_euler_multicomponent_1d.jl b/src/equations/compressible_euler_multicomponent_1d.jl
index da434579f7..ec13b1f84e 100644
--- a/src/equations/compressible_euler_multicomponent_1d.jl
+++ b/src/equations/compressible_euler_multicomponent_1d.jl
@@ -618,4 +618,12 @@ end
     return SVector{ncomponents(equations), real(equations)}(u[i + 2] * v
                                                             for i in eachcomponent(equations))
 end
+
+@inline function velocity(u, equations::CompressibleEulerMulticomponentEquations1D)
+    rho = density(u, equations)
+    v1 = u[1] / rho
+    return v1
+end
+
+
 end # @muladd
diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl
index 2424ad0301..6e2099b2ab 100644
--- a/src/equations/compressible_euler_multicomponent_2d.jl
+++ b/src/equations/compressible_euler_multicomponent_2d.jl
@@ -825,4 +825,27 @@ end
     return SVector{ncomponents(equations), real(equations)}(u[i + 3] * v
                                                             for i in eachcomponent(equations))
 end
+
+@inline function velocity(u, equations::CompressibleEulerMulticomponentEquations2D)
+    rho = density(u, equations)
+    v1 = u[1] / rho
+    v2 = u[2] / rho
+    return SVector(v1, v2)
+end
+
+@inline function velocity(u, orientation::Int, equations::CompressibleEulerMulticomponentEquations2D)
+    rho = density(u, equations)
+    v = u[orientation] / rho
+    return v
+end
+
+@inline function velocity(u, normal_direction::AbstractVector,
+                          equations::CompressibleEulerMulticomponentEquations2D)
+                          rho = density(u, equations)
+    v1 = u[1] / rho
+    v2 = u[2] / rho
+    v = v1 * normal_direction[1] + v2 * normal_direction[2]
+    return v
+end
+
 end # @muladd
diff --git a/src/equations/ideal_glm_mhd_1d.jl b/src/equations/ideal_glm_mhd_1d.jl
index 3253e4410f..cbfe37daa7 100644
--- a/src/equations/ideal_glm_mhd_1d.jl
+++ b/src/equations/ideal_glm_mhd_1d.jl
@@ -326,7 +326,7 @@ function flux_hllc(u_ll, u_rr, orientation::Integer,
         Sdiff = SsR - SsL
 
         # Compute HLL values for vStar, BStar
-        # These correspond to eq. (28) and (30) from the referenced paper 
+        # These correspond to eq. (28) and (30) from the referenced paper
         # and the classic HLL intermediate state given by (2)
         rho_HLL = (SsR * rho_rr - SsL * rho_ll - (f_rr[1] - f_ll[1])) / Sdiff
 
@@ -394,7 +394,7 @@ function flux_hllc(u_ll, u_rr, orientation::Integer,
                          B1_rr * (B1_rr * v1_rr + B2_rr * v2_rr + B3_rr * v3_rr))) /
                        SdiffStar # (23)
 
-            # Classic HLLC update (32)           
+            # Classic HLLC update (32)
             f1 = f_rr[1] + SsR * (densStar - u_rr[1])
             f2 = f_rr[2] + SsR * (mom_1_Star - u_rr[2])
             f3 = f_rr[3] + SsR * (mom_2_Star - u_rr[3])
@@ -555,6 +555,12 @@ end
     return rho
 end
 
+@inline function velocity(u, equations::IdealGlmMhdEquations1D)
+    rho = u[1]
+    v1 = u[2] / rho
+    return v1
+end
+
 @inline function pressure(u, equations::IdealGlmMhdEquations1D)
     rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3 = u
     p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho
diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl
index 7a6d19facd..801f89afab 100644
--- a/src/equations/ideal_glm_mhd_2d.jl
+++ b/src/equations/ideal_glm_mhd_2d.jl
@@ -1119,6 +1119,28 @@ end
     return rho
 end
 
+@inline function velocity(u, equations::IdealGlmMhdEquations2D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    return SVector(v1, v2)
+end
+
+@inline function velocity(u, orientation::Int, equations::IdealGlmMhdEquations2D)
+    rho = u[1]
+    v = u[orientation + 1] / rho
+    return v
+end
+
+@inline function velocity(u, normal_direction::AbstractVector,
+                          equations::IdealGlmMhdEquations2D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    v = v1 * normal_direction[1] + v2 * normal_direction[2]
+    return v
+end
+
 @inline function pressure(u, equations::IdealGlmMhdEquations2D)
     rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u
     p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho
diff --git a/src/equations/ideal_glm_mhd_3d.jl b/src/equations/ideal_glm_mhd_3d.jl
index e922a2e6fd..afbddbb78b 100644
--- a/src/equations/ideal_glm_mhd_3d.jl
+++ b/src/equations/ideal_glm_mhd_3d.jl
@@ -1010,6 +1010,31 @@ end
     return u[1]
 end
 
+@inline function velocity(u, equations::IdealGlmMhdEquations3D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    v3 = u[4] / rho
+    return SVector(v1, v2, v3)
+end
+
+@inline function velocity(u, orientation::Int, equations::IdealGlmMhdEquations3D)
+    rho = u[1]
+    v = u[orientation + 1] / rho
+    return v
+end
+
+@inline function velocity(u, normal_direction::AbstractVector,
+                          equations::IdealGlmMhdEquations3D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    v3 = u[4] / rho
+    v = v1 * normal_direction[1] + v2 * normal_direction[2] + v3 * normal_direction[3]
+    return v
+end
+
+
 @inline function pressure(u, equations::IdealGlmMhdEquations3D)
     rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u
     p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho
diff --git a/src/equations/polytropic_euler_2d.jl b/src/equations/polytropic_euler_2d.jl
index 571d4723f6..e060b389ff 100644
--- a/src/equations/polytropic_euler_2d.jl
+++ b/src/equations/polytropic_euler_2d.jl
@@ -389,6 +389,19 @@ end
     return rho
 end
 
+@inline function velocity(u, equations::PolytropicEulerEquations2D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    return SVector(v1, v2)
+end
+
+@inline function velocity(u, orientation, equations::PolytropicEulerEquations2D)
+    rho = u[1]
+    v = u[orientation + 1] / rho
+    return v
+end
+
 @inline function pressure(u, equations::PolytropicEulerEquations2D)
     rho, rho_v1, rho_v2 = u
     p = equations.kappa * rho^equations.gamma
diff --git a/src/equations/shallow_water_2d.jl b/src/equations/shallow_water_2d.jl
index 4ecaf3b6e1..3f5cc082dd 100644
--- a/src/equations/shallow_water_2d.jl
+++ b/src/equations/shallow_water_2d.jl
@@ -283,7 +283,7 @@ Further details are available in the papers:
   shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry
   [DOI: 10.1016/j.jcp.2017.03.036](https://doi.org/10.1016/j.jcp.2017.03.036)
 - Patrick Ersing, Andrew R. Winters (2023)
-  An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on 
+  An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on
   curvilinear meshes
   [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699)
 """
@@ -840,6 +840,21 @@ end
     return SVector(v1, v2)
 end
 
+@inline function velocity(u, orientation::Int, equations::ShallowWaterEquations2D)
+    h = u[1]
+    v = u[orientation + 1] / h
+    return v
+end
+
+@inline function velocity(u, normal_direction::AbstractVector,
+                          equations::ShallowWaterEquations2D)
+    h = u[1]
+    v1 = u[2] / h
+    v2 = u[3] / h
+    v = v1 * normal_direction[1] + v2 * normal_direction[2]
+    return v
+end
+
 # Convert conservative variables to primitive
 @inline function cons2prim(u, equations::ShallowWaterEquations2D)
     h, _, _, b = u
@@ -903,11 +918,11 @@ end
                        equations::ShallowWaterEquations2D)
 
 Calculate Roe-averaged velocity `v_roe` and wavespeed `c_roe = sqrt{g * h_roe}` depending on direction.
-See for instance equation (62) in 
+See for instance equation (62) in
 - Paul A. Ullrich, Christiane Jablonowski, and Bram van Leer (2010)
   High-order finite-volume methods for the shallow-water equations on the sphere
   [DOI: 10.1016/j.jcp.2010.04.044](https://doi.org/10.1016/j.jcp.2010.04.044)
-Or [this slides](https://faculty.washington.edu/rjl/classes/am574w2011/slides/am574lecture20nup3.pdf), 
+Or [this slides](https://faculty.washington.edu/rjl/classes/am574w2011/slides/am574lecture20nup3.pdf),
 slides 8 and 9.
 """
 @inline function calc_wavespeed_roe(u_ll, u_rr, orientation::Integer,
diff --git a/test/test_type.jl b/test/test_type.jl
index d507cded95..9787bdd956 100644
--- a/test/test_type.jl
+++ b/test/test_type.jl
@@ -167,6 +167,7 @@ isdir(outdir) && rm(outdir, recursive = true)
             @test eltype(@inferred cons2entropy(u, equations)) == RealT
             @test eltype(@inferred entropy2cons(u, equations)) == RealT
             @test typeof(@inferred density(u, equations)) == RealT
+            @test typeof(@inferred velocity(u, equations)) == RealT
             @test typeof(@inferred pressure(u, equations)) == RealT
             @test typeof(@inferred density_pressure(u, equations)) == RealT
             @test typeof(@inferred entropy(cons, equations)) == RealT
@@ -218,6 +219,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                 end
             end
 
+            @test eltype(@inferred velocity(u, normal_direction, equations)) == RealT
             @test eltype(@inferred flux(u, normal_direction, equations)) == RealT
             @test eltype(@inferred flux_shima_etal(u_ll, u_rr, normal_direction, equations)) ==
                   RealT
@@ -257,6 +259,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                   RealT
 
             for orientation in orientations
+                @test eltype(@inferred velocity(u, orientation, equations)) == RealT
                 @test eltype(@inferred flux(u, orientation, equations)) == RealT
                 @test eltype(@inferred flux_shima_etal(u_ll, u_rr, orientation, equations)) ==
                       RealT
@@ -300,6 +303,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                       RealT
             end
 
+            @test eltype(@inferred velocity(u, equations)) == RealT
             @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT
             @test eltype(@inferred cons2prim(u, equations)) == RealT
             @test eltype(@inferred prim2cons(u, equations)) == RealT
@@ -367,6 +371,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                 end
             end
 
+            @test eltype(@inferred velocity(u, normal_direction, equations)) == RealT
             @test eltype(@inferred flux(u, normal_direction, equations)) == RealT
             @test eltype(@inferred flux_shima_etal(u_ll, u_rr, normal_direction, equations)) ==
                   RealT
@@ -392,6 +397,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                                                           equations)) == RealT
 
             for orientation in orientations
+                @test eltype(@inferred velocity(u, orientation, equations)) == RealT
                 @test eltype(@inferred flux(u, orientation, equations)) == RealT
                 @test eltype(@inferred flux_shima_etal(u_ll, u_rr, orientation, equations)) ==
                       RealT
@@ -421,6 +427,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                                                               equations)) == RealT
             end
 
+            @test eltype(@inferred velocity(u, equations)) == RealT
             @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT
             @test eltype(@inferred cons2prim(u, equations)) == RealT
             @test eltype(@inferred prim2cons(u, equations)) == RealT
@@ -456,6 +463,7 @@ isdir(outdir) && rm(outdir, recursive = true)
             @test eltype(@inferred source_terms_convergence_test(u, x, t, equations)) ==
                   RealT
 
+            @test eltype(@inferred velocity(u, equations)) == RealT
             @test eltype(@inferred flux(u, orientation, equations)) == RealT
             @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation,
                                                       equations)) == RealT
@@ -499,11 +507,13 @@ isdir(outdir) && rm(outdir, recursive = true)
             @test eltype(@inferred source_terms_convergence_test(u, x, t, equations)) ==
                   RealT
 
+            @test eltype(@inferred velocity(u, normal_direction, equations)) == RealT
             @test eltype(@inferred flux(u, normal_direction, equations)) == RealT
             @test eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction,
                                                 equations)) == RealT
 
             for orientation in orientations
+                @test eltype(@inferred velocity(u, orientation, equations)) == RealT
                 @test eltype(@inferred flux(u, orientation, equations)) == RealT
                 @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation,
                                                           equations)) == RealT
@@ -515,6 +525,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                       RealT
             end
 
+            @test eltype(@inferred velocity(u, equations)) == RealT
             @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT
             @test eltype(@inferred cons2prim(u, equations)) == RealT
             @test eltype(@inferred prim2cons(u, equations)) == RealT
@@ -1094,6 +1105,7 @@ isdir(outdir) && rm(outdir, recursive = true)
             @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) ==
                   RealT
 
+            @test eltype(@inferred velocity(u, equations)) == RealT
             @test eltype(@inferred flux(u, orientation, equations)) == RealT
             @test eltype(@inferred flux_hllc(u_ll, u_rr, orientation, equations)) == RealT
             @test eltype(@inferred flux_derigs_etal(u_ll, u_rr, orientation, equations)) ==
@@ -1165,6 +1177,7 @@ isdir(outdir) && rm(outdir, recursive = true)
             @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) ==
                   RealT
 
+            @test eltype(@inferred velocity(u, normal_direction, equations)) == RealT
             @test eltype(@inferred flux(u, normal_direction, equations)) == RealT
             @test eltype(@inferred flux_nonconservative_powell(u_ll, u_rr,
                                                                normal_direction_ll,
@@ -1183,6 +1196,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                                                           equations)) == RealT
 
             for orientation in orientations
+                @test eltype(@inferred velocity(u, orientation, equations)) == RealT
                 @test eltype(@inferred flux(u, orientation, equations)) == RealT
                 @test eltype(@inferred flux_nonconservative_powell(u_ll, u_rr, orientation,
                                                                    equations)) == RealT
@@ -1221,6 +1235,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                                                               equations)) == RealT
             end
 
+            @test eltype(@inferred velocity(u, equations)) == RealT
             @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT
             @test eltype(@inferred cons2prim(u, equations)) == RealT
             @test eltype(@inferred prim2cons(u, equations)) == RealT
@@ -1283,6 +1298,7 @@ isdir(outdir) && rm(outdir, recursive = true)
             @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) ==
                   RealT
 
+            @test eltype(@inferred velocity(u, normal_direction, equations)) == RealT
             @test eltype(@inferred flux(u, normal_direction, equations)) == RealT
             @test eltype(@inferred flux_nonconservative_powell(u_ll, u_rr,
                                                                normal_direction_ll,
@@ -1301,6 +1317,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                                                           equations)) == RealT
 
             for orientation in orientations
+                @test eltype(@inferred velocity(u, orientation, equations)) == RealT
                 @test eltype(@inferred flux(u, orientation, equations)) == RealT
                 @test eltype(@inferred flux_nonconservative_powell(u_ll, u_rr, orientation,
                                                                    equations)) == RealT
@@ -1319,6 +1336,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                                                               equations)) == RealT
             end
 
+            @test eltype(@inferred velocity(u, equations)) == RealT
             @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT
             @test eltype(@inferred cons2prim(u, equations)) == RealT
             @test eltype(@inferred prim2cons(u, equations)) == RealT
@@ -1425,6 +1443,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                   RealT
 
             for orientation in orientations
+                @test eltype(@inferred velocity(u, orientation, equations)) == RealT
                 @test eltype(@inferred flux(u, orientation, equations)) == RealT
                 @test eltype(@inferred flux_nonconservative_powell(u_ll, u_rr, orientation,
                                                                    equations)) ==
@@ -2086,6 +2105,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                 @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) ==
                       RealT
 
+                @test eltype(@inferred velocity(u, normal_direction, equations)) == RealT
                 @test eltype(@inferred flux(u, normal_direction, equations)) == RealT
                 @test eltype(@inferred flux_winters_etal(u_ll, u_rr, normal_direction,
                                                          equations)) ==
@@ -2101,6 +2121,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                       RealT
 
                 for orientation in orientations
+                    @test eltype(@inferred velocity(u, orientation, equations)) == RealT
                     @test eltype(@inferred flux(u, orientation, equations)) == RealT
                     @test eltype(@inferred flux_winters_etal(u_ll, u_rr, orientation,
                                                              equations)) ==
@@ -2113,6 +2134,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                           RealT
                 end
 
+                @test eltype(@inferred velocity(u, equations)) == RealT
                 @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT
                 @test eltype(@inferred cons2prim(u, equations)) == RealT
                 @test eltype(@inferred prim2cons(u, equations)) == RealT
@@ -2240,6 +2262,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                                                                 surface_flux_function,
                                                                 equations)) == RealT
 
+            @test eltype(@inferred velocity(u, normal_direction, equations)) == RealT
             @test eltype(@inferred flux(u, normal_direction, equations)) == RealT
             @test eltype(@inferred flux_nonconservative_wintermeyer_etal(u_ll, u_rr,
                                                                          normal_direction_ll,
@@ -2285,6 +2308,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                           RealT
                 end
 
+                @test eltype(@inferred velocity(u, orientation, equations)) == RealT
                 @test eltype(@inferred flux(u, orientation, equations)) == RealT
                 @test eltype(@inferred flux_nonconservative_wintermeyer_etal(u_ll, u_rr,
                                                                              orientation,