Skip to content

Commit

Permalink
Docs for spherical distance
Browse files Browse the repository at this point in the history
  • Loading branch information
milankl committed Oct 18, 2024
1 parent 6e8baae commit 1587415
Showing 1 changed file with 41 additions and 0 deletions.
41 changes: 41 additions & 0 deletions docs/src/orography.md
Original file line number Diff line number Diff line change
Expand Up @@ -171,6 +171,47 @@ as illustrated in the spectral representation or the surface geopotential
here. This is because scalar fields do not use that last degree,
see [One more degree for spectral fields](@ref).

## Spherical distance

In the example above we have defined the "Super Hawaii orography" as

```julia
orography=(λ,φ) -> H*exp((--λ₀)^2 --φ₀)^2)/2σ^2)
```

note however, that this approximates the distance on the sphere with
Cartesian coordinates which is here not too bad as we are not too far north
(where longitude distance become considerably shorter) and also as
we are far away from the prime meridian. If ``\lambda_0 = 0`` in the example
above then calculating ``\lambda - \lambda_0`` for ``\lambda = 359˚E``
is really far even though ``\lambda`` is actually relatively close to
the prime meridian. To avoid this problem, we define a function called
`spherical_distance` (inputs in degrees, output in meters) to actually
calculate the [great-circle distance](https://en.wikipedia.org/wiki/Great-circle_distance)
or spherical distance. Compare

```@example orography
λ₀ = 0 # move "Super Hawaii" mountain onto the prime meridian
set!(model, orography=(λ,φ) -> H*exp((-(λ-λ₀)^2 - (φ-φ₀)^2)/2σ^2))
heatmap(model.orography.orography, title="Mountain on prime meridian, cartesian coordinates")
save("mountain_cartesian.png", ans) # hide
nothing # hide
```
![Mountain in cartesian coordinates](mountain_cartesian.png)

which clearly shows that the mountain is only on the Eastern hemisphere
(the signal on the western hemisphere is due to the interpolation during plotting)
-- probably not what you wanted. Rewrite this as

```@example orography
λ₀ = 0 # move "Super Hawaii" mountain onto the prime meridian
set!(model, orography=(λ,φ) -> H*exp((-spherical_distance((λ,φ), (λ₀,φ₀), radius=1)^2/2σ^2)
heatmap(model.orography.orography, title="Mountain on prime meridian, spherical distance")
save("mountain_spherical.png", ans) # hide
nothing # hide
```
![Mountain in spherical coordinates](mountain_spherical.png)

## Defining a new orography type

You can also define a new orography like we defined `ZonalRidge` or `EarthOrography`.
Expand Down

0 comments on commit 1587415

Please sign in to comment.