Implement segmented ridges with flexible positioning

[adbayonao]

I added a new feature to the compute_thermal_structure function by creating an additional function that allows for the creation of segmented ridges. Previously, only a single ridge segment could be placed at the boundary of the box. Now, it is possible to position the ridge at any location within the domain, and the ridge can consist of multiple segments. Each segment has its own thermal profile.

The new functionality works as follows: First, the coordinates of the ridge segments are used, where each segment has two parts: the start and end coordinates. Using these coordinates, delimiters are created, where the end of one segment aligns with the start of the next. These delimiters divide the domain into distinct regions, with each region corresponding to a specific ridge segment. To assign coordinates to these regions, the function determines which side of each delimiter a given point lies on, considering the direction of the delimiter. Once the domain is divided, the perpendicular distance from each point to the corresponding ridge segment is calculated, which in turn generates the thermal profile for that region. Each region is influenced by only one ridge segment.

This new functionality was implemented using multiple dispatch. Depending on the parameters provided by the user, either the original or the new version of the function is executed. This allows the function to handle both single ridge segment configurations and multiple segment configurations, providing greater flexibility in modeling ridge structures.

Additionally, I modified the add_box function to integrate this new feature, enabling the inclusion of segmented ridges within the box configuration.

[boriskaus]

thanks for your contribution - can you make the title of this PR a bit clearer so it reflects what this is about?

Also, if it is not yet ready to be reviewed you can mark it as draft (upper right corner)

[boriskaus]

if you want this routine to be accessible in the GMG package, this file needs two be loaded in one of the other files with include("Polygon3D.jl")

[adbayonao]

I am going to review this.

[boriskaus]

Note that we have a function add_polygon! already with test defined in test/test_setup_geometry.jl.

Which functionality does this PR add that is not yet possible with the existing routine?
And is it perhaps possible to extend the existing routine (and test file)?

[adbayonao]

Note that we have a function add_polygon! already with test defined in test/test_setup_geometry.jl.

Which functionality does this PR add that is not yet possible with the existing routine? And is it perhaps possible to extend the existing routine (and test file)?

I am going to try to extend the existing routine to add the new functionality

[adbayonao]

thanks for your contribution - can you make the title of this PR a bit clearer so it reflects what this is about?

Also, if it is not yet ready to be reviewed you can mark it as draft (upper right corner)

Ok, thank you. I changed the title.

[github-actions]

Your PR requires formatting changes to meet the project's style guidelines.
Please consider running Runic (git runic main) to apply these changes.

Click here to view the suggested changes.

