




      subroutine tanvecs(ijk,u,v,ru,rv,ruu,rvv,ruv)
c*******************************************************
c purpose: to evaluate the 2 tangent vectors and their derivatives
c          at a point.
c inputs:  u - value of u
c          v - value of v
c outputs: ru - u tangent vector (not normalized) (dr/du)
c          rv - v tangent vector (not normalized) (dr/dv)
c          ruu- 2nd order derivative of r(u,v) wrt u (d(dr/du)/du)
c          rvv- 2nd order derivative of r(u,v) wrt v (d(dr/dv)/dv)
c          ruv- 2nd order derivative of r(u,v) wrt u and v
c                note: d(dr/dv)/du = d(dr/du)/dv
c*******************************************************

      integer ijk
      real r(3),ru(3),rv(3),ruu(3),rvv(3),ruv(3)
      real u, v

      include 'ctparm.inc'
      include 'ctgeom.inc'
      include 'ctgeom2.inc'

c.. to speed things up, if the surface is a nurb, compute all
c   surface derivatives at once
      if(ntype(ijk) .eq. ep_pattypnurb) then
        call gethotans(iparms(ijk),u,v,r,ru,rv,ruu,rvv,ruv)
      else
        call tanvec_uandv(ijk,u,v,ru,rv)
        call hotanvec(ijk,u,v,ruu,rvv,ruv)
      endif

      return
      end

      subroutine hotanvec(ijk,u,v,ruu,rvv,ruv)
c*******************************************************************
c purpose:
c*******************************************************************

      include 'ctparm.inc'
      include 'ctgeom.inc'
      include 'ctgeom2.inc'

      integer ijk, i
      real ruu(3),rvv(3),ruv(3),x(3),z(3)
      real rdum(3),ru(3),rv(3)
      real u, v, r, rp, rpp, delri
      real a, b, c, al, rm, ack

c use the az & el angles for the circle orientations to fix
c the x- & z- axes...
      call uw_scamp(boxor(1,ijk),boxor(2,ijk),x,z)

      if( ntype(ijk) .eq. ep_pattypcirct ) then		!circular, t
        call get_abcl(ijk,a,b,c,al)
        r = sqrt(a*u**2 + b*u + c)
        rp = (a*u + b/2.) / r		    !note: srcrcs has special cases
        rpp = (a*c - b**2/4.) / r**3	    !for a=0, r small, etc. not yet
      	do i = 1,3			    !implemented here.
          ruu(i) = rpp * ( x(i)*cos(v) + z(i)*sin(v) )
          rvv(i) = -r  * ( x(i)*cos(v) + z(i)*sin(v) )
          ruv(i) = rp  * (-x(i)*sin(v) + z(i)*cos(v) )
        enddo
      elseif( ntype(ijk) .eq. ep_pattypcircpsi) then		!circular, psi
        rm  = psicoef(4,iparms(ijk))
        ack = psicoef(5,iparms(ijk))
        do i = 1,3
          delri = psidr(i,iparms(ijk))	  !(normalized delta-r)
          ruu(i) = ack*rm*cos(u)*delri - rm*sin(u)*( x(i)*cos(v) +
     %                                                 z(i)*sin(v) )
          rvv(i) = -rm*sin(u)*( x(i)*cos(v) + z(i)*sin(v) )
          ruv(i) =  rm*cos(u)*(-x(i)*sin(v) + z(i)*cos(v) )
        enddo
      elseif(ntype(ijk) .eq. ep_pattypbic) then         ! bicubic
        call bictanvecs(ijk,iparms(ijk),u,v,ruu,3)
        call bictanvecs(ijk,iparms(ijk),u,v,rvv,4)
        call bictanvecs(ijk,iparms(ijk),u,v,ruv,5)
      elseif(ntype(ijk) .eq. ep_pattypwire) then         ! wire
        do i=1,3
          ruu(i) = 0.
          rvv(i) = 0.
          ruv(i) = 0.
        end do
      elseif(ntype(ijk) .eq. ep_pattypnurb) then	!nurb
        call gethotans(iparms(ijk), u, v, rdum, ru, rv, ruu, rvv, ruv)
      elseif(ntype(ijk) .eq. ep_pattypflat) then
        do i=1,3
          ruu(i) = 0.
          rvv(i) = 0.
          ruv(i) = 0.
        end do
      else
        write(6,*) ' *** internal error *** in hotanvec. '
        write(6,*) '     -  patch/wire #',ijk,' is not a known type.'
      endif
  
      return
      end


      subroutine tridata(ijk,i2,i3,i4,n1,n2,n3)
c******************************************************
c purpose: to get the edges & nodes if the ijk-th patch.
c******************************************************

      include 'ctparm.inc'
      include 'ctgeom.inc'

      integer ijk, i2, i3, i4, n1, n2, n3

c-- get edges
      i2 = nbound(1,ijk)
      i3 = nbound(2,ijk)
      i4 = nbound(3,ijk)
