SUBROUTINE AREALSINK ( Xvert,Yvert,nvert,   
     &                         Xquery,Yquery,nquery,gamma0,
! output arrays:
     &                         PHI,Vx,Vy,A)
	!DEC$ ATTRIBUTES DLLEXPORT,ALIAS:'_arealsink_' :: ArealSink 

	IMPLICIT NONE
	
	INTEGER nvert,nquery
        DOUBLE PRECISION Xvert(nvert),Yvert(nvert)
        DOUBLE PRECISION Xquery(nquery),Yquery(nquery)
        DOUBLE PRECISION gamma0

C       Output arrays (double, allocated in C wrapper):

	DOUBLE PRECISION PHI(nquery),Vx(nquery),Vy(nquery)
        DOUBLE PRECISION A

C       Local variables

	COMPLEX*16 I / (0.0D0, 1.0D0) /
        INTEGER q,j
        COMPLEX*16 z,W,nuj,Zj,Cj,smallfj,Fj
	DOUBLE PRECISION Lj2,Yj,pi
 
        pi = 2.D0*AIMAG(LOG(I))
        A = 0.D0

C       Loop on all query points

        DO q = 1,nquery

          z = Xquery(q)+I*Yquery(q)
          PHI(q) = 0.D0
          W = 0.D0 + I*0.D0
          
          DO j=1,nvert

            IF(j.lt.nvert)THEN
              nuj = Xvert(j+1)-Xvert(j)+I*(Yvert(j+1)-Yvert(j))
            ELSE
              nuj = Xvert(1)-Xvert(j)+I*(Yvert(1)-Yvert(j))
            ENDIF

            Lj2 = ABS(nuj)
            Lj2 = Lj2*Lj2

! Compute local coordinate of current query point z,
! using Eq. (8.562) p.336 defining (Zj)+1 hence Zj 
            Zj = 2.0*(z-Xvert(j)-I*Yvert(j))/nuj - 1.0
            Yj = AIMAG(Zj);
! Constant Cj Eq. (8.571) p.337
            Cj = -2.0*LOG(0.5*nuj)
! Function f(Zj) (small f) Eq. (8.564) p.336
            smallfj = LOG((Zj-1.0)/(Zj+1.0))
! Function F(Zj) (capital F) Eq. (8.585) p.341
            Fj = Zj*smallfj-LOG(Zj-1.0)-LOG(Zj+1.0)+2.0+Cj
! Update sum in Eq. (8.576) p.338 defining W, with shortcuts:
!    (a) conj(Zj)-Zj = -2*I*Im{Zj} = -2*I*Yj
!    (b) F(Zj)+F(conj(Zj)) = F(Zj) + conj(F(Zj)) = 2*Re{F(Zj)}
            W = W + CONJG(nuj)*(-2.0*I*Yj*smallfj-2.0*DBLE(Fj))
! Update sum in Eq. (8.581) p.338 defining Phi, real-valued potential
            PHI(q) = PHI(q) + 2.0*Yj*Lj2*DBLE(Fj)

            IF(q.EQ.1)THEN
              A = A + 0.25*Yj*Lj2
            ENDIF

          ENDDO

! Apply multiplicating factor for Phi, Eq. (8.581)
          PHI(q) = -(gamma0/(32.*pi))*PHI(q)
! Apply multiplicating factor for W, Eq. (8.576)
          W = -(gamma0/(16.*pi*I))*W
! Get real-valued velocity components of V = conj(W)
          Vx(q) = DBLE(W)
          Vy(q) = -AIMAG(W)

        ENDDO

        END


        SUBROUTINE INPOLYGON( Xvert,Yvert,nvert,   
     &                        Xquery,Yquery,nquery,
! Output array
     &                        RES)
	!DEC$ ATTRIBUTES DLLEXPORT,ALIAS:'_inpolygon_' :: InPolygon 

	IMPLICIT NONE
	
	INTEGER nvert,nquery
        DOUBLE PRECISION Xvert(nvert),Yvert(nvert)
        DOUBLE PRECISION Xquery(nquery),Yquery(nquery)
!        INTEGER RES(nquery)       
        DOUBLE PRECISION RES(nquery)       

C       Local variables

	COMPLEX*16 I / (0.0D0, 1.0D0) /
        INTEGER q,j
        COMPLEX*16 z,zv1,zvnext,dz,dznext
        DOUBLE PRECISION pi,sumArg

        pi = 2.D0*AIMAG(LOG(I))
        zv1 = Xvert(1)+I*Yvert(1)

        DO q=1,nquery

          z = Xquery(q)+I*Yquery(q)
          dz = z-zv1
          sumArg = 0.D0

          DO j=1,nvert
              
            IF(j.LT.nvert)THEN
              zvnext = Xvert(j+1)+I*Yvert(j+1)
            ELSE
              zvnext = Xvert(1)+I*Yvert(1)
            ENDIF

            dznext = z-zvnext
!           Use of Eq. (8.552) p.334 with shortcut
!           ln((z-z_{j+1})/(z-z_j)) - conj(ln((z-z_{j+1})/(z-z_j))) = 2*Im{ln(dznext/dz)}
            sumArg = sumArg + AIMAG(LOG(dznext/dz))  
            dz = dznext
    
          ENDDO   
   
!          RES(q) = IDINT(sumArg/2.D0/pi)        
          RES(q) = sumArg/2.D0/pi      

        ENDDO
  
        END