diff --git a/ext/GMT_utils.jl b/ext/GMT_utils.jl
index 98998b7..14885d9 100644
--- a/ext/GMT_utils.jl
+++ b/ext/GMT_utils.jl
@@ -158,12 +158,12 @@ function import_GeoTIFF(fname::String; fieldname = :layer1, negative = false, is
         if length(size(G.image)) == 3
             data = permutedims(G.image, [2, 1, 3])
         elseif length(size(G.image)) == 2
-            if size(G.image)==(nx,ny)
-              data[:,:,1] = G.image
-            elseif size(G.image)==(ny,nx)
-              data[:,:,1] = G.image'
+            if size(G.image) == (nx, ny)
+                data[:, :, 1] = G.image
+            elseif size(G.image) == (ny, nx)
+                data[:, :, 1] = G.image'
             else
-              error("unknown size; ")
+                error("unknown size; ")
             end
         end
 
diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl
index b33c0a1..d796006 100644
--- a/src/Setup_geometry.jl
+++ b/src/Setup_geometry.jl
@@ -243,13 +243,15 @@ julia> write_paraview(Grid, "Ridge_Thermal_Structure")  # Save model to Paraview
 1-element Vector{String}:
  "Ridge_Thermal_Structure.vts"
 """
-function add_box!(Phase, Temp, Grid::AbstractGeneralGrid;       # required input
-        xlim::Tuple = (20,100), ylim=nothing, zlim::Tuple = (10,80),     # limits of the box
-        Origin=nothing, StrikeAngle=0, DipAngle=0,      # origin & dip/strike
+function add_box!(
+        Phase, Temp, Grid::AbstractGeneralGrid;       # required input
+        xlim::Tuple = (20, 100), ylim = nothing, zlim::Tuple = (10, 80),     # limits of the box
+        Origin = nothing, StrikeAngle = 0, DipAngle = 0,      # origin & dip/strike
         phase = ConstantPhase(1),                       # Sets the phase number(s) in the box
-        T=nothing,                              # Sets the thermal structure (various functions are available)
-        segments=nothing,                       # Allows defining multiple ridge segments  
-        cell=false )                            # if true, Phase and Temp are defined on cell centers
+        T = nothing,                              # Sets the thermal structure (various functions are available)
+        segments = nothing,                       # Allows defining multiple ridge segments
+        cell = false
+    )                            # if true, Phase and Temp are defined on cell centers
 
 
     # Retrieve 3D data arrays for the grid
@@ -1148,31 +1150,31 @@ Note: the thermal age at the mid oceanic ridge is set to 1 year to avoid divisio
     MORside = "left"   # side of box where the MOR is located
     SpreadingVel = 3   # spreading velocity [cm/yr]
     AgeRidge = 0       # Age of the ridge [Myrs
-    maxAge  = 60       # maximum thermal age of plate [Myrs]
+    maxAge = 60       # maximum thermal age of plate [Myrs]
 end
 
 function compute_thermal_structure(Temp, X, Y, Z, Phase, s::SpreadingRateTemp)
-    @unpack Tsurface, Tmantle, Adiabat, MORside, SpreadingVel, AgeRidge, maxAge  = s
-
-    kappa       =   1e-6;
-    SecYear     =   3600*24*365
-    dz          =   Z[end]-Z[1];
-
-    MantleAdiabaticT    =   Tmantle .+ Adiabat*abs.(Z);   # Adiabatic temperature of mantle
-           
-    if MORside=="left"
-        Distance = X .- X[1,1,1];
-    elseif MORside=="right"
-        Distance = X[end,1,1] .- X;
-    elseif MORside=="front"
-        Distance = Y .- Y[1,1,1];
-    elseif MORside=="back"
-        Distance = Y[1,end,1] .- Y;
+    @unpack Tsurface, Tmantle, Adiabat, MORside, SpreadingVel, AgeRidge, maxAge = s
+
+    kappa = 1.0e-6
+    SecYear = 3600 * 24 * 365
+    dz = Z[end] - Z[1]
+
+    MantleAdiabaticT = Tmantle .+ Adiabat * abs.(Z)    # Adiabatic temperature of mantle
+
+    if MORside == "left"
+        Distance = X .- X[1, 1, 1]
+    elseif MORside == "right"
+        Distance = X[end, 1, 1] .- X
+    elseif MORside == "front"
+        Distance = Y .- Y[1, 1, 1]
+    elseif MORside == "back"
+        Distance = Y[1, end, 1] .- Y
 
     else
         error("unknown side")
     end
-        
+
     for i in eachindex(Temp)
         ThermalAge = abs(Distance[i] * 1.0e3 * 1.0e2) / SpreadingVel + AgeRidge * 1.0e6    # Thermal age in years
         if ThermalAge > maxAge * 1.0e6
@@ -1187,7 +1189,7 @@ function compute_thermal_structure(Temp, X, Y, Z, Phase, s::SpreadingRateTemp)
         Temp[i] = (Tsurface .- Tmantle) * erfc((abs.(Z[i]) * 1.0e3) ./ (2 * sqrt(kappa * ThermalAge))) + MantleAdiabaticT[i]
 
     end
-        
+
     return Temp
 end
 
@@ -1218,14 +1220,14 @@ The thermal age is capped at `maxAge` years, and the temperature is adjusted bas
 
 function compute_thermal_structure(Temp, X, Y, Z, Phase, s::SpreadingRateTemp, segments::Vector{Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}})
     @unpack Tsurface, Tmantle, Adiabat, SpreadingVel, AgeRidge, maxAge = s
-    kappa = 1e-6;
+    kappa = 1.0e-6
     SecYear = 3600 * 24 * 365
-    dz = Z[end]-Z[1];
+    dz = Z[end] - Z[1]
 
     MantleAdiabaticT = Tmantle .+ Adiabat * abs.(Z)
 
     #Create delimiters
-    delimiters = [(segments[i][2], segments[i + 1][1]) for i in 1:length(segments) - 1] 
+    delimiters = [(segments[i][2], segments[i + 1][1]) for i in 1:(length(segments) - 1)]
 
     for I in eachindex(X)
         px, py, pz = X[I], Y[I], Z[I]
@@ -1241,18 +1243,18 @@ function compute_thermal_structure(Temp, X, Y, Z, Phase, s::SpreadingRateTemp, s
         Distance = perpendicular_distance_to_segment(px, py, x1, y1, x2, y2)
 
         # Calculate thermal age
-        ThermalAge = abs(Distance * 1e5) / SpreadingVel + AgeRidge * 1e6  # Thermal age in years
-        if ThermalAge > maxAge * 1e6
-            ThermalAge = maxAge * 1e6
+        ThermalAge = abs(Distance * 1.0e5) / SpreadingVel + AgeRidge * 1.0e6  # Thermal age in years
+        if ThermalAge > maxAge * 1.0e6
+            ThermalAge = maxAge * 1.0e6
         end
 
         ThermalAge = ThermalAge * SecYear  # Convert to seconds
         if ThermalAge == 0
-            ThermalAge = 1e-6  # Avoid zero
+            ThermalAge = 1.0e-6  # Avoid zero
         end
 
         # Calculate temperature
-        Temp[I] = (Tsurface - Tmantle) * erfc(abs(pz) * 1e3 / (2 * sqrt(kappa * ThermalAge))) + MantleAdiabaticT[I]
+        Temp[I] = (Tsurface - Tmantle) * erfc(abs(pz) * 1.0e3 / (2 * sqrt(kappa * ThermalAge))) + MantleAdiabaticT[I]
     end
 
     return Temp
@@ -1271,7 +1273,7 @@ end
 # Function to determine the side of a point with respect to a line (adjusted for segment direction)
 function side_of_line(x, y, x1, y1, x2, y2, direction)
     side = (x2 - x1) * (y - y1) - (y2 - y1) * (x - x1)
-    direction == :left ? side > 0 : side < 0
+    return direction == :left ? side > 0 : side < 0
 end
 
 # Function to determine in which region a point lies (based on delimiters)
diff --git a/test/runtests.jl b/test/runtests.jl
index 91cb814..0a4b346 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -75,7 +75,7 @@ using Test
     end
 
     @testset "Ridge Thermal Structure Tests" begin
-    	include("test_ridge_segments.jl")
+        include("test_ridge_segments.jl")
     end
 
     @testset "Waterflow" begin
diff --git a/test/test_ridge_segments.jl b/test/test_ridge_segments.jl
index 36b530a..d394c38 100644
--- a/test/test_ridge_segments.jl
+++ b/test/test_ridge_segments.jl
@@ -4,9 +4,9 @@ using Test
 @testset "Ridge Thermal Structure Tests" begin
     # Grid parameters
     nx, ny, nz = 512, 512, 128
-    x = range(-1000, 1000, length=nx)
-    y = range(-1000, 1000, length=ny)
-    z = range(-660, 0, length=nz)
+    x = range(-1000, 1000, length = nx)
+    y = range(-1000, 1000, length = ny)
+    z = range(-660, 0, length = nz)
     Grid = CartData(xyz_grid(x, y, z))
 
     # Phases and temperature
@@ -17,12 +17,12 @@ using Test
     segments = [
         ((-500.0, -1000.0), (-500.0, 0.0)),  # Segment 1
         ((-250.0, 0.0), (-250.0, 200.0)),    # Segment 2
-        ((-750.0, 200.0), (-750.0, 1000.0))  # Segment 3
+        ((-750.0, 200.0), (-750.0, 1000.0)),  # Segment 3
     ]
 
     # Lithospheric phases
-    lith = LithosphericPhases(Layers=[15 55], Phases=[1 2], Tlab=1250)
-    add_box!(Phases, Temp, Grid; xlim=(-1000.0,0.0), ylim=(-500.0, 500.0), zlim=(-80.0, 0.0), phase = lith, T=SpreadingRateTemp(SpreadingVel=3), segments=segments)
+    lith = LithosphericPhases(Layers = [15 55], Phases = [1 2], Tlab = 1250)
+    add_box!(Phases, Temp, Grid; xlim = (-1000.0, 0.0), ylim = (-500.0, 500.0), zlim = (-80.0, 0.0), phase = lith, T = SpreadingRateTemp(SpreadingVel = 3), segments = segments)
 
     # Add and save results
     Grid = addfield(Grid, (; Phases, Temp))
@@ -30,10 +30,10 @@ using Test
 
     @test minimum(Temp) >= 0.0  # Minimum temperature
     @test maximum(Temp) <= 1350.0  # Maximum temperature
-    @test all(≥(0),    Temp) # Ensure no negative temperatures
+    @test all(≥(0), Temp) # Ensure no negative temperatures
     @test all(≤(1350), Temp) # Ensure no temperatures above max
 
     # Check if phases are correctly assigned in expected regions
     @test first(Phases) == 2  # Example: Verify a point's phase
-    @test last(Phases)  == 2  # Example: Verify another point's phase
+    @test last(Phases) == 2  # Example: Verify another point's phase
 end