stereo: adjusted parametric law and initial point distribution

oceanmesh/mesh_generator.py

@@ -698,21 +698,37 @@ 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 +.5) / 180 * np.pi
+    y7 = ((90 - lat) * 128+.25) / 180 * np.pi
+    y8 = ((90 - lat) * 252+.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]])
+        bbox = np.array([[-180, 180], [-89, 90]])
     p = np.mgrid[
-        tuple(slice(min, max + min_edge_length, min_edge_length) for min, max in bbox)
+        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