anchor.diff:  END SUBROUTINE
   
   
  +!DIJ is the integral of dG/dn(XI,XJ) on the wake
  +!  XI is on the body
  +!  XJ is on the wake
  +double precision function DIJ(XI, XJM, XJP)
  +      IMPLICIT NONE
  +      real(kind=8), dimension(2), intent(IN) :: XI, XJM, XJP
  +      real(kind=8) :: XM, XP, YS, RM, RP
  +      real(kind=8) :: AYS, S
  +      real(kind=8), parameter :: PI = dble(4)*datan(1.D0)
  +      call LocFR(XI, XJM, XJP, XM, XP, YS, RM, RP)
  +      AYS = dabs(YS)
  +      S = dsign(1.D0,YS)
  +      DIJ = (datan2(S*XP,AYS) - datan2(S*XM,AYS))/(dble(2)*PI)
  +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)
  
positive.code: 
  !Calculate the matrix B and C
  !  N        Figure of elements  (IN)
  !  B        Matrix B            (OUT)
  !  C        Matrix C            (OUT)
  !  XNODE    Nodes vector        (IN)
  SUBROUTINE SRFMATBC(N, XNODE, B, C)
        IMPLICIT NONE
        INTEGER, INTENT(IN) :: N
  
negative.code:       INTEGER :: I, J, NOE
        DO I = 1, N
         !Midpoint of I element
         XI = (XNODE(NOE(N,I,0),:) + XNODE(NOE(N,I,1),:))/2.
         C(I,I) = 0.
         DO J = 1, N
          CALL BCIJ(XI, XNODE(NOE(N,J,0),:), XNODE(NOE(N,J,1),:),BIJ,CIJ)
          B(I,J) = BIJ
          IF (I .NE. J) THEN
            C(I,J) = CIJ
  