c-- get nodes
      n1 = nbound(4,ijk)
      n2 = nbound(5,ijk)
      n3 = nbound(6,ijk)

      return
      end



      subroutine ac_trans(ipat,xsi,eta,u,v)

c purpose: to convert a point in the local area coordinates (xsi,eta)
c for patch ipat into (u,v) coordinates.
    
      include 'ctparm.inc'
      include 'ctgeom2.inc'

      integer ipat
      real xsi, eta, u, v

      u = tzpt(6,ipat) + xsi * ( tzpt(2,ipat) - tzpt(6,ipat) ) +
     %                   eta * ( tzpt(4,ipat) - tzpt(6,ipat) )

      v = tzpt(5,ipat) + xsi * ( tzpt(1,ipat) - tzpt(5,ipat) ) +
     %                   eta * ( tzpt(3,ipat) - tzpt(5,ipat) )

      return
      end



       subroutine tanvec_uandv(ijk, u, v, dru, drv)
c---------------------------------------------------------------------
c purpose: to compute the surface tangent vector along u (dru) for
c          patch ijk at point (u,v). the actual equation for
c          the tangent vector depends on how the surface was
c          parameterized.
c---------------------------------------------------------------------

      include 'ctparm.inc'
      include 'ctgeom.inc'
      include 'ctgeom2.inc'
      include 'ctgeomwire.inc'

      integer ijk, i, k
      real dru(3),drv(3),x(3),z(3), u, v
      real rm, ack, delri, rp, r
      real r1(3), r2(3), r3(3), rjnk(3)
      real a, b, c, al


c use the az & el angles for the circle orientations to fix
c the x- & z- axes...
       call uw_scamp(boxor(1,ijk),boxor(2,ijk),x,z)

       if (ntype(ijk).eq.ep_pattypcircpsi) then		! psi - phi parameters
         rm = psicoef(4,iparms(ijk)) 
         ack = psicoef(5,iparms(ijk))
         do i = 1,3
          delri = psidr(i,iparms(ijk))
          dru(i) = rm*ack*sin(u)*delri + rm*cos(u)*
     %				( x(i)*cos(v) + z(i)*sin(v) )
          drv(i) = rm * sin(u) * ( -x(i)*sin(v) + z(i)*cos(v) )
         enddo

       elseif (ntype(ijk).eq.ep_pattypcirct) then	! t - phi parameters
         call get_abcl(ijk,a,b,c,al)
         r = sqrt( a*u*u + b*u + c )
         rp = ( a*u + (b/2.) ) / sqrt( a*u**2 + b*u + c )
         do i = 1,3
          delri = tr0(i+3,iparms(ijk))
          dru(i) = al*delri + rp*( x(i)*cos(v) + z(i)*sin(v) )
          drv(i) = r * (-x(i)*sin(v) + z(i)*cos(v) )
         enddo

       elseif(ntype(ijk) .eq. ep_pattypbic) then
        call bictanvecs(ijk,iparms(ijk),u,v,dru,1)
        call bictanvecs(ijk,iparms(ijk),u,v,drv,2)

       elseif(ntype(ijk) .eq. ep_pattypwire) then           ! wire
         do i=1,3
           dru(i) = datnod(i,nconn_wire(2,ijk)) - 
     &              datnod(i,nconn_wire(1,ijk)) 
         end do

       elseif(ntype(ijk) .eq. ep_pattypnurb) then	!nurb
         call gettans(iparms(ijk), u, v, rjnk, dru, drv)

       elseif(ntype(ijk) .eq. ep_pattypflat) then	!flat
         do k=1,3
           r1(k) = datnod(k,nbound(4,ijk))
           r2(k) = datnod(k,nbound(5,ijk))
           r3(k) = datnod(k,nbound(6,ijk))
         end do
         call vecdif(r2,r3,dru)
         call vecdif(r1,r3,drv)

      else
        write(6,*) ' *** internal error *** in tanvec_uandv. '
        write(6,*) '     -  patch/wire #',ijk,' is not a known type.'
       endif


       return
       end



       subroutine get_abcl(ipat,a,b,c,al)
c-----------------------------------------------------------------------
c purpose: to get the values of a, b, c & l for patch ipat. this
c          assumes t-parameterization.
c-----------------------------------------------------------------------

      include 'ctparm.inc'
      include 'ctgeom2.inc'

      integer ipat
      real a, b, c, al
c note: a,b,c & al are real*4

       a = tcoef(1,iparms(ipat))
       b = tcoef(2,iparms(ipat))
       c = tcoef(3,iparms(ipat))
       al= tcoef(4,iparms(ipat))

       return
       end

       real function get_r(ipat_s,u_s,v_s,ipat_f,u_f,v_f)
