geom.f

Go to the documentation of this file.

C> @file geom.f
C> Basic geometrical calculations
C GRIDMAN, grid managment library. Author: Vladislav Kotov, v.kotov@fz-juelich.de

!    Copyright (c) 2017  Forschungszentrum Juelich GmbH
!                        Vladislav Kotov 
!
!    This file is part of GRIDMAN.
!
!    GRIDMAN is free software: you can redistribute it and/or modify
!    it under the terms of the GNU General Public License as published by
!    the Free Software Foundation, either version 3 of the License, or
!    (at your option) any later version.
!
!    GRIDMAN is distributed in the hope that it will be useful,
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU General Public License for more details.
!
!    You should have received a copy of the GNU General Public License
!    along with GRIDMAN.  If not, see <http://www.gnu.org/licenses/>.

C> Area of a triangle
       FUNCTION gridman_triaarea(X1,Y1,X2,Y2,X3,Y3)

       USE gridman,ONLY:gridman_dp
       IMPLICIT NONE
       INTRINSIC abs

       REAL(GRIDMAN_DP),INTENT(IN) :: X1,Y1,X2,Y2,X3,Y3
       REAL(GRIDMAN_DP) :: GRIDMAN_TRIAAREA,S

       s=x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2)
       gridman_triaarea=0.5*abs(s)

       END FUNCTION gridman_triaarea


C> Volume of the solid of revolution obtained by
!! rotating of a triangle around axis X=0
      FUNCTION gridman_triavol(X1,Y1,X2,Y2,X3,Y3)

       USE gridman,ONLY:gridman_dp,gridman_2pi
       IMPLICIT NONE
 
       REAL(GRIDMAN_DP),INTENT(IN) :: X1,Y1,X2,Y2,X3,Y3
       REAL(GRIDMAN_DP) :: GRIDMAN_TRIAVOL,GRIDMAN_TRIAAREA
       REAL(GRIDMAN_DP),PARAMETER :: PI3=gridman_2pi/3.0

       gridman_triavol=pi3*(x1+x2+x3)*
     *                 gridman_triaarea(x1,y1,x2,y2,x3,y3)

       END FUNCTION gridman_triavol


C> Distance between two points on a plane
       FUNCTION gridman_dist2d(X1,Y1,X2,Y2)

       USE gridman,ONLY:gridman_dp
       IMPLICIT NONE
       INTRINSIC sqrt

       REAL(GRIDMAN_DP),INTENT(IN) :: X1,Y1,X2,Y2
       REAL(GRIDMAN_DP) :: GRIDMAN_DIST2D

       gridman_dist2d=sqrt((x1-x2)**2+(y1-y2)**2)

       END FUNCTION gridman_dist2d


C> Vector (XY,YN) normal to (X0,Y0)
       SUBROUTINE gridman_norm2d(X0,Y0,XN,YN)

       USE gridman,ONLY:gridman_dp
       IMPLICIT NONE

       REAL(GRIDMAN_DP),INTENT(IN) :: X0,Y0
       REAL(GRIDMAN_DP),INTENT(OUT) :: XN,YN

         xn=-y0
         yn=x0

       END SUBROUTINE gridman_norm2d


C> Area of the conical surface obtained by rotation 
!! of the interval (X1,Y1),(X2,Y2) around X=0
       FUNCTION gridman_conarea(X1,Y1,X2,Y2)

       USE gridman,ONLY:gridman_dp,gridman_pi
       IMPLICIT NONE
       INTRINSIC sqrt

       REAL(GRIDMAN_DP),INTENT(IN) :: X1,Y1,X2,Y2
       REAL(GRIDMAN_DP) :: GRIDMAN_CONAREA

        gridman_conarea=gridman_pi*(x1+x2)*sqrt((x1-x2)**2+(y1-y2)**2)

       END FUNCTION gridman_conarea

C> Intersection of "particle trajectory" with an interval (2D)
C>
C> Find a point where a "particle" with velocity (VX0,VY0) which starts
!! from point (X0,Y0) intersects the interval [(X1,Y1)(X2,Y2)]. 
!! The function returns .TRUE. if intersection exist 
!! and .FALSE. otherwise
       FUNCTION gridman_cross2d(X0,Y0,V,U,X1,Y1,X2,Y2,XI,YI)
       USE gridman,ONLY:gridman_dp
       IMPLICIT NONE
       INTRINSIC abs,tiny

