inpoly.f

Go to the documentation of this file.

!                   ***********************
                    LOGICAL FUNCTION inpoly
!                   ***********************
!
     &( x , y , xsom , ysom , nsom )
!
!***********************************************************************
! BIEF   V6P1                                   21/08/2010
!***********************************************************************
!
!brief    INDICATES IF A POINT WITH COORDINATES X AND Y IS
!+                IN A POLYGON WITH GIVEN VERTICES.
!code
!+    PRINCIPLE: TAKES HALF A LINE STARTING FROM THE POINT AND COUNTS THE
!+               NUMBER OF TIMES IT INTERSECTS WITH THE POLYGON
!+
!+    ALSO WORKS IF THE POLYGON IS NOT CONVEX
!+
!+    INTERSECTIONS ARE IDENTIFIED USING THE LINES PARAMETRIC EQUATIONS :
!+
!+
!+    X + A * MU = XDEP + (XARR-XDEP) * LAMBDA
!+    Y + B * MU = YDEP + (YARR-YDEP) * LAMBDA
!+
!+    THE HALF-LINE IS CHARACTERISED BY THE CHOICE OF A AND B, AND THE
!+    SIGN OF MU. THERE IS INTERSECTION IF MU > 0 AND 0 < LAMBDA < 1
!
!warning  THE POLYGON VERTICES MUST BE DISTINCT (NO DUPLICATE NODES)
!
!history  E. DAVID (LHF)
!+
!+
!+   ORIGINAL IDEA AND CODE
!
!history  J.-M. HERVOUET (LNH)
!+        18/06/96
!+        V5P2
!+
!
!history  JEAN-PHILIPPE RENAUD (CSN BRISTOL)
!+        27/07/99
!+
!+   CORRECTION FOR A SPECIAL CASE
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        13/07/2010
!+        V6P0
!+   Translation of French comments within the FORTRAN sources into
!+   English comments
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        21/08/2010
!+        V6P0
!+   Creation of DOXYGEN tags for automated documentation and
!+   cross-referencing of the FORTRAN sources
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| NSOM           |-->| NUMBER OF APICES OF POLYGON
!| X              |-->| ABSCISSA OF POINT
!| Y              |-->| ORDINATE OF POINT
!| XSOM           |-->| ABSCISSAE OF POLYGON APICES
!| YSOM           |-->| ORDINATES OF POLYGON APICES
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: NSOM
      DOUBLE PRECISION, INTENT(IN) :: X,Y
      DOUBLE PRECISION, INTENT(IN) :: XSOM(nsom),YSOM(nsom)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER N,NSECT
!
      DOUBLE PRECISION A,B,ANGLE,XDEP,YDEP,XARR,YARR,DET,MU,LAMBDA,EPS
!
      DOUBLE PRECISION PI
      INTRINSIC cos,sin,abs,mod
!
!-----------------------------------------------------------------------
!
      eps = 1.d-9
      angle = -1.d0
!
! CHOOSES A AND B SUCH AS TO AVOID SPECIAL CASES
!
1000  CONTINUE
      angle = angle + 1.d0
      IF(angle.GT.360.d0) THEN
!       SPECIAL CASE OF A POINT ON THE CONTOUR
        inpoly=.true.
        RETURN
      ENDIF
      pi = 4.d0 * atan( 1.d0 )
      a = cos( angle*pi/180.d0 )
      b = sin( angle*pi/180.d0 )
      nsect=0
!
! LOOP ON ALL THE SEGMENTS OF THE POLYGON
!
      DO n=1,nsom
!
!     DEP : 1ST POINT OF THE SEGMENT    ARR : 2ND POINT
!
      xdep=xsom(n)
      ydep=ysom(n)
      IF(n.LT.nsom) THEN
        xarr=xsom(n+1)
        yarr=ysom(n+1)
      ELSE
        xarr=xsom(1)
        yarr=ysom(1)
      ENDIF
!
!     CASE WHERE TWO SUCCESSIVE POINTS ARE DUPLICATES
!
      IF(abs(xdep-xarr)+abs(ydep-yarr).LT.eps) THEN
        WRITE(lu,*) ' '
        WRITE(lu,*) ' '
        WRITE(lu,*) 'INPOLY: SUPERIMPOSED POINTS IN THE POLYGON'
        WRITE(lu,*) 'AT POINT: ',xdep,'  AND  ',ydep,' WITH NUMBER ',n
        IF(n.EQ.nsom) THEN
          WRITE(lu,*) 'THE LAST POINT MUST NOT BE EQUAL TO THE FIRST'
          WRITE(lu,*) 'FOR EXAMPLE, GIVE 3 POINTS FOR A TRIANGLE'
        ENDIF
        WRITE(lu,*) 'INPOLY IS PROBABLY CALLED BY FILPOL'
        WRITE(lu,*) ' '
        WRITE(lu,*) ' '
        CALL plante(1)
        stop
      ENDIF
!
!     CASE WHERE THE POINT IS A VERTEX
!     (THE GENERAL ALGORITHM WOULD DEFINE IT AS EXTERNAL WITH 2 INTERSECTIONS)
!
      IF(abs(x-xdep).LE.eps.AND.abs(y-ydep).LE.eps) THEN
        nsect=1
        EXIT
      ENDIF
!
!     DETERMINANT OF THE KRAMER SYSTEM
!
      det = a*(ydep-yarr)-b*(xdep-xarr)
      IF(abs(det).LT.eps) GO TO 1000
!
      mu     = ( (xdep-x)*(ydep-yarr)-(ydep-y)*(xdep-xarr) ) / det
      lambda = (    a    *(ydep-y   )-    b   *(xdep-x   ) ) / det
!
!-------------------------------------------------------
! JP RENAUD (CSN BRISTOL) CORRECTION TO AVOID THAT THE INTERSECTION
! POINT BE ONE OF THE VRTICES
!
! IF THE INTERSECTION POINT IS A VERTEX, INCREASES THE ANGLE
! OTHERWISE THE POINT WOULD BE COUNTED TWICE INSTEAD OF JUST ONCE
!
      IF ((abs(x+a*mu-xdep).LE.eps.AND.abs(y+b*mu-ydep).LE.eps)
     &    .OR.
     &    (abs(x+a*mu-xarr).LE.eps.AND.abs(y+b*mu-yarr).LE.eps))
     &    GOTO 1000
!
! END OF JP RENAUD CORRECTION
!-------------------------------------------------------
!
      IF(mu.GE.-eps.AND.lambda.GE.-eps.AND.lambda.LE.1.d0+eps) THEN
        nsect=nsect+1
      ENDIF
!
      ENDDO ! N
!
      inpoly=(mod(nsect,2).EQ.1)
!
!-----------------------------------------------------------------------
!
      RETURN
      END