2023-10-09 16:29:25 (UTC-03:00)

Marcel Rodrigues <marcelgmr@gmail.com>

map: add Eckert IV projection

diff --git a/lib/anim/map.lua b/lib/anim/map.lua
index 048b442..e69d2ba 100644
--- a/lib/anim/map.lua
+++ b/lib/anim/map.lua
@@ -132,6 +132,48 @@ function Mercator:model()
     return {type="sphere", r=self.r}
 end
 
+-- Eckert IV Projection for the Spherical Earth.
+
+local EckertIV = {}
+EckertIV.__index = EckertIV
+
+local e4_cx = 2 / sqrt(4 * pi + pi * pi)
+local e4_cy = 2 * sqrt(pi / (4 + pi))
+local e4_c = 2 + pi / 2
+
+function EckertIV:map(lon, lat)
+    lon, lat = rad(lon), rad(lat)
+    local theta = lat / 2
+    local sinlat = sin(lat)
+    local delta = 1
+    local costheta, sintheta
+    while abs(delta) > 1e-6 do
+        costheta = cos(theta)
+        sintheta = sin(theta)
+        local num = theta + sintheta * costheta + 2 * sintheta - e4_c * sinlat
+        local den = 2 * costheta * (1 + costheta)
+        delta = - num / den
+        theta = theta + delta
+    end
+    local x = e4_cx * self.r * corangle(lon - self.lon) * (1 + costheta)
+    local y = e4_cy * self.r * sintheta
+    return x, y
+end
+
+function EckertIV:inv(x, y)
+    local theta = asin(y / (e4_cy * self.r))
+    local costheta = cos(theta)
+    local sintheta = sin(theta)
+    local lon = self.lon + x / (e4_cx * self.r * (1 + costheta))
+    local lat = asin((theta + sintheta * costheta + 2 * sintheta) / e4_c)
+    lon, lat = deg(lon), deg(lat)
+    return lon, lat
+end
+
+function EckertIV:model()
+    return {type="sphere", r=self.r}
+end
+
 -- Lambert Azimuthal Equal-Area Projection for the Spherical Earth.
 
 local AzimuthalEqualArea = {}
@@ -174,6 +216,7 @@ end
 
 local projs = {
     Mercator=Mercator,
+    EckertIV=EckertIV,
     AzimuthalEqualArea=AzimuthalEqualArea,
 }