@@ -243,13 +243,15 @@ julia> write_paraview(Grid, "Ridge_Thermal_Structure") # Save model to Paraview
243
243
1-element Vector{String}:
244
244
"Ridge_Thermal_Structure.vts"
245
245
"""
246
- function add_box! (Phase, Temp, Grid:: AbstractGeneralGrid ; # required input
247
- xlim:: Tuple = (20 ,100 ), ylim= nothing , zlim:: Tuple = (10 ,80 ), # limits of the box
248
- Origin= nothing , StrikeAngle= 0 , DipAngle= 0 , # origin & dip/strike
246
+ function add_box! (
247
+ Phase, Temp, Grid:: AbstractGeneralGrid ; # required input
248
+ xlim:: Tuple = (20 , 100 ), ylim = nothing , zlim:: Tuple = (10 , 80 ), # limits of the box
249
+ Origin = nothing , StrikeAngle = 0 , DipAngle = 0 , # origin & dip/strike
249
250
phase = ConstantPhase (1 ), # Sets the phase number(s) in the box
250
- T= nothing , # Sets the thermal structure (various functions are available)
251
- segments= nothing , # Allows defining multiple ridge segments
252
- cell= false ) # if true, Phase and Temp are defined on cell centers
251
+ T = nothing , # Sets the thermal structure (various functions are available)
252
+ segments = nothing , # Allows defining multiple ridge segments
253
+ cell = false
254
+ ) # if true, Phase and Temp are defined on cell centers
253
255
254
256
255
257
# Retrieve 3D data arrays for the grid
@@ -817,20 +819,22 @@ julia> write_paraview(Grid, "Plate") # Save model to Paraview
817
819
"Plate.vts"
818
820
"""
819
821
820
- function add_plate! (Phase, Temp, Grid:: AbstractGeneralGrid ;
821
- xlim= (), ylim= (), zlim:: Tuple = (0.0 ,0.8 ),
822
- phase = ConstantPhase (1 ),
823
- T= nothing , segments= nothing , cell= false )
822
+ function add_plate! (
823
+ Phase, Temp, Grid:: AbstractGeneralGrid ;
824
+ xlim = (), ylim = (), zlim:: Tuple = (0.0 , 0.8 ),
825
+ phase = ConstantPhase (1 ),
826
+ T = nothing , segments = nothing , cell = false
827
+ )
824
828
825
829
xlim_ = collect (xlim)
826
830
ylim_ = collect (ylim)
827
831
zlim_ = collect (zlim)
828
832
829
- X, Y, Z = coordinate_grids (Grid, cell= cell)
833
+ X, Y, Z = coordinate_grids (Grid, cell = cell)
830
834
ind = zeros (Bool, size (X))
831
835
ind_slice = zeros (Bool, size (X[:, :, 1 ]))
832
836
833
- for k = 1 : size (Z, 3 )
837
+ for k in 1 : size (Z, 3 )
834
838
if zlim_[1 ] <= Z[1 , 1 , k] <= zlim_[2 ]
835
839
inpolygon! (ind_slice, xlim_, ylim_, X[:, :, k], Y[:, :, k])
836
840
@views ind[:, :, k] = ind_slice
@@ -949,7 +953,7 @@ function add_volcano!(
949
953
# @views Temp[ind .== false] .= 0.0
950
954
if ! isempty (ind_flat)
951
955
if ! isnothing (T)
952
- Temp[ind_flat] = compute_thermal_structure (Temp[ind_flat], Grid. x. val[ind], Grid. y. val[ind], depth[ind], Phases[ind_flat], T)
956
+ Temp[ind_flat] = compute_thermal_structure (Temp[ind_flat], Grid. x. val[ind], Grid. y. val[ind], depth[ind], Phases[ind_flat], T)
953
957
end
954
958
end
955
959
@@ -1235,26 +1239,26 @@ Note: the thermal age at the mid oceanic ridge is set to 1 year to avoid divisio
1235
1239
MORside = " left" # side of box where the MOR is located
1236
1240
SpreadingVel = 3 # spreading velocity [cm/yr]
1237
1241
AgeRidge = 0 # Age of the ridge [Myrs]
1238
- maxAge = 60 # maximum thermal age of plate [Myrs]
1242
+ maxAge = 60 # maximum thermal age of plate [Myrs]
1239
1243
end
1240
1244
1241
1245
function compute_thermal_structure (Temp, X, Y, Z, Phase, s:: SpreadingRateTemp )
1242
- @unpack Tsurface, Tmantle, Adiabat, MORside, SpreadingVel, AgeRidge, maxAge = s
1246
+ @unpack Tsurface, Tmantle, Adiabat, MORside, SpreadingVel, AgeRidge, maxAge = s
1243
1247
1244
- kappa = 1e-6 ;
1245
- SecYear = 3600 * 24 * 365
1246
- dz = Z[end ]- Z[1 ];
1248
+ kappa = 1.0e-6
1249
+ SecYear = 3600 * 24 * 365
1250
+ dz = Z[end ] - Z[1 ]
1247
1251
1248
- MantleAdiabaticT = Tmantle .+ Adiabat* abs .(Z); # Adiabatic temperature of mantle
1252
+ MantleAdiabaticT = Tmantle .+ Adiabat * abs .(Z) # Adiabatic temperature of mantle
1249
1253
1250
- if MORside== " left"
1251
- Distance = X .- X[1 ,1 , 1 ];
1252
- elseif MORside== " right"
1253
- Distance = X[end ,1 , 1 ] .- X;
1254
- elseif MORside== " front"
1255
- Distance = Y .- Y[1 ,1 , 1 ];
1256
- elseif MORside== " back"
1257
- Distance = Y[1 ,end ,1 ] .- Y;
1254
+ if MORside == " left"
1255
+ Distance = X .- X[1 , 1 , 1 ]
1256
+ elseif MORside == " right"
1257
+ Distance = X[end , 1 , 1 ] .- X
1258
+ elseif MORside == " front"
1259
+ Distance = Y .- Y[1 , 1 , 1 ]
1260
+ elseif MORside == " back"
1261
+ Distance = Y[1 , end , 1 ] .- Y
1258
1262
1259
1263
else
1260
1264
error (" unknown side" )
@@ -1305,14 +1309,14 @@ The thermal age is capped at `maxAge` years, and the temperature is adjusted bas
1305
1309
1306
1310
function compute_thermal_structure (Temp, X, Y, Z, Phase, s:: SpreadingRateTemp , segments:: Vector{Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}} )
1307
1311
@unpack Tsurface, Tmantle, Adiabat, SpreadingVel, AgeRidge, maxAge = s
1308
- kappa = 1e-6 ;
1312
+ kappa = 1.0e-6
1309
1313
SecYear = 3600 * 24 * 365
1310
- dz = Z[end ]- Z[1 ];
1314
+ dz = Z[end ] - Z[1 ]
1311
1315
1312
1316
MantleAdiabaticT = Tmantle .+ Adiabat * abs .(Z)
1313
1317
1314
1318
# Create delimiters
1315
- delimiters = [(segments[i][2 ], segments[i + 1 ][1 ]) for i in 1 : length (segments) - 1 ]
1319
+ delimiters = [(segments[i][2 ], segments[i + 1 ][1 ]) for i in 1 : ( length (segments) - 1 ) ]
1316
1320
1317
1321
for I in eachindex (X)
1318
1322
px, py, pz = X[I], Y[I], Z[I]
@@ -1328,18 +1332,18 @@ function compute_thermal_structure(Temp, X, Y, Z, Phase, s::SpreadingRateTemp, s
1328
1332
Distance = perpendicular_distance_to_segment (px, py, x1, y1, x2, y2)
1329
1333
1330
1334
# Calculate thermal age
1331
- ThermalAge = abs (Distance * 1e5 ) / SpreadingVel + AgeRidge * 1e6 # Thermal age in years
1332
- if ThermalAge > maxAge * 1e6
1333
- ThermalAge = maxAge * 1e6
1335
+ ThermalAge = abs (Distance * 1.0e5 ) / SpreadingVel + AgeRidge * 1.0e6 # Thermal age in years
1336
+ if ThermalAge > maxAge * 1.0e6
1337
+ ThermalAge = maxAge * 1.0e6
1334
1338
end
1335
1339
1336
1340
ThermalAge = ThermalAge * SecYear # Convert to seconds
1337
1341
if ThermalAge == 0
1338
- ThermalAge = 1e -6 # Avoid zero
1342
+ ThermalAge = 1.0e -6 # Avoid zero
1339
1343
end
1340
1344
1341
1345
# Calculate temperature
1342
- Temp[I] = (Tsurface - Tmantle) * erfc (abs (pz) * 1e3 / (2 * sqrt (kappa * ThermalAge))) + MantleAdiabaticT[I]
1346
+ Temp[I] = (Tsurface - Tmantle) * erfc (abs (pz) * 1.0e3 / (2 * sqrt (kappa * ThermalAge))) + MantleAdiabaticT[I]
1343
1347
end
1344
1348
1345
1349
return Temp
@@ -1358,7 +1362,7 @@ end
1358
1362
# Function to determine the side of a point with respect to a line (adjusted for segment direction)
1359
1363
function side_of_line (x, y, x1, y1, x2, y2, direction)
1360
1364
side = (x2 - x1) * (y - y1) - (y2 - y1) * (x - x1)
1361
- direction == :left ? side > 0 : side < 0
1365
+ return direction == :left ? side > 0 : side < 0
1362
1366
end
1363
1367
1364
1368
# Function to determine in which region a point lies (based on delimiters)
0 commit comments