      subroutine sigint(ijk,u_f,v_f,rint)
c*************************************************************************
c purpose: to evaluate the singular integral used in both the vector & scalar 
c potential expressions for a curved patch ijk:
c
c     i = int(1/rs) du dv 
c
c at field point u_f,v_f.
c modification 11/91: if u & v are not orthogonal, perform a change of 
c variable to make them so.
c*************************************************************************

      include 'ctparm.inc'
      include 'ctgeom2.inc'
       real dru(3),drv(3)
       data eps/1.e-5/		!roundoff error

c compute the tangent vectors at the field point
       call tanvec_uandv(ijk,u_f,v_f,dru,drv)

c determine whether u & v are orthogonal. if they aren't, then need to 
c do a coordinate transform to form 2 orthogonal vectors
       a = dotp(dru,dru)
       b = dotp(drv,drv)
       c = 2.*dotp(dru,drv)
       if ( abs(c) .le. eps ) then
         a = sqrt(a)
         b = sqrt(b)
c change of coordinates for 3 vertices and field point
         x1 = a * tzpt(2,ijk)
         y1 = b * tzpt(1,ijk)
         x2 = a * tzpt(4,ijk)
         y2 = b * tzpt(3,ijk)
         x3 = a * tzpt(6,ijk)
         y3 = b * tzpt(5,ijk)
         x =  a * u_f
         y =  b * v_f
         area = abs(x1*(y2-y3) - y1*(x2-x3) + x2*y3 - y2*x3)/2.
         call intglcha(x1,y1,0.,x2,y2,0.,x3,y3,0.,x,y,0.,pot,area)
         rint = pot / (a*b)

       else			

c description of the treatment of non-orthogonal (u,v) is given in
c note by d.wilkes dated nov 12, 1991
         if (abs(b).lt.eps) then	!b=0 (can't also have a=0)
            bcap = 0.
            acap = -c/(2.*a)
            sqr_alp = sqrt(a)
            sqr_bet = sqrt(abs( b - c**2/(4.*a) ))
         else
            acap = 0.
            bcap = -c/(2.*b)
            sqr_bet = sqrt(b)
            sqr_alp = sqrt(abs( a - c**2/(4.*b) ))
         endif
c change of coordinates for 3 vertices
         u1 = tzpt(2,ijk)
         v1 = tzpt(1,ijk)
         u2 = tzpt(4,ijk)
         v2 = tzpt(3,ijk)
         u3 = tzpt(6,ijk)
         v3 = tzpt(5,ijk)
         x1 = sqr_alp * ( (u_f - u1) - acap*(v_f - v1) )
         x2 = sqr_alp * ( (u_f - u2) - acap*(v_f - v2) )
         x3 = sqr_alp * ( (u_f - u3) - acap*(v_f - v3) )
         y1 = sqr_bet * ( (v_f - v1) - bcap*(u_f - u1) )
         y2 = sqr_bet * ( (v_f - v2) - bcap*(u_f - u2) )
         y3 = sqr_bet * ( (v_f - v3) - bcap*(u_f - u3) )
         x =  0.
         y =  0.
         area = abs(x1*(y2-y3) - y1*(x2-x3) + x2*y3 - y2*x3)/2.
         call intglcha(x1,y1,0.,x2,y2,0.,x3,y3,0.,x,y,0.,pot,area)
         rint = pot / (sqr_alp*sqr_bet)

       endif
 
      return
      end

