Changeset 2262 in lmdz_wrf

Timestamp:

Dec 19, 2018, 8:13:23 PM (6 years ago)

Author:

lfita

Message:

Adding polygon-vertex functionalities

File:

• trunk/tools/module_scientific.f90 (modified) (6 diffs)

trunk/tools/module_scientific.f90

@@ -14 +14 @@
 !   to a map of integer polygons filtrered by a minimal area
 ! clean_polygons: Subroutine to clean polygons from non-used paths, polygons only left as path since they are inner path of a hole
+! distanceRK: Function to provide the distance between two points
 ! FindMinimumR_K*: Function returns the location of the minimum in the section between Start and End.
 ! fill3DI_2Dvec: Subroutine to fill a 3D integer matrix from a series of indices from a 2D matrix
@@ -19 +20 @@
 ! gridpoints_InsidePolygon: Subroutine to determine if a series of grid points are inside a polygon 
 !   following ray casting algorithm
+! Heron_area_triangle: Subroutine to compute area of a triangle using Heron's formula
+! intersectfaces: Subroutine to provide if two faces of two polygons intersect
+! intersection_2Dlines: Subroutine to provide the intersection point between two lines on the plane using Cramer's method
+! join_polygon: Subroutine to join two polygons
 ! look_clockwise_borders: Subroutine to look clock-wise for a next point within a collection of borders 
 !   (limits of a region)
@@ -25 +30 @@
 ! path_properties: Subroutine to determine the properties of a path
 ! percentiles_R_K*D: Subroutine to compute the percentiles of a *D R_K array along given set of axis
+! point_inside: Function to determine if a given point is inside a polygon providing its sorted vertices
+! poly_has_point: Function to determine if a polygon has already a given point as one of its vertex
 ! polygon_properties: Subroutine to determine the properties of a polygon (as .TRUE. matrix)
 ! polygons: Subroutine to search the polygons of a border field. FORTRAN based. 1st = 1!
 ! polygons_t: Subroutine to search the polygons of a temporal series of boolean fields. FORTRAN based. 1st = 1!
 ! poly_overlap_tracks: Subroutine to determine tracks of a series of consecutive 2D field with polygons using 
-!   maximum overlaping/coincidence! PrintQuantilesR_K: Subroutine to print the quantiles of values REAL(r_k)
+!   maximum overlaping/coincidence
+! PrintQuantilesR_K: Subroutine to print the quantiles of values REAL(r_k)
 ! poly_overlap_tracks_area: Subroutine to determine tracks of a series of consecutive 2D field with polygons using 
 !   maximum overlaping/coincidence filtrered by a minimal area
@@ -43 +51 @@
 !   and y coordinates
 ! runmean_F1D: Subroutine fo computing the running mean of a given set of float 1D values
+! shoelace_area_polygon: Computing the area of a polygon using sholace formula
+! sort_polygon: Subroutine to sort a polygon using its center as average of the coordinates
 ! SortR_K*: Subroutine receives an array x() r_K and sorts it into ascending order.
 ! StatsR_K: Subroutine to provide the minmum, maximum, mean, the quadratic mean, and the standard deviation of a 
@@ -3474 +3484 @@
 ! Subroutine swaps the values of its two formal arguments.
 
-      IMPLICIT  NONE
+      IMPLICIT NONE
 
       REAL(r_k), INTENT(INOUT)                           :: a, b
@@ -4044 +4054 @@
   END SUBROUTINE percentiles_R_K4D
 