c---------------------------------------------------------------------
c purpose: to compute the magnitude of the distance vector between the 
c          source & field points. 
c inputs: source patch # ipat_s with source point (u_s,v_s) (local coord)
c         field  patch # ipat_f with field  point (u_f,v_f) (local coord)
c---------------------------------------------------------------------

      integer ipat_s, ipat_f
      real u_s, v_s, u_f, v_f
      real r_s(3),r_f(3)

c find the vectors (in cartesian coord) to the field & source points
       call rvec(ipat_s, u_s, v_s, r_s)
       call rvec(ipat_f, u_f, v_f, r_f)

c compute r = |r-r'|
       get_r = sqrt( (r_f(1) - r_s(1))**2 + (r_f(2) - r_s(2))**2 +
     %               (r_f(3) - r_s(3))**2 )

       return
       end

       real function get_rs(ipat,u_s,v_s,u_f,v_f)
c---------------------------------------------------------------------
c purpose: to compute the approximate distance vector between the 
c          source & field points. 
c assumption: both the source point (u_s,v_s) and the field point
c             (u_f,v_f) lie in the same patch (ipat)
c---------------------------------------------------------------------

      integer ipat
      real u_s, u_f, v_s, v_f
      real dru(3),drv(3)
      real tmpa, tmpb, tmpc
      real dotp

c find the tangent vectors (in cartesian coord) at the field point
       call tanvec_uandv(ipat, u_f, v_f, dru, drv)

c compute rs 
       tmpa = dotp(dru,dru) * (u_f - u_s)**2
       tmpb = dotp(drv,drv) * (v_f - v_s)**2
       tmpc = 2. * dotp(dru,drv) * (u_f - u_s) * (v_f - v_s)

       get_rs = sqrt( tmpa + tmpb + tmpc )

       return
       end



       subroutine uw_scamp(el,az,u_sc,w_sc)
c----------------------------------------------------------------------
c purpose:
c----------------------------------------------------------------------
      real el, az
      real u_sc(3),w_sc(3)

       u_sc(1) = -sin(el)*sin(az)
       u_sc(2) = -sin(el)*cos(az)
       u_sc(3) =  cos(el)
 
       w_sc(1) = cos(az)
       w_sc(2) =-sin(az)
       w_sc(3) = 0.

       return
       end

       subroutine rvec(ipat,u,v,r)
c----------------------------------------------------------------------
c purpose: to compute the vector r (in cartesian coord) to the point
c          (u,v) in patch ipat.
c----------------------------------------------------------------------

      include 'ctparm.inc'
      include 'ctgeom.inc'
      include 'ctgeom2.inc'
      include 'ctgeomwire.inc'
      real r(3),x(3),z(3)
      real u, v, rm, ack, vr0i, delri
      real a, b, c, al, rr
      real r1(3), r2(3), r3(3)
      integer ipat, i, k

       call uw_scamp(boxor(1,ipat),boxor(2,ipat),x,z)

       if (ntype(ipat).eq.ep_pattypcircpsi) then		 ! psi - phi parameters
         rm  = psicoef(4,iparms(ipat)) 
         ack = psicoef(5,iparms(ipat))
         do i = 1,3
           vr0i = psicoef(i,iparms(ipat))
           delri = psidr(i,iparms(ipat))
           r(i) = vr0i - rm*ack*cos(u)*delri + rm*sin(u)*
     %		      ( cos(v)*x(i) + sin(v)*z(i) )
         enddo
       elseif(ntype(ipat).eq.ep_pattypcirct) then	      	! t - phi parameters
         call get_abcl(ipat,a,b,c,al)
         rr = sqrt( a*u**2 + b*u + c )
         do i = 1,3
           vr0i = tr0(i,iparms(ipat))
           delri = tr0(3+i,iparms(ipat))
           r(i) = vr0i + u*al*delri + rr*(cos(v)*x(i) + sin(v)*z(i))
         enddo
       elseif(ntype(ipat) .eq. ep_pattypbic) then	   ! bicubic
         call bictanvecs(ipat,iparms(ipat),u,v,r,0)
       elseif(ntype(ipat) .eq. ep_pattypwire) then         ! wire
         do i=1,3
           r(i) = datnod(i,nconn_wire(1,ipat)) + 
     &        u*( datnod(i,nconn_wire(2,ipat)) - 
     *            datnod(i,nconn_wire(1,ipat)) )
         end do
       elseif(ntype(ipat) .eq. ep_pattypnurb) then 	   ! nurb
         call getvertex(iparms(ipat), u, v, r)
       elseif(ntype(ipat) .eq. ep_pattypflat) then	   ! flat
         do k=1,3
           r1(k) = datnod(k,nbound(4,ipat))
           r2(k) = datnod(k,nbound(5,ipat))
           r3(k) = datnod(k,nbound(6,ipat))
         end do
         do k=1,3
           r(k) = r3(k) + u * (r2(k)-r3(k)) + v * (r1(k)-r3(k))
         end do
       else
         write(6,*) ' *** internal error *** in rvec. '
         write(6,*) '     -  patch/wire #',ijk,' is not a known type.'
       endif

       return
       end
