c    **************************************************************
c    *                                                            *
c    *	 this program computes velocity and pressure coefficient  *
c    *	 around an airfoil, whose contour is approximated by	  *
c    *	 vortex panels of linearly varying strength.              *
c    *                                                            *
c    **************************************************************
c
c
      parameter (npts = 13)
      dimension xb(npts),yb(npts),x(npts),y(npts-1),s(npts-1),
     .          sine(npts-1),cosine(npts-1),theta(npts-1),v(npts-1),
     .          cp(npts-1),gama(npts),rhs(npts),cn1(npts-1,npts-1),
     .          cn2(npts-1,npts-1),ct1(npts-1,npts-1),
     .          ct2(npts-1,npts-1),an(npts,npts),at(npts-1,npts),
     .          indx(npts)
c
c    Specify coordinates (xb,yb) of boundary points on airfoil
c    surface.  The last point coincides with the first.
c
c
c  Nominal data
c
      data xb
     1 /1.0000,.933, .750, .500, .250, .067, 0. , .067, .250, .500,
     2	 .750, .933, 1.00/
      data yb
     1 /.0000,-.005,-.017,-.033,-.042,-.033,0.,.045,.076,
     2 .072,.044,.013, 0./
c
c 18 panel Wortmann airfoil data
c
c     data xb/
c    .  0.1000E+01,  0.8889E+00,  0.7778E+00,  0.6667E+00,  0.5556E+00,
c    .  0.4444E+00,  0.3333E+00,  0.2222E+00,  0.1111E+00,  0.0000E+00,
c    .  0.1111E+00,  0.2222E+00,  0.3333E+00,  0.4444E+00,  0.5556E+00,
c    .  0.6667E+00,  0.7778E+00,  0.8889E+00,  0.1000E+01/
c     data yb/
c    .  0.0000E+00,  0.2364E-01,  0.2745E-01,  0.2148E-01,  0.9967E-02,
c    . -0.3737E-02, -0.1508E-01, -0.2161E-01, -0.2209E-01,  0.0000E+00,
c    .  0.8570E-01,  0.1112E+00,  0.1211E+00,  0.1199E+00,  0.1095E+00,
c    .  0.9087E-01,  0.6605E-01,  0.3765E-01,  0.0000E+00/
c
      m = npts - 1
      mp1 = npts
      pi = 4.*atan(1.)
      alpha = 8.*pi/180.
c
c    Coordinates (x,y) of control point and panel length s are
c    computed for each of the vortex panels.   rhs represents
c    the right-hand side of eq. (47).
c
      do 1 i = 1,m
         ip1 = i+1
         x(i) = .5*(xb(i)+xb(ip1))
         y(i) = .5*(yb(i)+yb(ip1))
         s(i) = sqrt((xb(ip1)-xb(i))**2+(yb(ip1)-yb(i))**2)
         theta(i) = atan2(yb(ip1)-yb(i),xb(ip1)-xb(i))
         sine(i) = sin(theta(i))
         cosine(i) = cos(theta(i))
    1    rhs(i) = sin(theta(i)-alpha)
      do 3 i = 1,m
         do 3 j = 1,m
            if (i.eq.j) go to 2
            a = -(x(i)-xb(j))*cosine(j)-(y(i)-yb(j))*sine(j)
            b = (x(i)-xb(j))**2+(y(i)-yb(j))**2
            c = sin(theta(i)-theta(j))
            d = cos(theta(i)-theta(j))
            e = (x(i)-xb(j))*sine(j)-(y(i)-yb(j))*cosine(j)
            f = alog(1.+s(j)*(s(j)+2.*a)/b)
            g = atan2(e*s(j),b+a*s(j))
            p = (x(i)-xb(j))*sin(theta(i)-2.*theta(j))
     .        +(y(i)-yb(j))*cos(theta(i)-2.*theta(j))
            q = (x(i)-xb(j))*cos(theta(i)-2.*theta(j))
     .        -(y(i)-yb(j))*sin(theta(i)-2.*theta(j))
            cn2(i,j) = d+.5*q*f/s(j)-(a*c+d*e)*g/s(j)
            cn1(i,j) = .5*d*f+c*g-cn2(i,j)
            ct2(i,j) = c+.5*p*f/s(j)+(a*d-c*e)*g/s(j)
            ct1(i,j) = .5*c*f-d*g-ct2(i,j)
            go to 3
    2       cn1(i,j) = -1.0
            cn2(i,j) = 1.0
            ct1(i,j) = .5*pi
            ct2(i,j) = .5*pi
    3 continue
c
c    compute influence coefficients in eqs. (47) and (49),
c    respectively.
c
      do 4 i = 1,m
         an(i,1) = cn1(i,1)
         an(i,mp1) = cn2(i,m)
         at(i,1) = ct1(i,1)
         at(i,mp1) = ct2(i,m)
         do 4 j = 2,m
            an(i,j) = cn1(i,j)+cn2(i,j-1)
    4       at(i,j) = ct1(i,j)+ct2(i,j-1)
      an(mp1,1) = 1.0
      an(mp1,mp1) = 1.0
      do 5 j = 2,m
    5    an(mp1,j) = 0.
      rhs(mp1) = 0.
c
c    Solve eq.(47) for dimensionless strengths gama using
c    LU decomposition.  Then compute and print dimensionless
c    velocity and pressure coefficient at control points.
c
      write(6,55)m,alpha*180./pi
   55 format('1'/30x,i2,' panels,  alpha =',f5.2,' deg'/)
      write(6,6)
    6 format(11x,'i',4x,'x(i)',4x,'y(i)',4x,'theta(i)',
     1       3x,'s(i)',3x,'gama(i)',3x,'v(i)',6x,'cp(i)'/10x,
     2       '---',3x,'----',4x,'----',4x,'--------',3x,'----',
     3       3x,'-------',3x,'----',6x,'-----')
