inpoly.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\inpoly.f
00002 !
00081                      LOGICAL FUNCTION INPOLY
00082 !                    ***********************
00083 !
00084      &( X , Y , XSOM , YSOM , NSOM )
00085 !
00086 !***********************************************************************
00087 ! BIEF   V6P1                                   21/08/2010
00088 !***********************************************************************
00089 !
00090 !
00091 !
00092 !
00093 !
00094 !
00095 !
00096 !
00097 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00098 !| NSOM           |-->| NUMBER OF APICES OF POLYGON
00099 !| X              |-->| ABSCISSA OF POINT
00100 !| Y              |-->| ORDINATE OF POINT
00101 !| XSOM           |-->| ABSCISSAE OF POLYGON APICES
00102 !| YSOM           |-->| ORDINATES OF POLYGON APICES
00103 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00104 !
00105       IMPLICIT NONE
00106       INTEGER LNG,LU
00107       COMMON/INFO/LNG,LU
00108 !
00109 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00110 !
00111       INTEGER, INTENT(IN) :: NSOM
00112       DOUBLE PRECISION, INTENT(IN) :: X,Y
00113       DOUBLE PRECISION, INTENT(IN) :: XSOM(NSOM),YSOM(NSOM)
00114 !
00115 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00116 !
00117       INTEGER N,NSECT
00118 !
00119       DOUBLE PRECISION A,B,ANGLE,XDEP,YDEP,XARR,YARR,DET,MU,LAMBDA,EPS
00120 !
00121       INTRINSIC COS,SIN,ABS,MOD
00122 !
00123 !-----------------------------------------------------------------------
00124 !
00125       EPS = 1.D-9
00126       ANGLE = -1.D0
00127 !
00128 ! CHOOSES A AND B SUCH AS TO AVOID SPECIAL CASES
00129 !
00130 1000  CONTINUE
00131       ANGLE = ANGLE + 1.D0
00132       IF(ANGLE.GT.360.D0) THEN
00133 !       SPECIAL CASE OF A POINT ON THE CONTOUR
00134         INPOLY=.TRUE.
00135         RETURN
00136       ENDIF
00137       A = COS(ANGLE*3.141592653D0/180.D0)
00138       B = SIN(ANGLE*3.141592653D0/180.D0)
00139       NSECT=0
00140 !
00141 ! LOOP ON ALL THE SEGMENTS OF THE POLYGON
00142 !
00143       DO N=1,NSOM
00144 !
00145 !     DEP : 1ST POINT OF THE SEGMENT    ARR : 2ND POINT
00146 !
00147       XDEP=XSOM(N)
00148       YDEP=YSOM(N)
00149       IF(N.LT.NSOM) THEN
00150         XARR=XSOM(N+1)
00151         YARR=YSOM(N+1)
00152       ELSE
00153         XARR=XSOM(1)
00154         YARR=YSOM(1)
00155       ENDIF
00156 !
00157 !     CASE WHERE TWO SUCCESSIVE POINTS ARE DUPLICATES
00158 !
00159       IF(ABS(XDEP-XARR)+ABS(YDEP-YARR).LT.EPS) THEN
00160         WRITE(LU,*) ' '
00161         WRITE(LU,*) ' '
00162         IF(LNG.EQ.1) THEN
00163           WRITE(LU,*) 'INPOLY : POINTS CONFONDUS DANS LE POLYGONE'
00164           WRITE(LU,*) 'AU POINT DE COORDONNEES : ',XDEP,'  ET  ',YDEP
00165           WRITE(LU,*) 'DE NUMERO ',N
00166           IF(N.EQ.NSOM) THEN
00167             WRITE(LU,*) 'LE DERNIER POINT NE DOIT PAS ETRE EGAL AU'
00168             WRITE(LU,*) 'PREMIER (POUR UN TRIANGLE, PAR EXEMPLE,'
00169             WRITE(LU,*) 'IL FAUT DONNER TROIS POINTS ET NON QUATRE)'
00170           ENDIF
00171           WRITE(LU,*) 'INPOLY EST PROBABLEMENT CALLED BY FILPOL'
00172         ENDIF
00173         IF(LNG.EQ.2) THEN
00174           WRITE(LU,*) 'INPOLY: SUPERIMPOSED POINTS IN THE POLYGON'
00175           WRITE(LU,*) 'AT POINT: ',XDEP,'  AND  ',YDEP,' WITH NUMBER ',N
00176           IF(N.EQ.NSOM) THEN
00177             WRITE(LU,*) 'THE LAST POINT MUST NOT BE EQUAL TO THE FIRST'
00178             WRITE(LU,*) 'FOR EXAMPLE, GIVE 3 POINTS FOR A TRIANGLE'
00179           ENDIF
00180           WRITE(LU,*) 'INPOLY IS PROBABLY CALLED BY FILPOL'
00181         ENDIF
00182         WRITE(LU,*) ' '
00183         WRITE(LU,*) ' '
00184         CALL PLANTE(1)
00185         STOP
00186       ENDIF
00187 !
00188 !     CASE WHERE THE POINT IS A VERTEX
00189 !     (THE GENERAL ALGORITHM WOULD DEFINE IT AS EXTERNAL WITH 2 INTERSECTIONS)
00190 !
00191       IF(ABS(X-XDEP).LE.EPS.AND.ABS(Y-YDEP).LE.EPS) THEN
00192         NSECT=1
00193         EXIT
00194       ENDIF
00195 !
00196 !     DETERMINANT OF THE KRAMER SYSTEM
00197 !
00198       DET = A*(YDEP-YARR)-B*(XDEP-XARR)
00199       IF(ABS(DET).LT.EPS) GO TO 1000
00200 !
00201       MU     = ( (XDEP-X)*(YDEP-YARR)-(YDEP-Y)*(XDEP-XARR) ) / DET
00202       LAMBDA = (    A    *(YDEP-Y   )-    B   *(XDEP-X   ) ) / DET
00203 !
00204 !-------------------------------------------------------
00205 ! JP RENAUD (CSN BRISTOL) CORRECTION TO AVOID THAT THE INTERSECTION
00206 ! POINT BE ONE OF THE VRTICES
00207 !
00208 ! IF THE INTERSECTION POINT IS A VERTEX, INCREASES THE ANGLE
00209 ! OTHERWISE THE POINT WOULD BE COUNTED TWICE INSTEAD OF JUST ONCE
00210 !
00211       IF ((ABS(X+A*MU-XDEP).LE.EPS.AND.ABS(Y+B*MU-YDEP).LE.EPS)
00212      &    .OR.
00213      &    (ABS(X+A*MU-XARR).LE.EPS.AND.ABS(Y+B*MU-YARR).LE.EPS))
00214      &    GOTO 1000
00215 !
00216 ! END OF JP RENAUD CORRECTION
00217 !-------------------------------------------------------
00218 !
00219       IF(MU.GE.-EPS.AND.LAMBDA.GE.-EPS.AND.LAMBDA.LE.1.D0+EPS) THEN
00220         NSECT=NSECT+1
00221       ENDIF
00222 !
00223       ENDDO ! N
00224 !
00225       INPOLY=(MOD(NSECT,2).EQ.1)
00226 !
00227 !-----------------------------------------------------------------------
00228 !
00229       RETURN
00230       END

Generated on Fri Aug 31 2013 18:12:58 by S.E.Bourban (HRW) using doxygen 1.7.0