stereo: adjusted parametric law and initial point distribution

oceanmesh/mesh_generator.py

@@ -698,21 +698,39 @@ def _stereo_distortion_dist(lat):
 
 def _parametric(lat):
     ones = np.ones(lat.shape)
-    res = ((90 - lat) * 2 + 18) / 180 * np.pi
-    return np.minimum(res, ones)
+    y1 = ((90 - lat) * 2 + 16) / 180 * np.pi
+    y2 = ((90 - lat) * 4 + 8) / 180 * np.pi
+    y3 = ((90 - lat) * 8 + 4) / 180 * np.pi
+    y4 = ((90 - lat) * 16 + 2) / 180 * np.pi
+    y5 = ((90 - lat) * 32 + 1) / 180 * np.pi
+    y6 = ((90 - lat) * 64 + 0.5) / 180 * np.pi
+    y7 = ((90 - lat) * 128 + 0.25) / 180 * np.pi
+    y8 = ((90 - lat) * 252 + 0.125) / 180 * np.pi
+    y = np.minimum(y1, ones)
+    y = np.minimum(y2, y)
+    y = np.minimum(y3, y)
+    y = np.minimum(y4, y)
+    y = np.minimum(y5, y)
+    y = np.minimum(y6, y)
+    y = np.minimum(y7, y)
+    y = np.minimum(y8, y)
+    return y
 
 
 def _generate_initial_points(min_edge_length, geps, bbox, fh, fd, pfix, stereo=False):
     """Create initial distribution in bounding box (equilateral triangles)"""
     if stereo:
-        bbox = np.array([[-180, 180], [-89, 89]])
-    p = np.mgrid[
-        tuple(slice(min, max + min_edge_length, min_edge_length) for min, max in bbox)
-    ].astype(float)
+        bbox = np.array([[-180, 180], [-89, 90]])
+    p = np.mgrid[tuple(slice(min, max, min_edge_length) for min, max in bbox)].astype(
+        float
+    )
     if stereo:
         # for global meshes in stereographic projections,
         # we need to reproject the points from lon/lat to stereo projection
         # then, we need to rectify their coordinates to lat/lon for the sizing function
+        p += (
+            np.random.rand(*p.shape) * min_edge_length / 2
+        )  # randomise the distribution
         p0 = p.reshape(2, -1).T
         x, y = to_stereo(p0[:, 0], p0[:, 1])
         p = np.asarray([x, y]).T