c
c    Solve the system of equations with an LU decomposition and
c    back substitution
c
      call ludcmp(an,mp1,mp1,indx,d)
      call lubksb(an,mp1,mp1,indx,rhs)
c
      do 105 i = 1 , mp1
  105    gama(i) = rhs(i)
c
      do 8 i = 1,m
         v(i) = cos(theta(i)-alpha)
         do 7 j = 1,mp1
            v(i) = v(i)+at(i,j)*gama(j)
    7       cp(i) = 1.-v(i)**2
    8    write(6,9)i,x(i),y(i),theta(i),s(i),gama(i),v(i),cp(i)
c
    9 format(10x,i2,f8.4,f9.4,f10.4,f8.4,2f9.4,f10.4)
      write(6,10)mp1,gama(mp1)
   10 format(10x,i2,35x,f9.4/)
c
c
      cfx=0.0
      cfy=0.0
      cm=0.0
c
      do 95 i=1,m
         dx=xb(i+1)-xb(i)
         dy=yb(i+1)-yb(i)
         cfx=cfx+cp(i)*dy
         cfy=cfy-cp(i)*dx
         cm=cm+cp(i)*(dx*x(i)+dy*y(i))
   95 continue
c
      cd=cfx*cos(alpha)+cfy*sin(alpha)
      cl=cfy*cos(alpha)-cfx*sin(alpha)
c
      write(6,110)cd,cl,cm
  110 format('   cd =',f8.5,'   cl =',f8.5,'   cm =',f8.5)
c
c
      stop
      end
c
c
c                  lu decomposition and backsubstitution
c                      taken from "numerical recipes"
c                              june 6, 1990
c
c ......................................................................
c  "ludcmp()"; given an "n" by "n" matrix "a", with physical dimension
c  "np", this routine replaces it by the lu decomposition of a rowwise
c  permutation of itself.  "a" and "n" are input.  "a" is output; "indx"
c  is an output vector which records the row permutation effected by the
c  partial pivoting; "d" is output as + or - 1 depending on whether the
c  number of row interchanges was even or odd, respectively.  this
c  routine is used in combination with "lubksb()" to solve linear
c  equations or invert a matrix.
c
      subroutine ludcmp(a, n, np, indx, d)
      parameter (nmax = 100, tiny = 1.0e-20)
      dimension a(np, np), indx(n), vv(nmax)
      d = 1.0
      do 12 i = 1, n
        aamax = 0.0
        do 11 j = 1, n
           if (abs(a(i,j)) .gt. aamax) aamax = abs(a(i,j))
   11   continue
        if (aamax .eq. 0.0) pause 'singular matrix'
           vv(i) = 1.0 / aamax
   12 continue
      do 19 j = 1, n
         do 14 i = 1, j - 1
            sum = a(i,j)
            do 13 k = 1, i - 1
               sum = sum - a(i,k) * a(k,j)
   13       continue
         a(i,j) = sum
   14    continue
         aamax = 0.0
         do 16 i = j, n
            sum = a(i,j)
            do 15 k = 1, j - 1
                  sum = sum - a(i,k) * a(k,j)
   15       continue
            a(i,j) = sum
            dum = vv(i) * abs(sum)
            if (dum .ge. aamax) then
                imax = i
                aamax = dum
                endif
   16    continue
         if (j .ne. imax) then
             do 17 k = 1, n
                dum = a(imax,k)
                a(imax,k) = a(j,k)
                a(j,k) = dum
   17        continue
             d = -d
             vv(imax) = vv(j)
         endif
         indx(j) = imax
         if (a(j,j) .eq. 0.0) a(j,j) = tiny
         if (j .ne. n) then
             dum = 1.0 / a(j,j)
             do 18 i = j + 1, n
                a(i,j) = a(i,j) * dum
   18        continue
         endif
   19 continue
      return
      end
c
c .....................................................................
c  "lubksb()" solves the set of "n" linear equations [a] * {x} = {b}.
c  here "a" is input, not as the matrix "a" but rather as its lu
c  decomposition, determined by the routine "ludcmp()".  "indx" is
c  input as the permutation vector returned by "ludcmp()".  "b" is
c  input as the right-hand side vector, and returns with the solution
c  vector {x}.  "a", "n", "np" and "indx" are not modified by this
c  routine and can be left in place for successive calls with different
c  right-hand sides "b".  this routine takes into account the possibility
c  that "b" will begin with many zero elements, so it is efficient for
c  use in matrix inversion.
c
      subroutine lubksb(a, n, np, indx, b)
      dimension a(np, np), indx(n), b(n)
      ii = 0
      do 12 i = 1, n
         ll = indx(i)
         sum = b(ll)
         b(ll) = b(i)
         if (ii .ne. 0) then
            do 11 j = ii, i - 1
               sum = sum - a(i,j) * b(j)
   11       continue
         else if (sum .ne. 0.0) then
            ii = i
         endif
         b(i) = sum
   12 continue
      do 14 i = n, 1, -1
         sum = b(i)
         if (i .lt. n) then
             do 13 j = i + 1, n
                sum = sum - a(i,j) * b(j)
   13        continue
         endif
         b(i) = sum / a(i,i)
   14 continue
      return
      end