anchor.diff:        DO I = 1, NX
          XI = X(I,:)
          DO J = 1, N
  -        B(I,J) = BIJ(XI, XNODE(NODOFEL(N,J,1),:), XNODE(NODOFEL(N,J,2),:))
  -        C(I,J) = CIJ(XI, XNODE(NODOFEL(N,J,1),:), XNODE(NODOFEL(N,J,2),:))
  +        CALL BCIJ(XI, XNODE(NODOFEL(N,J,1),:), XNODE(NODOFEL(N,J,2),:), BIJ, CIJ)
  +        B(I,J) = BIJ
  +        C(I,J) = CIJ
          END DO
         END DO
   END SUBROUTINE
  
positive.code: SUBROUTINE BCONDVEL(N, XNODE, U, CHI)
        IMPLICIT NONE
        INTEGER, INTENT(IN) :: N
        REAL, DIMENSION(N,2), INTENT(IN) :: XNODE
        REAL, DIMENSION(2), INTENT(IN) :: U
        REAL, DIMENSION(N), INTENT(OUT) :: CHI
        REAL, DIMENSION(2) :: T
        INTEGER :: I, NODOFEL
        DO I = 1, N
         T = XNODE(NODOFEL(N,I,2),:) - XNODE(NODOFEL(N,I,1),:)
         T = T/SQRT(T(1)**2 + T(2)**2)
         CHI(I) = U(1)*T(2) - U(2)*T(1)
        END DO
  
negative.code: END FUNCTION
  
  
  !BIJ is the integral of G(XI,XJ) from XJM to XJP, nodes of the J element
  !CIJ is the integral of dG/dn(XI,XJ) from XJM to XJP, nodes of the J element
  SUBROUTINE BCIJ(XI, XJM, XJP, BIJ, CIJ)
        IMPLICIT NONE
        REAL, DIMENSION(2), INTENT(IN) :: XI, XJM, XJP
        REAL, INTENT(OUT) :: BIJ, CIJ
        REAL, DIMENSION(2) :: EX, VXP, VXM
        REAL :: XM, XP, YS, RM, RP
        REAL, PARAMETER :: PI = 4.*ATAN(1.)
        EX = XJP - XJM
        EX = EX/SQRT(EX(1)**2 + EX(2)**2)
  !     Vector X_+
        VXP = (XI - XJP)
  !     Vector X_-
        VXM = (XI - XJM)
        XP = -VXP(1)*EX(1) - VXP(2)*EX(2)
        XM = -VXM(1)*EX(1) - VXM(2)*EX(2)
        YS = -VXP(1)*EX(2) + VXP(2)*EX(1)
        RP = SQRT(XP**2 + YS**2)
        RM = SQRT(XM**2 + YS**2)
        BIJ = (XP*LOG(RP) - XP + YS*ATAN2(XP,YS) - XM*LOG(RM) + XM - YS*ATAN2(XM,YS))/(2.*PI)
        CIJ = (ATAN2(XP,YS) - ATAN2(XM,YS))/(2.*PI)      
  END SUBROUTINE
  
  
  !Calculate the matrix B and C
  