Skip to content

Commit

Permalink
solve edge problem
Browse files Browse the repository at this point in the history
  • Loading branch information
Liu committed Jan 3, 2024
1 parent 6034967 commit 2ef7bf3
Showing 1 changed file with 6 additions and 3 deletions.
9 changes: 6 additions & 3 deletions src/emperical_denudation.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,10 +10,13 @@ function calculate_slope(elevation::Matrix{Float64}, cellsize::Float64)
nrows, ncols = size(elevation)
slope = similar(elevation)

padded_elevation = zeros(Float64, nrows + 2, ncols + 2)
padded_elevation[2:end-1, 2:end-1] = elevation

for i in 2:nrows-1
for j in 2:ncols-1
dzdx = ((elevation[i - 1, j + 1] + 2*elevation[i, j + 1] + elevation[i + 1, j + 1]) - (elevation[i - 1, j - 1] + 2*elevation[i, j - 1] + elevation[i + 1, j - 1])) ./ (8 * cellsize)
dzdy = ((elevation[i + 1, j - 1] + 2*elevation[i + 1, j] + elevation[i + 1, j + 1]) - (elevation[i - 1, j - 1] + 2*elevation[i - 1, j] + elevation[i - 1, j + 1]))/ (8 * cellsize)
dzdx = ((padded_elevation[i - 1, j + 1] + 2*padded_elevation[i, j + 1] + padded_elevation[i + 1, j + 1]) - (padded_elevation[i - 1, j - 1] + 2*padded_elevation[i, j - 1] + padded_elevation[i + 1, j - 1])) ./ (8 * cellsize)
dzdy = ((padded_elevation[i + 1, j - 1] + 2*padded_elevation[i + 1, j] + padded_elevation[i + 1, j + 1]) - (padded_elevation[i - 1, j - 1] + 2*padded_elevation[i - 1, j] + padded_elevation[i - 1, j + 1]))/ (8 * cellsize)
slope[i, j] = atan.(sqrt.(dzdx.^2 + dzdy.^2)) * (180 / π)
end
end
Expand All @@ -28,7 +31,7 @@ function calculate_D(precip::Float64, elevation::Matrix{Float64}, cellsize::Floa
# function format
for i in 2:nrows-1
for j in 2:ncols-1
D[i,j] = (9.1363 ./ (1 .+ exp.(-0.008519.*(precip .- 580.51)))) .* (9.0156 ./ (1 .+ exp.(-0.1245.*(slope[i,j] .- 4.91086))))
D[i,j] = (9.1363 ./ (1 .+ exp.(-0.008519.*(precip .- 580.51)))) .* (9.0156 ./ (1 .+ exp.(-0.1245.*(slope[i,j] .- 4.91086)))) # using logistic function
end
end
return D./1000 #m/kyr
Expand Down

1 comment on commit 2ef7bf3

@EmiliaJarochowska
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Currently these are magic numbers, e.g. 9.1363, -0.008519 etc. They should be at least introduced as named parameters, e.g. slope coefficient etc. and their values should be set for a given run. This way users can at least understand where these numbers come from and, in the future, they can be tweaked if we want to use other erosion parameters, e.g. because this model does not yield realistic results.
An even better approach would be to introduce a function describing the relationship so a whole model can be substituted.
A tiny detail is that single-letter names likeD are very cryptic, better would be to call it denudation.

Please sign in to comment.