diff --git a/src/02_acoustics/03_propagation/tracing/02_beam_pressure.jl b/src/02_acoustics/03_propagation/tracing/02_beam_pressure.jl
index 612f20f..6e5212d 100644
--- a/src/02_acoustics/03_propagation/tracing/02_beam_pressure.jl
+++ b/src/02_acoustics/03_propagation/tracing/02_beam_pressure.jl
@@ -12,13 +12,13 @@ function beam_pressure(::ModelName{:Gaussian},
     c::Union{<:Real, <:AbstractArray{<:Real}},
     r::Union{<:Real, <:AbstractArray{<:Real}},
     A::Union{<:Real, <:AbstractArray{<:Real}},
-    φ::Union{<:Real, <:AbstractArray{<:Real}},
+    ϕ::Union{<:Real, <:AbstractArray{<:Real}},
     τ::Union{<:Real, <:AbstractArray{<:Real}},
     p::Union{<:Real, <:AbstractArray{<:Real}},
     q::Union{<:Real, <:AbstractArray{<:Real}},
     θ::Union{<:Real, <:AbstractArray{<:Real}}
 )
-    return @. complex(A, φ) * sqrt(
+    return @. complex(A, ϕ) * sqrt(
         c / (r * q) * q₀ * f * cos(θ₀)
     ) * cispi(
         -2f * (
@@ -37,7 +37,7 @@ function beam_pressure(model::Model, s::ArcLengthType, n::NormalDisplacementType
     c::Function,
     r::Function,
     A::Function,
-    φ::Function,
+    ϕ::Function,
     τ::Function,
     p::Function,
     q::Function,
@@ -50,7 +50,7 @@ function beam_pressure(model::Model, s::ArcLengthType, n::NormalDisplacementType
     c′ = c(s)
     r′ = r(s)
     A′ = A(s)
-    φ′ = φ(s)
+    φ′ = ϕ(s)
     p′ = p(s)
     q′ = q(s)
     τ′ = τ(s)
@@ -58,6 +58,6 @@ function beam_pressure(model::Model, s::ArcLengthType, n::NormalDisplacementType
     return @. beam_pressure(model, 
         n, f,
         c(0), θ(0), q(0),
-        c(s), r(s), A(s), φ(s), τ(s), p(s), q(s), θ(s)
+        c(s), r(s), A(s), ϕ(s), τ(s), p(s), q(s), θ(s)
     )
 end
diff --git a/src/02_acoustics/03_propagation/tracing/03_reflection.jl b/src/02_acoustics/03_propagation/tracing/03_reflection.jl
index 1a57f15..a75ce57 100644
--- a/src/02_acoustics/03_propagation/tracing/03_reflection.jl
+++ b/src/02_acoustics/03_propagation/tracing/03_reflection.jl
@@ -11,7 +11,7 @@ function reflect!(ntg::ODEIntegrator, vars::NamedTuple, pars::NamedTuple, ctx::T
     ξ = ntg.u[vars.ξ]
     ζ = ntg.u[vars.ζ]
     A = ntg.u[vars.A]
-    φ = ntg.u[vars.φ]
+    ϕ = ntg.u[vars.ϕ]
     p = ntg.u[vars.p_re] + im * ntg.u[vars.p_im]
     q = ntg.u[vars.q_re] + im * ntg.u[vars.q_im]
 
@@ -57,7 +57,7 @@ function reflect!(ntg::ODEIntegrator, vars::NamedTuple, pars::NamedTuple, ctx::T
     R = R_ntf(r, (θ_inc_deg - θ_bnd_deg) |> abs)
 
     ntg.u[[vars.ξ, vars.ζ]] .= cossind(θ_rfl_deg) ./ c
-    ntg.u[[vars.A, vars.φ]] .= [A * abs(R), φ + angle(R)]
+    ntg.u[[vars.A, vars.ϕ]] .= [A * abs(R), ϕ + angle(R)]
     ntg.u[[vars.p_re, vars.p_im]] .= reim(p + q * N)
 
     nothing
diff --git a/src/02_acoustics/03_propagation/tracing/04_beam.jl b/src/02_acoustics/03_propagation/tracing/04_beam.jl
index aecc52c..29f31b0 100644
--- a/src/02_acoustics/03_propagation/tracing/04_beam.jl
+++ b/src/02_acoustics/03_propagation/tracing/04_beam.jl
@@ -30,7 +30,7 @@ function AcousticTracingODESystem(
         θ(s)::Real = θ₀
         c(s)::Real # = c_ocn(r₀, z₀)
         A(s)::Real = 1.0
-        φ(s)::Real = 0
+        ϕ(s)::Real = 0
         τ(s)::Real = 0.0
 		r(s)::Real = r₀
 		z(s)::Real = z₀
@@ -56,7 +56,7 @@ function AcousticTracingODESystem(
         θ ~ atan(ζ, ξ)
         c ~ c_ocn(r, z)
         Ds(A) ~ 0
-        Ds(φ) ~ 0
+        Ds(ϕ) ~ 0
         Ds(τ) ~ 1 / c_ocn(r, z)
         Ds(r) ~ c_ocn(r, z) * ξ
         Ds(z) ~ c_ocn(r, z) * ζ
@@ -75,7 +75,7 @@ function AcousticTracingODESystem(
 
     surface_reflection = [z ~ z_ati(r)] => (
         reflect!,
-        [r, z, ξ, ζ, A, φ, p_re, p_im, q_re, q_im],
+        [r, z, ξ, ζ, A, ϕ, p_re, p_im, q_re, q_im],
         [],
         [],
         (:surface, c_ocn, z_ati, R_srf)
@@ -83,7 +83,7 @@ function AcousticTracingODESystem(
 
     bottom_reflection = [z ~ z_bty(r)] => (
         reflect!,
-        [r, z, ξ, ζ, A, φ, p_re, p_im, q_re, q_im],
+        [r, z, ξ, ζ, A, ϕ, p_re, p_im, q_re, q_im],
         [],
         [],
         (:bottom, c_ocn, z_bty, R_btm)
@@ -127,7 +127,7 @@ struct Beam{
     # ξ::Function
     # ζ::Function
     # A::Function
-    # φ::Function
+    # ϕ::Function
     # τ::Function
     # p::Function
     # q::Function
@@ -149,8 +149,8 @@ struct Beam{
         θ(s::AbstractVector{<:Real}) = sol(s, idxs = sys.θ) |> collect
         A(s::Real) = sol(s, idxs = sys.A)
         A(s::AbstractVector{<:Real}) = sol(s, idxs = sys.A) |> collect
-        φ(s::Real) = sol(s, idxs = sys.φ)
-        φ(s::AbstractVector{<:Real}) = sol(s, idxs = sys.φ) |> collect
+        ϕ(s::Real) = sol(s, idxs = sys.ϕ)
+        ϕ(s::AbstractVector{<:Real}) = sol(s, idxs = sys.ϕ) |> collect
         τ(s::Real) = sol(s, idxs = sys.τ)
         τ(s::AbstractVector{<:Real}) = sol(s, idxs = sys.τ) |> collect
 
@@ -166,7 +166,7 @@ struct Beam{
             beam_pressure(model, s, n;
                 (
                     var => getproperty(beam, var)
-                    for var in (:c, :f, :r, :A, :φ, :τ, :p, :q, :θ)
+                    for var in (:c, :f, :r, :A, :ϕ, :τ, :p, :q, :θ)
                 )...
             )
         end