bugfix: Replaced the faulty function for checking if a point is insid…

…e a polygon defined on a sphere with a working approach. Also, fixed predicting the number of stars in the FOV by properly applying the vignetting and extinction correction.

RMS/Math.py

@@ -2,7 +2,7 @@
 
 import numpy as np
 
-
+from RMS.Routines.SphericalPolygonCheck import sphericalPolygonCheck
 
 def lineFunc(x, m, k):
     """ Linear function.
@@ -276,29 +276,24 @@ def sphericalToCartesian(r, theta, phi):
 
 def pointInsideConvexPolygonSphere(points, vertices):
     """
+    LEGACY FUNCTION.
+
     Polygon must be convex
     https://math.stackexchange.com/questions/4012834/checking-that-a-point-is-in-a-spherical-polygon
 
 
     Arguments:
-        points: [array] points with dimension (npoints, 2). The two dimensions are ra and dec
-        vertices: [array] vertices of convex polygon with dimension (nvertices, 2)
+        points: [array] Points with dimension (npoints, 2). The two dimensions are ra and dec in degrees.
+        vertices: [array] Vertices of convex polygon with dimension (nvertices, 2). The two dimensions are 
+            ra and dec in degrees.
         
     Return:
         filter: [array of bool] Array of booleans on whether a given point is inside the polygon on 
             the sphere.
     """
-    # convert ra dec to spherical
-    points = points[:, ::-1]
-    vertices = vertices[:, ::-1]
-    points[:, 0] = 90 - points[:, 0]
-    vertices[:, 0] = 90 - vertices[:, 0]
-    points = np.array(sphericalToCartesian(*np.hstack((np.ones((len(points), 1)), np.radians(points))).T))
-    vertices = np.array(sphericalToCartesian(*np.hstack((np.ones((len(vertices), 1)), np.radians(vertices))).T))
-
-    great_circle_normal = np.cross(vertices, np.roll(vertices, 1, axis=1), axis=0)
-    dot_prod = np.dot(great_circle_normal.T, points)
-    return np.sum(dot_prod < 0, axis=0, dtype=int) == 0  # inside if n . p < 0 for no n
+
+    # Call the new function to check if the points are inside the polygon
+    return sphericalPolygonCheck(vertices, points)
 
 
 ##############################################################################################################
@@ -441,10 +436,106 @@ def testAngSeparationDeg():
     return True
 
 
+def testPointInsideConvexPolygonSphere():
+
+    # # Create a polygon in the sky (has to be convex)
+    # perimiter_polygon = [
+    #     # RA, Dec
+    #     [ 46.00, 79.00],
+    #     [136.00, 84.00],
+    #     [216.55, 75.82],
+    #     [245.85, 61.53],
+    #     [319.86, 56.92],
+    #     [336.00, 70.64],
+    #     [0.00, 77.71],
+    #     ]
+
+    perimiter_polygon = [
+    (0.55, 36.68),
+    (347.46, 43.15),
+    (331.99, 45.77),
+    (315.97, 44.23),
+    (301.11, 38.12),
+    (283.91, 59.51),
+    (239.95, 68.72),
+    (195.93, 59.46),
+    (178.34, 38.02),
+    (163.94, 43.93),
+    (148.00, 45.49),
+    (132.59, 42.86),
+    (119.58, 36.30),
+    (99.46, 56.19),
+    (60.00, 64.18),
+    (20.60, 56.25),
+    (0.55, 36.68),
+]
+
+    perimiter_polygon = np.array(perimiter_polygon)
+
+    # test_points = [
+    #     # RA, Dec
+    #     [47.11, 83.83], # Inside
+    #     [50.23, 74.42], # Outside
+    #     [255.01, 77.33], # Inside
+    #     [185.01, 69.45], # Outside
+    #     [316.75, 64.33] # Outside
+    # ]
+
+    # test_points = np.array(test_points)
+
+    # Sample test points between 0 and 360 degrees in RA and 0 - 90 Declination
+    ra_samples = np.linspace(0, 360, 40)
+    dec_samples = np.linspace(np.min(perimiter_polygon[:, 1]), 90, 20)
+    test_points = np.array(np.meshgrid(ra_samples, dec_samples)).T.reshape(-1, 2)
+
+
+    # Sort the vertices by RA
+    #perimiter_polygon = perimiter_polygon[np.argsort(perimiter_polygon[:, 0])][::-1]
+
+    # # Sort the perimiter by Dec
+    # perimiter_polygon = perimiter_polygon[np.argsort(perimiter_polygon[:, 1])]
+
+
+    # Compute the points that are inside the polygon
+    inside = pointInsideConvexPolygonSphere(test_points, perimiter_polygon)
+
+    # Make the plot in polar coordiantes
+    import matplotlib.pyplot as plt
+
+    fig = plt.figure()
+    ax = fig.add_subplot(111, polar=True)
+
+    # Plot the test points
+    test_points = np.array(test_points)
+    ax.scatter(np.radians(test_points[:, 0]), np.radians(90 - test_points[:, 1]), c='k')
+
+    # Mark the points inside the polygon with an empry green circle
+    ax.scatter(np.radians(test_points[inside, 0]), np.radians(90 - test_points[inside, 1]), edgecolors='g', facecolors='none', s=100, label='Inside')
+
+    # Add the first point to close the polygon
+    ra_dec = np.vstack((perimiter_polygon, perimiter_polygon[0]))
+
+    # Plot the perimeter of the polygon as a continours curve
+    plt.plot(np.radians(ra_dec[:, 0]), np.radians(90 - ra_dec[:, 1]), 'k-')
+
+    # Mark the vertices of the polygon with numbers
+    for i in range(len(perimiter_polygon)):
+        ax.text(np.radians(perimiter_polygon[i, 0]), np.radians(90 - perimiter_polygon[i, 1]), str(i), fontsize=12)
+
+    plt.legend()
+
+
+    ax.set_theta_zero_location('N')
+    ax.set_theta_direction(-1)
+
+    plt.show()
+
+
 def tests():
 
     function_to_test = [["dimHypot", testDimHypot()],
-                         ["angSeparationDeg", testAngSeparationDeg()]]
+                         ["angSeparationDeg", testAngSeparationDeg()],
+                         ["pointInsideConvexPolygonSphere", testPointInsideConvexPolygonSphere()]]
 
     for func_name, func in function_to_test:
         if func: