C***************************************************************
C       Subroutine VG2LC by boundary-element-method.com                  
C***************************************************************
C
C VG2LC.FOR Version 1, Sept 2015
C by Stephen Kirkup
C http://www.researchgate.net/profile/Stephen_Kirkup
C  School of Engineering
C  University of Central Lancashire

C A test on the geometrical input to *2LC routines in the 2D
C boundary element method.
C Source of the code:
C  http://www.boundary-element-methods.com/fortran/VG2LC.FOR
C Source of the user-guide:
C  http://www.boundary-element-methods.com/fortran/VG2LC_FOR.pdf
C Source of the test code:
C  http://www.boundary-element-methods.com/fortran/VG2LC_TESTS.FOR

      LOGICAL FUNCTION VG2LC(P,VECP,QA,QB,LPONEL,EGEOM,IU)
      REAL*8 P(2),VECP(2),QA(2),QB(2)
      LOGICAL LPONEL
      INTEGER IU
      REAL*8 EGEOM
      REAL*8 QBMA(2),PMQA(2),PMQB(2),NORMQ(2)
      REAL*8 QLEN,PQALEN,PQBLEN
      REAL*8 DIFF(2),QSTAR(2)
      REAL*8 XSTAR

      REAL*8 SIZE2,SSIZE2,DIST2
      CALL VEC2(QA,QB,QBMA)
      QLEN=SIZE2(QBMA)
      CALL VEC2(P,QA,PMQA)
      CALL VEC2(P,QB,PMQB)
      PQALEN=SIZE2(PMQA)
      PQBLEN=SIZE2(PMQB)

      VG2LC=.TRUE.

C Verify EGEOM
C  Verify that EGEOM is positive
      IF (EGEOM.LE.0.0D0) THEN
        VG2LC=.FALSE.
        WRITE(IU,*) 'ERROR(*2LC) - EGEOM must be positive'
      END IF 
C  Verify that EGEOM is not too large
      IF (EGEOM.GT.(QLEN+PQALEN+PQBLEN)*0.001) THEN
        VG2LC=.FALSE.
        WRITE(IU,*) 'ERROR(*2LC) - EGEOM is set too large'
      END IF



C Verify that QA and QB are not the same point
      IF (QLEN<EGEOM*(SIZE2(QB)+SIZE2(QA)+2)) THEN
        VG2LC=.FALSE.
        WRITE(IU,*) 'ERROR(*2LC) - The points QA and QB must be'
        WRITE(IU,*) ' distinct'
      END IF     

C Verify VECP
C  Verify that VECP has unit magnitude
      IF (ABS(SSIZE2(VECP)-1.0D0).GT.EGEOM) THEN
        VG2LC=.FALSE.
        WRITE(IU,*) 'ERROR(*2LC) - VECP must be a unit vector'
      END IF 
C  Verify that VECP is equal to the unit normal to the panel when     
C    the point P lies on the panel (LPONEL=.TRUE.)
      IF (LPONEL) THEN
        CALL NORM2(QA,QB,NORMQ)
        CALL VEC2(NORMQ,VECP,DIFF)
        IF (ABS(SIZE2(DIFF)).GT.2*EGEOM) THEN
          VG2LC=.FALSE.
          WRITE(IU,*) 'ERROR(*2LC) - If LPONEL=.TRUE then VECP be the' 
          WRITE(IU,*) ' unit leftward normal to panel QA-QB'
        END IF
       END IF  
             
C Verify LPONEL and that P does not coincide with QA or QB
C  Set QSTAR; the point on the panel QA-QB that is nearest to P
      
      IF (VG2LC) THEN
        XSTAR=-(QBMA(1)*PMQA(1)+QBMA(2)*PMQA(2))/(QLEN**2)
        IF (XSTAR.LE.0.0D0) THEN
          QSTAR(1)=QA(1)
          QSTAR(2)=QA(2)
        ELSE IF (XSTAR.GE.1.0D0) THEN
          QSTAR(1)=QB(1)
          QSTAR(2)=QB(2)
        ELSE
          QSTAR(1)=QA(1)+XSTAR*QBMA(1)  
          QSTAR(2)=QA(2)+XSTAR*QBMA(2)  
        END IF

C  If LPONEL=.TRUE. then QSTAR must be the same point as P. 
        IF (LPONEL) THEN
          IF (DIST2(P,QSTAR).GT.EGEOM*(SIZE2(QSTAR)+SIZE2(P)+2)) THEN
            VG2LC=.FALSE.
            WRITE(IU,*) 'ERROR(*2LC) - If LPONEL=.TRUE. then the '
            WRITE(IU,*) 'point P must lie on QA-QB'
          END IF
C  If LPONEL=.FALSE. then QSTAR must not be the same point as P
        ELSE
          IF (DIST2(P,QSTAR).LT.EGEOM*(SIZE2(QSTAR)+SIZE2(P)+2)) THEN
            VG2LC=.FALSE.
            WRITE(IU,*) 'ERROR(*2LC) - If LPONEL=FALSE. then the '
            WRITE(IU,*) 'point P must not lie on QA-QB'
          END IF
        END IF
      END IF


C  Verify that P is not coincident with QA
      IF (DIST2(P,QA).LT.EGEOM*(SIZE2(P)+SIZE2(QA)+2)) THEN
        VG2LC=.FALSE.
        WRITE(IU,*) 'ERROR(*2LC) - The points P and QA must be'
        WRITE(IU,*) ' distinct'
      END IF 
            
C  Verify that P is not coincident with QB
      IF (DIST2(P,QB).LT.EGEOM*(SIZE2(P)+SIZE2(QB)+2)) THEN
        VG2LC=.FALSE.
        WRITE(IU,*) 'ERROR(*2LC) - The points P and QB must be'
        WRITE(IU,*) ' distinct'
      END IF       

C Warning if LPONEL=.FALSE. and the point P is too close to QA-QB
      IF (.NOT.LPONEL.AND.VG2LC) THEN
        IF (DIST2(P,QSTAR).LT.QLEN/10) THEN
          WRITE(IU,*) 'WARNING(*2LC) - If LPONEL=.FALSE. then if the'
          WRITE(IU,*) ' observaton point P is close to QA-QB'
          WRITE(IU,*) ' the expected error in solution likely to be'
          WRITE(IU,*) ' magnified.'
        END IF
      END IF
       
C Warning if LPONEL=.TRUE. and the point P is too close to QA or QB
      IF (LPONEL.AND.VG2LC) THEN
        IF ((DIST2(P,QA).LT.QLEN/10).OR.(DIST2(P,QB).LT.QLEN/10)) THEN
          WRITE(IU,*) 'WARNING(*2LC) - If LPONEL=.TRUE. then if the'
          WRITE(IU,*) ' observaton point P is close to QA or QB'
          WRITE(IU,*) ' the expected error in solution likely to be'
          WRITE(IU,*) ' magnified.'
        END IF
      END IF  
 
 
      END


