Add default Ω_Earth and R_Earth for FPlane and BetaPlane (#550)

Add default Ω_Earth and R_Earth for FPlane and BetaPlane

Author: ali-ramadhan

src/coriolis.jl

@@ -1,5 +1,12 @@
 using Oceananigans.Operators: ℑxyᶜᶠᵃ, ℑxyᶠᶜᵃ
 
+#####
+##### Physical constants
+#####
+
+const Ω_Earth = 7.292115e-5 # [s⁻¹] https://en.wikipedia.org/wiki/Earth%27s_rotation#Angular_speed
+const R_Earth = 6371.0e3    # Mean radius of the Earth [m] https://en.wikipedia.org/wiki/Earth
+
 #####
 ##### Functions for non-rotating models
 #####
@@ -23,26 +30,32 @@ struct FPlane{FT} <: AbstractRotation
 end
 
 """
-    FPlane([FT=Float64;] f=nothing, rotation_rate=nothing, latitude=nothing)
+    FPlane([FT=Float64;] f=nothing, rotation_rate=Ω_Earth, latitude=nothing)
 
 Returns a parameter object for constant rotation at the angular frequency
 `f/2`, and therefore with background vorticity `f`, around a vertical axis.
 If `f` is not specified, it is calculated from `rotation_rate` and
 `latitude` according to the relation `f = 2*rotation_rate*sind(latitude).
 
+By default, `rotation_rate` is assumed to be Earth's.
+
 Also called `FPlane`, after the "f-plane" approximation for the local effect of
-Earth's rotation in a planar coordinate system tangent to the Earth's surface.
+a planet's rotation in a planar coordinate system tangent to the planet's surface.
 """
-function FPlane(FT::DataType=Float64; f=nothing, rotation_rate=nothing, latitude=nothing)
+function FPlane(FT::DataType=Float64; f=nothing, rotation_rate=Ω_Earth, latitude=nothing)
 
-    if f == nothing && rotation_rate != nothing && latitude != nothing
-        return FPlane{FT}(2rotation_rate*sind(latitude))
-    elseif f != nothing && rotation_rate == nothing && latitude == nothing
+    use_f = !isnothing(f)
+    use_planet_parameters = !isnothing(latitude)
+
+    if !xor(use_f, use_planet_parameters)
+        throw(ArgumentError("Either both keywords rotation_rate and latitude must be " *
+                            "specified, *or* only f must be specified."))
+    end
+
+    if use_f
         return FPlane{FT}(f)
-    else
-        throw(ArgumentError("Either both keywords rotation_rate and
-                             latitude must be specified, *or* only f
-                             must be specified."))
+    elseif use_planet_parameters
+        return FPlane{FT}(2rotation_rate*sind(latitude))
     end
 end
 
@@ -57,8 +70,7 @@ end
 """
     BetaPlane{T} <: AbstractRotation
 
-A parameter object for meridionally increasing Coriolis
-parameter (`f = f₀ + βy`).
+A parameter object for meridionally increasing Coriolis parameter (`f = f₀ + βy`).
 """
 struct BetaPlane{T} <: AbstractRotation
     f₀ :: T
@@ -67,26 +79,28 @@ end
 
 """
     BetaPlane([T=Float64;] f₀=nothing, β=nothing,
-                           rotation_rate=nothing, latitude=nothing, radius=nothing)
+                           rotation_rate=Ω_Earth, latitude=nothing, radius=R_Earth)
 
 A parameter object for meridionally increasing Coriolis parameter (`f = f₀ + βy`).
 
-The user may specify both `f₀` and `β`, or the three parameters
-`rotation_rate`, `latitude`, and `radius` that specify the rotation rate and radius
-of a planet, and the central latitude at which the `β`-plane approximation is to be made.
+The user may specify both `f₀` and `β`, or the three parameters `rotation_rate`,
+`latitude`, and `radius` that specify the rotation rate and radius of a planet, and
+the central latitude at which the `β`-plane approximation is to be made.
+
+By default, the `rotation_rate` and planet `radius` is assumed to be Earth's.
 """
 function BetaPlane(T=Float64; f₀=nothing, β=nothing,
-                              rotation_rate=nothing, latitude=nothing, radius=nothing)
+                              rotation_rate=Ω_Earth, latitude=nothing, radius=R_Earth)
 
-    f_and_β = f₀ != nothing && β != nothing
-    planet_parameters = rotation_rate != nothing && latitude != nothing && radius != nothing
+    use_f_and_β = !isnothing(f₀) && !isnothing(β)
+    use_planet_parameters = !isnothing(latitude)
 
-    (f_and_β && all(Tuple(p === nothing for p in (rotation_rate, latitude, radius)))) ||
-    (planet_parameters && all(Tuple(p === nothing for p in (f₀, β)))) ||
-        throw(ArgumentError("Either both keywords f₀ and β must be specified,
-                            *or* all of rotation_rate, latitude, and radius."))
+    if !xor(use_f_and_β, use_planet_parameters)
+        throw(ArgumentError("Either both keywords f₀ and β must be specified, " *
+                            "*or* all of rotation_rate, latitude, and radius."))
+    end
 
-    if planet_parameters
+    if use_planet_parameters
         f₀ = 2rotation_rate * sind(latitude)
          β = 2rotation_rate * cosd(latitude) / radius
      end

test/test_coriolis.jl

@@ -1,43 +1,42 @@
-function instantiate_fplane_1(T)
-    coriolis = FPlane(T, f=π)
-    return coriolis.f == T(π)
+function instantiate_fplane_1(FT)
+    coriolis = FPlane(FT, f=π)
+    return coriolis.f == FT(π)
 end
 
-function instantiate_fplane_2(T)
-    coriolis = FPlane(T, rotation_rate=2, latitude=30)
-    return coriolis.f == T(2)
+function instantiate_fplane_2(FT)
+    coriolis = FPlane(FT, rotation_rate=2, latitude=30)
+    return coriolis.f == FT(2)
 end
 
-function instantiate_betaplane_1(T)
-    coriolis = BetaPlane(T, f₀=π, β=2π)
-    return coriolis.f₀ == T(π)
+function instantiate_betaplane_1(FT)
+    coriolis = BetaPlane(FT, f₀=π, β=2π)
+    return coriolis.f₀ == FT(π)
 end
 
-function instantiate_betaplane_2(T)
-    coriolis = BetaPlane(T, latitude=70, radius=2π, rotation_rate=3π)
-    return coriolis.f₀ == T(6π * sind(70))
+function instantiate_betaplane_2(FT)
+    coriolis = BetaPlane(FT, latitude=70, radius=2π, rotation_rate=3π)
+    return coriolis.f₀ == FT(6π * sind(70))
 end
 
 @testset "Coriolis" begin
     println("Testing Coriolis...")
 
     @testset "Coriolis" begin
-        for T in float_types
-            @test instantiate_fplane_1(T)
-            @test instantiate_fplane_2(T)
-            @test instantiate_betaplane_1(T)
-            @test instantiate_betaplane_2(T)
+        for FT in float_types
+            @test instantiate_fplane_1(FT)
+            @test instantiate_fplane_2(FT)
+            @test instantiate_betaplane_1(FT)
+            @test instantiate_betaplane_2(FT)
+
             # Test that FPlane throws an ArgumentError
-            @test_throws ArgumentError FPlane(T, rotation_rate=7e-5)
-            @test_throws ArgumentError FPlane(T, latitude=40)
-            @test_throws ArgumentError FPlane(T, f=1, rotation_rate=7e-5)
-            @test_throws ArgumentError FPlane(T, f=1, latitude=40)
-            @test_throws ArgumentError FPlane(T, f=1, rotation_rate=7e-5, latitude=40)
+            @test_throws ArgumentError FPlane(FT, rotation_rate=7e-5)
+            @test_throws ArgumentError FPlane(FT, f=1, latitude=40)
+            @test_throws ArgumentError FPlane(FT, f=1, rotation_rate=7e-5, latitude=40)
+
             # Non-exhaustively test that BetaPlane throws an ArgumentError
-            @test_throws ArgumentError BetaPlane(T, f₀=1)
-            @test_throws ArgumentError BetaPlane(T, β=1)
-            @test_throws ArgumentError BetaPlane(T, latitude=70)
-            @test_throws ArgumentError BetaPlane(T, f₀=1e-4, β=1e-11, latitude=70)
+            @test_throws ArgumentError BetaPlane(FT, f₀=1)
+            @test_throws ArgumentError BetaPlane(FT, β=1)
+            @test_throws ArgumentError BetaPlane(FT, f₀=1e-4, β=1e-11, latitude=70)
         end
     end
 end