From 04d88639f6120154b9eb27992001dfd94f5eca05 Mon Sep 17 00:00:00 2001
From: tomsail <saillour.thomas@gmail.com>
Date: Thu, 4 Apr 2024 14:44:18 +0200
Subject: [PATCH] added parametric law for stereo distortion (north pole)

---
 oceanmesh/geodata.py        | 6 +++---
 oceanmesh/idw.py            | 1 +
 oceanmesh/mesh_generator.py | 8 +++++++-
 3 files changed, 11 insertions(+), 4 deletions(-)

diff --git a/oceanmesh/geodata.py b/oceanmesh/geodata.py
index d7268d6..9f0a514 100644
--- a/oceanmesh/geodata.py
+++ b/oceanmesh/geodata.py
@@ -794,9 +794,9 @@ def __init__(self, dem, crs="EPSG:4326", bbox=None, extrapolate=False):
                     topobathy = np.transpose(topobathy, (1, 0))
             # Ensure its a floating point array
             topobathy = topobathy.astype(np.float64)
-            topobathy[
-                topobathy == nodata_value
-            ] = np.nan  # set the no-data value to nan
+            topobathy[topobathy == nodata_value] = (
+                np.nan
+            )  # set the no-data value to nan
         elif not dem.exists():
             raise FileNotFoundError(f"File {dem} could not be located.")
 
diff --git a/oceanmesh/idw.py b/oceanmesh/idw.py
index f71b8fa..a165118 100644
--- a/oceanmesh/idw.py
+++ b/oceanmesh/idw.py
@@ -1,6 +1,7 @@
 """ invdisttree.py: inverse-distance-weighted interpolation using KDTree
     fast, solid, local
 """
+
 from __future__ import division
 
 import numpy as np
diff --git a/oceanmesh/mesh_generator.py b/oceanmesh/mesh_generator.py
index 06ae670..fca4a27 100644
--- a/oceanmesh/mesh_generator.py
+++ b/oceanmesh/mesh_generator.py
@@ -581,7 +581,7 @@ def _compute_forces(p, t, fh, min_edge_length, L0mult, opts):
         p1 = p[bars].sum(1) / 2
         x, y = to_lat_lon(p1[:, 0], p1[:, 1])
         p2 = np.asarray([x, y]).T
-        hbars = fh(p2) * _stereo_distortion_dist(y)
+        hbars = fh(p2) * _stereo_distortion_dist(y) * _parametric(y)
     else:
         hbars = fh(p[bars].sum(1) / 2)
     L0 = hbars * L0mult * (np.nanmedian(L) / np.nanmedian(hbars))
@@ -696,6 +696,12 @@ def _stereo_distortion_dist(lat):
     return res
 
 
+def _parametric(lat):
+    ones = np.ones(lat.shape)
+    res = ((90 - lat) * 2 + 18) / 180 * np.pi
+    return np.minimum(res, ones)
+
+
 def _generate_initial_points(min_edge_length, geps, bbox, fh, fd, pfix, stereo=False):
     """Create initial distribution in bounding box (equilateral triangles)"""
     if stereo: