table2d_intersect.F File Reference

#include "implicit_f.inc"

Go to the source code of this file.

Functions/Subroutines
subroutine	table2d_intersect (table, i1, i2, npt, xfac, xint, yint)

Function/Subroutine Documentation

table2d_intersect()

subroutine table2d_intersect	(	type(ttable), intent(in)	table,
		integer, intent(in)	i1,
		integer, intent(in)	i2,
		integer, intent(in)	npt,
		intent(in)	xfac,
		intent(out)	xint,
		intent(out)	yint )

Definition at line 30 of file table2d_intersect.F.

C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE table_mod
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
      INTEGER      ,INTENT(IN)  :: I1,I2,NPT
      my_real      ,INTENT(IN)  :: xfac
      TYPE(TTABLE) ,INTENT(IN)  :: TABLE
      my_real      ,INTENT(OUT) :: xint,yint
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: I,J1,J2,K
      my_real :: s1,s2,t1,t2,x1,x2,y1,y2,ax,ay,bx,by,cx,cy,dm,alpha,beta
c-----------------------------------------------
c     This routine checks if the functions in a 2 dim table do not intersect 
c     with respect to the second independent variable
C=======================================================================
c     Check segment intersections between 2 functions
 
      xint  = zero
      yint  = zero
      
c     I1 = index of first  strain rate 
c     I2 = index of second strain rate 
c
      j1 = (i1 - 1)*npt
      j2 = (i2 - 1)*npt
 
      DO k = 2,npt
        s1 = table%X(1)%VALUES(k-1)*xfac 
        s2 = table%X(1)%VALUES(k)  *xfac
        x1 = s1
        x2 = s2
        t1 = table%Y%VALUES(j1 + k-1) 
        t2 = table%Y%VALUES(j1 + k) 
        y1 = table%Y%VALUES(j2 + k-1) 
        y2 = table%Y%VALUES(j2 + k) 
c
        ax = x2 - x1
        ay = y2 - y1
        bx = s1 - s2
        by = t1 - t2
        dm = ay*bx - ax*by
        IF (dm /= zero) THEN  ! check if segments are not parallel
          cx = s1 - x1
          cy = t1 - y1
          alpha = (bx * cy - by * cx) / dm
          beta  = (ax * cy - ay * cx) / dm
          IF (alpha >= zero .and. alpha < one .and.
     .        beta  <= zero .and. beta  >-one .and. s1 > zero) THEN
            xint = x1 + alpha * ax
            yint = y1 + alpha * ay
            EXIT
          ENDIF
        ENDIF
      END DO
c-----------
      RETURN