Skip to content

Commit

Permalink
Add 'asinclamp' helper function
Browse files Browse the repository at this point in the history
  • Loading branch information
eliascarv committed Oct 28, 2024
1 parent 1b19bf6 commit 29dec85
Show file tree
Hide file tree
Showing 3 changed files with 10 additions and 4 deletions.
2 changes: 1 addition & 1 deletion src/crs/projected/albers.jl
Original file line number Diff line number Diff line change
Expand Up @@ -98,7 +98,7 @@ function backward(::Type{<:Albers{latₒ,lat₁,lat₂,Datum}}, x, y) where {lat
θ = n 0 ? atan(x, ρₒ - y) : atan(-x, y - ρₒ)
ρ′ = sqrt(x^2 + (ρₒ - y)^2)
α′ = (C - (ρ′^2 * n^2)) / n
β′ = asin(α′ / (1 - ((1 - e²) / (2 * e)) * log((1 - e) / (1 + e))))
β′ = asinclamp(α′ / (1 - ((1 - e²) / (2 * e)) * log((1 - e) / (1 + e))))

λ = θ / n
ϕ = auth2geod(β′, e²)
Expand Down
5 changes: 2 additions & 3 deletions src/crs/projected/orthographic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -160,18 +160,17 @@ Base.convert(::Type{Orthographic{Mode,latₒ}}, coords::CRS{Datum}) where {Mode,
# -----------------

function sphericalinv(x, y, λₒ, ϕₒ)
fix(x) = clamp(x, -one(x), one(x))
ρ = hypot(x, y)
if ρ < atol(x)
λₒ, ϕₒ
else
c = asin(fix(ρ))
c = asinclamp)
sinc = sin(c)
cosc = cos(c)
sinϕₒ = sin(ϕₒ)
cosϕₒ = cos(ϕₒ)
λ = atan(x * sinc, ρ * cosϕₒ * cosc - y * sinϕₒ * sinc)
ϕ = asin(fix(cosc * sinϕₒ + (y * sinc * cosϕₒ / ρ)))
ϕ = asinclamp(cosc * sinϕₒ + (y * sinc * cosϕₒ / ρ))
λ, ϕ
end
end
7 changes: 7 additions & 0 deletions src/utils.jl
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,13 @@ function atanpos(y, x)
ifelse zero(α), α, α + oftype(α, 2π))
end

"""
asinclamp(x)
Adjust `x` to the range `[-1,1]` to handle floating-point errors.
"""
asinclamp(x) = asin(clamp(x, -one(x), one(x)))

"""
checklat(lat)
Expand Down

0 comments on commit 29dec85

Please sign in to comment.