+  REAL(r_k) FUNCTION distanceRK(pointA, pointB)
+  ! Function to provide the distance between two points
+
+    IMPLICIT NONE
+
+    REAL(r_k), DIMENSION(2), INTENT(in)                  :: pointA, pointB
+
+!!!!!!! Variables
+! pointA, B: couple of points to compute the distance between them
+
+    fname = 'distanceRK'
+
+    distanceRK = SQRT( (pointB(1)-pointA(1))**2 + (pointB(2)-pointA(2))**2 )
+
+  END FUNCTION distanceRK
+
+  REAL(r_k) FUNCTION shoelace_area_polygon(Nvertex, poly)
+  ! Computing the area of a polygon using sholace formula
+  ! FROM: https://en.wikipedia.org/wiki/Shoelace_formula
+
+    IMPLICIT NONE
+
+      INTEGER, INTENT(in)                                :: Nvertex
+      REAL(r_k), DIMENSION(Nvertex,2), INTENT(in)        :: poly
+
+! Local
+      INTEGER                                            :: i
+      REAL(r_k)                                          :: areapos, areaneg
+
+!!!!!!! Variables
+! Nvertex: number of vertices of the polygon
+! poly: coordinates of the vertex of the polygon (sorted)
+
+    fname = 'shoelace_area_polygon'
+
+    areapos = 0.
+    areaneg = 0.
+
+    DO i=1, Nvertex-1
+      areapos = areapos + poly(i,1)*poly(i+1,2)
+      areaneg = areaneg + poly(i+1,1)*poly(i,2)
+    END DO
+
+    areapos = areapos + poly(Nvertex,1)*poly(1,2)
+    areaneg = areaneg + poly(1,1)*poly(Nvertex,2)
+
+    shoelace_area_polygon = 0.5*(areapos - areaneg)
+
+  END FUNCTION shoelace_area_polygon
+
+  SUBROUTINE intersection_2Dlines(lineA, lineB, intersect, ptintersect)
+  ! Subroutine to provide the intersection point between two lines on the plane using Cramer's method
+
+    IMPLICIT NONE
+
+    REAL(r_k), DIMENSION(2,2), INTENT(in)                  :: lineA, lineB
+    LOGICAL, INTENT(out)                                   :: intersect
+    REAL(r_k), DIMENSION(2), INTENT(out)                   :: ptintersect
+
+! Local
+    REAL(r_k), DIMENSION(2)                                :: segmentA, segmentB
+    REAL(r_k)                                              :: a11, a12, a21, a22, z1, z2
+    REAL(r_k)                                              :: det, detX, detY
+    LOGICAL                                                :: axisAx, axisBx, axisAy, axisBy
+
+!!!!!!! Variables
+! lineA: couple of coordinates for the line A
+! lineB: couple of coordinates for the line B
+! intersect: whether two lines intersect
+! ptintersect: point of intersection [(0,0) if they do not intersect]
+
+    fname = 'intersection_2Dlines'
+
+    axisAx = .FALSE.
+    axisAy = .FALSE.
+    axisBx = .FALSE.
+    axisBy = .FALSE.
+    ! Setting segment parameters y = A + B*x
+    IF (lineA(2,1) /= lineA(1,1)) THEN
+      segmentA(2) = (lineA(2,2)-lineA(1,2))/(lineA(2,1)-lineA(1,1))
+      IF ( (lineA(1,1)*segmentA(2) - lineA(1,2)) /= (lineA(2,1)*segmentA(2) - lineA(2,2)) ) THEN
+        PRINT *,'A = y1 - x2*B = ', lineA(1,2) - lineA(1,1)*segmentA(2)
+        PRINT *,'A = y2 - x2*B = ', lineA(2,2) - lineA(2,1)*segmentA(2)
+        msg = 'Wrong calculation of parameter A, for lineA'
+        CALL ErrMSg(msg, fname, -1)
+      END IF
+      segmentA(1) = lineA(1,2) - lineA(1,1)*segmentA(2)
+      a11 = segmentA(2)
+      a12 = -oneRK
+      z1 = -segmentA(1)
+      IF (lineA(2,2) == lineA(1,2)) axisAx = .TRUE.
+    ELSE
+      ! lineA || y-axis
+      axisAy = .TRUE.
+    END IF
+
+    IF (lineB(2,1) /= lineB(1,1)) THEN
+      segmentB(2) = (lineB(2,2)-lineB(1,2))/(lineB(2,1)-lineB(1,1))
+      IF ( (lineB(1,1)*segmentB(2) - lineB(1,2)) /= (lineB(2,1)*segmentB(2) - lineB(2,2)) ) THEN
+        PRINT *,'A = x1*B -y1 = ', lineB(1,1)*segmentB(2) - lineB(1,2)
+        PRINT *,'A = x2*B -y2 = ', lineB(2,1)*segmentB(2) - lineB(2,2)
+        msg = 'Wrong calculation of parameter A, for lineB'
+        CALL ErrMSg(msg, fname, -1)
+      END IF
+      segmentB(1) = lineB(1,2) - lineB(1,1)*segmentB(2)
+      a21 = segmentB(2)
+      a22 = -oneRK
+      z2 = -segmentB(1)
+      IF (lineB(2,2) == lineB(1,2)) axisBx = .TRUE.
+    ELSE
+      ! lineB || y-axis
+      axisBy = .TRUE.
+    END IF
+    ! Cramer's method
+    ! a11 = B1; a12 = -1
+    ! a21 = B2; a22 = -1
+    ! z1 = -A1
+    ! z2 = -A2
+    ! (a11 a12)(x) (z1)
+    ! (a21 a22)(y) (z2)
+    ! -------- ------ ----- ---- --- -- -
+    ! det = (a11*a22-a12*a21)
+    ! detX = (z1*a22-z2*a21)
+    ! detY = (a11*z1-a12*z2)
+    ! ptintercept = (detX/det, detY/det)
+
+    ! Cases when some of the lines are parallel to any given axis
+!    PRINT *,'          axisAx', axisAx, 'axisAy', axisAy, 'axisBx', axisBx, 'axisBy', axisBy
+    IF (axisAx .OR. axisAy .OR. axisBx .OR. axisBy) THEN
+      IF (axisAx) THEN
+        IF (axisBy) THEN
+          intersect = .TRUE.
+          ptintersect(1) = lineB(1,1)
+          ptintersect(2) = lineA(1,2)
+        ELSE
+          intersect = .TRUE.
+          ptintersect(1) = (lineA(1,2)-segmentB(1))/segmentB(2)
+          ptintersect(2) = lineA(1,2)
+        END IF
+      END IF
+      IF (axisAy) THEN
+        IF (axisBy) THEN
+          intersect = .TRUE.
+          ptintersect(1) = lineA(1,1)
+          ptintersect(2) = lineB(1,2)
+        ELSE
+          intersect = .TRUE.
+          ptintersect(1) = lineA(1,1)
+          ptintersect(2) = segmentB(1) + lineA(1,1)*segmentB(2)
+        END IF
+      END IF
+      IF (axisBx) THEN
+        IF (axisAy) THEN
+          intersect = .TRUE.
+          ptintersect(1) = lineA(1,1)
+          ptintersect(2) = lineB(1,2)
+        ELSE
+          intersect = .TRUE.
+          ptintersect(1) = (lineB(1,2)-segmentA(1))/segmentA(2)
+          ptintersect(2) = lineB(1,2)
+        END IF
+      END IF
+      IF (axisBy) THEN
+        IF (axisAx) THEN
+          intersect = .TRUE.
+          ptintersect(1) = lineB(1,1)
+          ptintersect(2) = lineA(1,2)
+        ELSE
+          intersect = .TRUE.
+          ptintersect(1) = lineB(1,1)
+          ptintersect(2) = segmentA(1) + lineB(1,1)*segmentA(2)
+        END IF
+      END IF
+    ELSE
+      det = (a11*a22-a12*a21)
+      ! Parallel lines !
+      IF (det == zeroRK) THEN
+        intersect = .FALSE.
+        ptintersect = zeroRK
+      ELSE
+        intersect = .TRUE.
+        detX = (z1*a22-z2*a12)
+        detY = (a11*z2-a21*z1)
+
+        ptintersect(1) = detX/det
+        ptintersect(2) = detY/det
+      END IF
+    END IF
+
+  END SUBROUTINE intersection_2Dlines
+
+!refs: 
+!https://www.mathopenref.com/heronsformula.html
+!https://math.stackexchange.com/questions/1406340/intersect-area-of-two-polygons-in-cartesian-plan
+!http://www.cap-lore.com/MathPhys/IP/
+!http://www.cap-lore.com/MathPhys/IP/IL.html
+!https://www.element84.com/blog/determining-the-winding-of-a-polygon-given-as-a-set-of-ordered-points
+!https://stackoverflow.com/questions/1165647/how-to-determine-if-a-list-of-polygon-points-are-in-clockwise-order
+!https://en.wikipedia.org/wiki/Shoelace_formula
+!https://en.wikipedia.org/wiki/Winding_number
+!https://en.wikipedia.org/wiki/Simple_polygon
+!https://en.wikipedia.org/wiki/Polygon#Properties
+!https://en.wikipedia.org/wiki/Convex_polygon
+!https://en.wikipedia.org/wiki/Jordan_curve_theorem
+!https://www.sangakoo.com/ca/temes/metode-de-cramer
+!https://www.geogebra.org/m/pw4QHFYT
+
+  SUBROUTINE intersectfaces(faceA, faceB, intersect, intersectpt)
+  ! Subroutine to provide if two faces of two polygons intersect
+  ! AFTER: http://www.cap-lore.com/MathPhys/IP/IL.html
+  !   A: faceA(1,:)
+  !   B: faceA(2,:)
+  !   C: faceB(1,:)
+  !   D: faceB(2,:)
+
+    IMPLICIT NONE
+
+    REAL(r_k), DIMENSION(2,2), INTENT(in)                :: faceA, faceB
+    INTEGER, INTENT(out)                                 :: intersect
+    REAL(r_k), DIMENSION(2), INTENT(out)                 :: intersectpt
+
+! Local
+    REAL(r_k)                                            :: Axmin, Aymin, Axmax, Aymax
+    REAL(r_k)                                            :: Bxmin, Bymin, Bxmax, Bymax
+    REAL(r_k)                                            :: areaABD, areaACD, areaBDC, areaDAB
+    REAL(r_k), DIMENSION(3,2)                            :: triangle
+    LOGICAL                                              :: Lintersect
+
+!!!!!!! Variables
+! faceA/B: coordinates of faces A and B to determine if they intersect
+! intersect: integer to say if they intersect (0, no-intersect, +/-1 intersect)
+! intersectpt: point where faces intersect [(0,0) otherwise]
+
+    fname = 'intersectfaces'
+
+!    PRINT *,'     ' // TRIM(fname) // ' ________'
+!    PRINT *,'            faceA:', faceA(1,:), ';',faceA(2,:)
+!    PRINT *,'            faceB:', faceB(1,:), ';',faceB(2,:)
+
+    Axmin = MINVAL(faceA(:,1))
+    Axmax = MAXVAL(faceA(:,1))
+    Aymin = MINVAL(faceA(:,2))
+    Aymax = MAXVAL(faceA(:,2))
+    Bxmin = MINVAL(faceB(:,1))
+    Bxmax = MAXVAL(faceB(:,1))
+    Bymin = MINVAL(faceB(:,2))
+    Bymax = MAXVAL(faceB(:,2))
+
+    ! No intersection
+    IF ( (Axmax <= Bxmin) .OR. (Axmin >= Bxmax) .OR. (Aymax <= Bymin) .OR. (Aymin >= Bymax) ) THEN
+      intersect = 0
+      intersectpt = zeroRK
+    ELSE
+      ! Triangle ABD
+      triangle(1,:) = faceA(1,:)
+      triangle(2,:) = faceA(2,:)
+      triangle(3,:) = faceB(2,:)
+      areaABD = shoelace_area_polygon(3, triangle)
+  
+      ! Triangle ACD
+      triangle(1,:) = faceA(1,:)
+      triangle(2,:) = faceB(1,:)
+      triangle(3,:) = faceB(2,:)
+      areaACD = shoelace_area_polygon(3, triangle)
+
+      ! Triangle BDC
+      triangle(1,:) = faceA(2,:)
+      triangle(2,:) = faceB(2,:)
+      triangle(3,:) = faceB(1,:)
+      areaBDC = shoelace_area_polygon(3, triangle)
+
+      ! Triangle DAB
+      triangle(1,:) = faceB(2,:)
+      triangle(2,:) = faceA(1,:)
+      triangle(3,:) = faceA(2,:)
+      areaDAB = shoelace_area_polygon(3, triangle)
+
+      IF (areaABD>zeroRK .AND. areaACD>zeroRK .AND. areaBDC>zeroRK .AND. areaDAB>zeroRK) THEN
+        intersect = INT(ABS(areaABD)/areaABD)
+        CALL intersection_2Dlines(faceA, faceB, Lintersect, intersectpt)
+      ELSE IF (areaABD<zeroRK .AND. areaACD<zeroRK .AND. areaBDC<zeroRK .AND. areaDAB<zeroRK) THEN
+        intersect = INT(ABS(areaABD)/areaABD)
+        CALL intersection_2Dlines(faceA, faceB, Lintersect, intersectpt)
+      ELSE
+        intersect = 0
+        intersectpt = zeroRK
+      END IF
+!      PRINT *,'     intersect faces: areaABD',areaABD, 'areaACD', areaACD, 'areaBDC',areaBDC, 'areaDAB',areaDAB, 'prod', &
+!        areaABD*areaACD*areaBDC*areaDAB, 'L:', areaABD*areaACD*areaBDC*areaDAB > zeroRK, 'I', intersect
+
+    END IF
+
+  END SUBROUTINE intersectfaces
+
+  LOGICAL FUNCTION poly_has_point(Nvertex, polygon, point)
+  ! Function to determine if a polygon has already a given point as one of its vertex
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: Nvertex
+    REAL(r_k), DIMENSION(Nvertex,2), INTENT(in)          :: polygon
+    REAL(r_k), DIMENSION(2), INTENT(in)                  :: point
+
+! Local
+    INTEGER                                              :: iv
+    REAL(r_k), DIMENSION(2)                              :: diff
+
+!!!!!!! Vertrex
+! Nvertex: number of vertexs of the polygon
+! polygon: vertexs of the polygon
+! point: point to look for its ownership into the polygon
+
+    fname = 'poly_has_point'
+
+    poly_has_point = .FALSE.
+    DO iv=1, Nvertex
+      diff = polygon(iv,:)-point
+      IF ( (diff(1) == zeroRK) .AND. (diff(2) == zeroRK)) THEN
+        poly_has_point = .TRUE.
+        EXIT
+      END IF
+    END DO
+
+  END FUNCTION poly_has_point
+
+  SUBROUTINE join_polygon(NvertexA, NvertexB, NvertexAB, polyA, polyB, Ncoinvertex, coinpoly)
+  ! Subroutine to join two polygons
+  ! AFTER: http://www.cap-lore.com/MathPhys/IP/ and http://www.cap-lore.com/MathPhys/IP/IL.html
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: NvertexA, NvertexB, NvertexAB
+    REAL(r_k), DIMENSION(NvertexA,2), INTENT(in)         :: polyA
+    REAL(r_k), DIMENSION(NvertexB,2), INTENT(in)         :: polyB
+    INTEGER, INTENT(out)                                 :: Ncoinvertex
+    REAL(r_k), DIMENSION(NvertexAB,2), INTENT(out)        :: coinpoly
+
+! Local
+    INTEGER                                              :: iA, iB, icoin, ii
+    REAL(r_k), DIMENSION(2,2)                            :: face1, face2
+    INTEGER                                              :: intersct
+    REAL(r_k), DIMENSION(2)                              :: ptintersct
+
+
+!!!!!!! variables
+! NvertexA: number of vertexs polygon A
+! NvertexB: number of vertexs polygon B
+! polyA: pairs of coordinates for the polygon A (clockwise)
+! polyB: pairs of coordinates for the polygon B (clockwise)
+! Ncoinvertex: number of vertexes for the coincident polygon
+! coinpoly: pairs of coordinates for the coincident polygon (clockwise)
+
+    fname = 'join_polygon'
+
+    icoin = 0
+    coinpoly = 0.
+
+    ! First, include that vertex which do not lay within any polygon
+    DO iA=1, NvertexA
+      PRINT *, '  iA:', iA, ':', polyA(iA,:), point_inside(polyA(iA,:), NvertexB, polyB)
+      IF (.NOT. point_inside(polyA(iA,:), NvertexB, polyB)) THEN
+        icoin = icoin + 1
+        coinpoly(icoin,:) = polyA(iA,:)
+      END IF
+    END DO
+
+    DO iB=1, NvertexB
+      PRINT *, '  iB:', iB, ':', polyB(iB,:), point_inside(polyB(iB,:), NvertexA, polyA)
+      IF (.NOT. point_inside(polyB(iB,:), NvertexA, polyA)) THEN
+        icoin = icoin + 1
+        coinpoly(icoin,:) = polyB(iB,:)
+      END IF
+    END DO
+
+    DO iA=1, NvertexA
+      ! Getting couple of vertexs from polyA and polyB
+      IF (iA /= NvertexA) THEN
+        face1(1,:) = polyA(iA,:)
+        face1(2,:) = polyA(iA+1,:)
+      ELSE
+        face1(1,:) = polyA(iA,:)
+        face1(2,:) = polyA(1,:)
+      END IF
+      DO iB=1, NvertexB
+        IF (iB /= NvertexB) THEN
+          face2(1,:) = polyB(iB,:)
+          face2(2,:) = polyB(iB+1,:)
+        ELSE
+          face2(1,:) = polyB(iB,:)
+          face2(2,:) = polyB(1,:)
+        END IF
+        
+        ! Compute areas of the four possible triangles. Introduce the coincident vertexs not included
+        CALL intersectfaces(face1, face2, intersct, ptintersct)
+        PRINT *,iA,':',face1(1,:),';',face1(2,:), '=', iB, face2(1,:),';',face2(2,:), '<>', intersct,':', ptintersct
+        IF (intersct == 1) THEN
+          IF (.NOT.poly_has_point(icoin,coinpoly(1:icoin,:),ptintersct) ) THEN
+            icoin = icoin + 1
+            coinpoly(icoin,:) = ptintersct
+          END IF
+        ELSE IF (intersct == -1) THEN
+          IF (.NOT.poly_has_point(icoin,coinpoly(1:icoin,:),ptintersct) ) THEN
+            icoin = icoin + 1
+            coinpoly(icoin,:) = ptintersct
+          END IF
+        END IF
+
+      END DO
+    END DO
+    Ncoinvertex = icoin
+
+  END SUBROUTINE join_polygon
+
+  SUBROUTINE sort_polygon(Nvertex, polygon, sortpoly)
+  ! Subroutine to sort a polygon using its center as average of the coordinates
+  !    Should be used the centroid instead, but by now let do it simple
+  !    https://en.wikipedia.org/wiki/Centroid
+
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: Nvertex
+    REAL(r_k), DIMENSION(Nvertex,2), INTENT(in)          :: polygon
+    REAL(r_k), DIMENSION(Nvertex,2), INTENT(out)         :: sortpoly
+
+! Local
+    INTEGER                                              :: iv, j
+    REAL(r_k)                                            :: vang
+    REAL(r_k), DIMENSION(2)                              :: center
+    REAL(r_k), DIMENSION(Nvertex)                        :: angles
+
+!!!!!!! Variables
+! Nvertex: number of vetexs
+! polygon: coordinates of the vertices of the polygon
+! sortpoly: sorted polygon (clockwise)
+
+    fname = 'sort_polygon'
+
+    ! To be substituted by centroid calculation (which requires already sorted vetexs...)
+    center(1) = SUM(polygon(:,1))/Nvertex
+    center(2) = SUM(polygon(:,2))/Nvertex
+
+    DO iv=1, Nvertex
+      angles(iv) = ATAN2(polygon(iv,2)-center(2),polygon(iv,1)-center(1))
+    END DO
+    CALL sortR_K(angles, Nvertex)
+
+    sortpoly = zeroRK
+    DO iv=1, Nvertex
+      DO j=1, Nvertex
+        vang = ATAN2(polygon(j,2)-center(2), polygon(j,1)-center(1))
+        IF (angles(iv) == vang) THEN
+          sortpoly(iv,:) = polygon(j,:)
+          EXIT
+        END IF
+      END DO
+    END DO
+
+  END SUBROUTINE sort_polygon
+
+  LOGICAL FUNCTION point_inside(point, Nvertex, polygon)
+  ! Function to determine if a given point is inside a polygon providing its sorted vertices
+  ! FROM: https://en.wikipedia.org/wiki/Point_in_polygon
+
+    IMPLICIT NONE
+
+    REAL(r_k), DIMENSION(2), INTENT(in)                  :: point
+    INTEGER, INTENT(in)                                  :: Nvertex
+    REAL(r_k), DIMENSION(Nvertex,2), INTENT(in)          :: polygon
+
+    ! Local
+    INTEGER                                              :: iv, Nintersect
+    INTEGER                                              :: cross
+    REAL(r_k)                                            :: xmin
+    REAL(r_k), DIMENSION(2)                              :: crosspoint
+    REAL(r_k), DIMENSION(2,2)                            :: face1, face2
+    REAL(r_k), DIMENSION(Nvertex)                        :: abovebelow
+
+!!!!!!! Variables
+! point: point to look for
+! Nvertrex: number of vertices of a polygon
+! polygon: vertices of a polygon
+
+    fname = 'point_inside'
+
+    xmin = MINVAL(polygon(:,1))
+
+    ! Looking for the intersection with the ray
+    Nintersect = 0
+    face1(1,:) = (/ xmin-0.5, point(2) /)
+    face1(2,:) = (/ point(1), point(2) /)
+
+    DO iv = 1, Nvertex
+      IF (iv /= Nvertex) THEN
+        face2(1,:) = polygon(iv,:)
+        face2(2,:) = polygon(iv+1,:)
+      ELSE
+        face2(1,:) = polygon(iv,:)
+        face2(2,:) = polygon(1,:)
+      END IF
+      CALL intersectfaces(face1, face2, cross, crosspoint)
+      IF (cross /= 0) THEN
+        Nintersect = Nintersect + 1
+        abovebelow(Nintersect) = iv
+      END IF
+    END DO
+
+    IF (MOD(Nintersect,2) == 0) THEN
+      point_inside = .FALSE.
+    ELSE
+      point_inside = .TRUE.
+    END IF
+
+  END FUNCTION point_inside
+
 END MODULE module_scientific