C> X-coordinate of the "particle" starting point 
       REAL(GRIDMAN_DP),INTENT(IN) :: X0
C> Y-coordinate of the "particle" starting point 
       REAL(GRIDMAN_DP),INTENT(IN) :: Y0      
C> X-component of the particle velocity
       REAL(GRIDMAN_DP),INTENT(IN) :: V
C> Y-component of the particle velocity
       REAL(GRIDMAN_DP),INTENT(IN) :: U    
C> X-coordinate of the first end of the interval
       REAL(GRIDMAN_DP),INTENT(IN) :: X1
C> Y-coordinate of the first end of the interval
       REAL(GRIDMAN_DP),INTENT(IN) :: Y1
C> X-coordinate of the second end of the interval
       REAL(GRIDMAN_DP),INTENT(IN) :: X2
C> Y-coordinate of the second end of the interval
       REAL(GRIDMAN_DP),INTENT(IN) :: Y2 
C> X-coordinate of the intersection point
       REAL(GRIDMAN_DP),INTENT(OUT) :: XI
C> Y-coordinate of the intersection point
       REAL(GRIDMAN_DP),INTENT(OUT) :: YI
       LOGICAL GRIDMAN_CROSS2D
      
       REAL(GRIDMAN_DP) :: D,T

        d=v*(y2-y1)+u*(x1-x2) 
        IF(abs(d).LT.10.*tiny(d)) THEN
         gridman_cross2d=.false.
         RETURN
        END IF 
        t=x1*y2-x2*y1+(y1-y2)*x0+(x2-x1)*y0
        t=t/d
        IF(.NOT.t.GT.0.) THEN !WRONG DIRECTION
         gridman_cross2d=.false.
         RETURN
        END IF
        xi=x0+v*t
        yi=y0+u*t

        IF( ( xi.GT.x1.AND.xi.GT.x2 ) .OR.   !INTERSECTION POINT IS
     f      ( xi.LT.x1.AND.xi.LT.x2 ) .OR.
     f      ( yi.GT.y1.AND.yi.GT.y2 ) .OR.   !OUTSIDE OF THE INTERVAL
     f      ( yi.LT.y1.AND.yi.LT.y2 ) ) THEN 
          gridman_cross2d=.false.
        ELSE
          gridman_cross2d=.true.
        END IF

       END FUNCTION gridman_cross2d

C> Intersection of two intervals on a plane
C>
C> Find intersection point of two intervals on a plane.
!! The function returns .TRUE. if intersection exist and .FALSE. otherwise
       FUNCTION gridman_intersect2d(X11,Y11,X12,Y12,
     f                              x21,y21,x22,y22,xi,yi)
       USE gridman,ONLY:gridman_dp
       IMPLICIT NONE
       INTRINSIC tiny,abs

C> X-coordinate of the first end of the first interval
       REAL(GRIDMAN_DP),INTENT(IN) :: X11
C> Y-coordinate of the first end of the first interval
       REAL(GRIDMAN_DP),INTENT(IN) :: Y11
C> X-coordinate of the second end of the first interval
       REAL(GRIDMAN_DP),INTENT(IN) :: X12
C> Y-coordinate of the second end of the first interval
       REAL(GRIDMAN_DP),INTENT(IN) :: Y12
C> X-coordinate of the first end of the second interval
       REAL(GRIDMAN_DP),INTENT(IN) :: X21
C> Y-coordinate of the first end of the second interval
       REAL(GRIDMAN_DP),INTENT(IN) :: Y21
C> X-coordinate of the second end of the second interval
       REAL(GRIDMAN_DP),INTENT(IN) :: X22
C> Y-coordinate of the second end of the second interval
       REAL(GRIDMAN_DP),INTENT(IN) :: Y22 
C> X-coordinate of the intersection point
       REAL(GRIDMAN_DP),INTENT(OUT) :: XI
C> Y-coordinate of the intersection point
       REAL(GRIDMAN_DP),INTENT(OUT) :: YI                
       LOGICAL :: GRIDMAN_INTERSECT2D
       
       REAL(GRIDMAN_DP) :: A1,B1,C1,A2,B2,C2,DET

        a1=y12-y11
        b1=x11-x12
        c1=y11*(x12-x11)-x11*(y12-y11)
        a2=y22-y21
        b2=x21-x22
        c2=y21*(x22-x21)-x21*(y22-y21)
        det=a1*b2-a2*b1
        IF(abs(det).GT.tiny(det)) THEN
         det=1/det
         xi=(b1*c2-b2*c1)*det
         yi=(c1*a2-c2*a1)*det
C        CHECK IF THE POINT IS OUT OF THE 1st INTERVAL
         IF(abs(a1).GT.abs(b1)) THEN
C         Y-PROJECTION OF THE 1st INTERVAL IS LARGER 
          IF( ( yi.GT.y11.AND.yi.GT.y12 ) .OR.
     f        ( yi.LT.y11.AND.yi.LT.y12 ) ) THEN
           gridman_intersect2d=.false.
           RETURN
          END IF    
         ELSE
C         X-PROJECTION OF THE 1st INTERVAL IS LARGER
          IF( ( xi.GT.x11.AND.xi.GT.x12 ) .OR.
     f        ( xi.LT.x11.AND.xi.LT.x12 ) ) THEN
           gridman_intersect2d=.false.
           RETURN
          END IF 
         END IF   
C        CHECK IF THE POINT IS OUT OF THE 2nd INTERVAL
         IF(abs(a2).GT.abs(b2)) THEN
C         Y-PROJECTION OF THE 2ndT INTERVAL IS LARGER
          IF( ( yi.GT.y21.AND.yi.GT.y22 ) .OR.
     f        ( yi.LT.y21.AND.yi.LT.y22 )) THEN
           gridman_intersect2d=.false.
           RETURN
          END IF    
         ELSE
C         X-PROJECTION OF THE 2nd INTERVAL IS LARGER
          IF( ( xi.GT.x21.AND.xi.GT.x22 ) .OR.
     f        ( xi.LT.x21.AND.xi.LT.x22 )) THEN
           gridman_intersect2d=.false.
           RETURN
          END IF
         END IF    
         gridman_intersect2d=.true.
        ELSE 
         gridman_intersect2d=.false. !PARALLEL LINES
        END IF

       END FUNCTION gridman_intersect2d

C> Find if point lies inside a closed polygon
C>
C> The function returns .TRUE. if the point (X0,Y0) lies inside the 
C> polygon (XP,YP) and .FALSE. otherwise
C>
C> It is assumed that the polygon is simple and closed.
C> That is, no intesection between polygon edges, first and
C> last vertex of the polygon are the same. 
C>
C> WARNING: this function does not check that the polygon is 
C>          simple and closed 
C>
C> Crossing number method -  even-odd test - is applied 
      FUNCTION gridman_point_in_polygon(XP,YP,NP,X0,Y0) 
       USE gridman,ONLY:gridman_dp,gridman_sp,gridman_tol
       IMPLICIT NONE
       INTRINSIC tiny,abs,mod

C>     Number of the polygon vertices - first and last vertex is counted twice  
       INTEGER(GRIDMAN_SP),INTENT(IN) :: NP
C>     X-coordinates of the vertices
       REAL(GRIDMAN_DP) :: XP(np)
C>     Y-coordinates of the vertices
       REAL(GRIDMAN_DP) :: YP(np)
C>     X-coordinate of the point to be checked 
       REAL(GRIDMAN_DP) :: X0
C>     Y-coordinate of the point to be checked 
       REAL(GRIDMAN_DP) :: Y0
       
       LOGICAL :: GRIDMAN_POINT_IN_POLYGON

       REAL(GRIDMAN_DP) :: Y1,Y2,DY,X1,X2,DX,XI
       INTEGER(GRIDMAN_SP) :: IP,CN

C       COUNT INTERSECTIONS OF POLYGON WITH THE HORIZONTAL RAY 
C       DRAWN FROM THE POINT (X0,Y0) TO THE RIGHT
        cn=0
        DO ip=1,np-1
         y1=yp(ip)
         y2=yp(ip+1)
         dy=y2-y1
         IF(.NOT.abs(dy).GT.10*tiny(dy)) cycle !JUST IN CASE ...
         IF(dy.GT.0.) THEN
C         POLYGON EDGE GOES UPWARDS       
          y1=y1-gridman_tol*abs(dy)
          IF(y1.LT.y0.AND.y0.LT.y2) THEN
           x1=xp(ip)
           x2=xp(ip+1)
           dx=x2-x1
           xi=x1+dx/dy*(y0-y1)
C          COUNT ONLY INTERSECTIONS TO THE RIGHT FROM THE CHECKED POINT
           IF(xi.GT.x0) cn=cn+1   
          END IF
         ELSEIF(dy.LT.0.) THEN
C         POLYGON EDGE GOES DOWNWARDS
          y2=y2-gridman_tol*abs(dy)
          IF(y2.LT.y0.AND.y0.LT.y1) THEN
           x1=xp(ip)
           x2=xp(ip+1)
           dx=x2-x1
           xi=x1+dx/dy*(y0-y1)
C          COUNT ONLY INTERSECTIONS TO THE RIGHT FROM THE CHECKED POINT
           IF(xi.GT.x0) cn=cn+1 
          END IF
         END IF  
        END DO  

        IF(mod(cn,2).EQ.0) THEN
         gridman_point_in_polygon=.false.
        ELSE 
         gridman_point_in_polygon=.true.
        END IF
 
      END FUNCTION gridman_point_in_polygon

C> Find if point is close to the interval within prescribed tolerance
C>
C> Function returns .TRUE. if the point (X,Y) lies on the interval
C> [(X1,Y1),(X2,Y2)]) withing relative tolerance TOL and .FALSE. otherwise
      FUNCTION gridman_close2interval2d(X1,Y1,X2,Y2,X,Y,TOL)
       USE gridman,ONLY:gridman_dp
       IMPLICIT NONE
       INTRINSIC tiny,abs

C> X-coordinate of the 1st end of the interval
       REAL(GRIDMAN_DP),INTENT(IN) :: X1
C> Y-coordinate of the 1st end of the interval
       REAL(GRIDMAN_DP),INTENT(IN) :: Y1
C> X-coordinate of the 2nd end of the interval
       REAL(GRIDMAN_DP),INTENT(IN) :: X2 
C> Y-coordinate of the 2nd end of the interval
       REAL(GRIDMAN_DP),INTENT(IN) :: Y2
C> X-coordinate of the point to be checked
       REAL(GRIDMAN_DP),INTENT(IN) :: X 
C> Y-coordinate of the point to be checked
       REAL(GRIDMAN_DP),INTENT(IN) :: Y
C> Tolerance parameter which defines relative accuracy, e.g. 1e-7
       REAL(GRIDMAN_DP),INTENT(IN) :: TOL
       LOGICAL :: GRIDMAN_CLOSE2INTERVAL2D

       REAL(GRIDMAN_DP) :: D,DX,DY,EPS

       dx=y2-y1
       dy=x1-x2
       d=abs(dx*x+dy*y+x2*y1-x1*y2)       
       eps=abs(dx*x)*tol+abs(dy*y)*tol+10.*tiny(x)
       IF(d.LT.eps) THEN
C       POINT (X,Y) IS CLOSE TO THE LINE DRAWN THROUGH (X1,Y1) AND(X2,Y2)
        IF( abs(dy).GT.abs(dx) ) THEN
C        X-PROJECTION OF [(X1,Y1),(X2,Y2)] IS LARGER THAN THE Y-PROJECTION
C          CHECK THAN X IS IN THE INTERVAL [X1,X2]
         eps=abs(x)*tol+10.*tiny(x)
         IF( (x.GT.x1-eps.AND.x.LT.x2+eps).OR.  
     f       (x.GT.x2-eps.AND.x.LT.x1+eps)) THEN
C         X IN THE INTERVAL [X1,X2]
          gridman_close2interval2d=.true.
         ELSE
          gridman_close2interval2d=.false.
         END IF
        ELSE
C        Y-PROJECTION OF [(X1,Y1),(X2,Y2)] IS LARGER THAN THE X-PROJECTION
C          CHECK THAN Y IS IN THE INTERVAL [Y1,Y2]
         eps=abs(y)*tol+10.*tiny(y)
         IF( (y.GT.y1-eps.AND.y.LT.y2+eps).OR.  
     f       (y.GT.y2-eps.AND.y.LT.y1+eps)) THEN
          gridman_close2interval2d=.true.
         ELSE
          gridman_close2interval2d=.false.
         END IF
        END IF
       ELSE
        gridman_close2interval2d=.false.
       END IF       

      END  FUNCTION gridman_close2interval2d