program flow
*
* driver for kirchhoff program kirch.
* computes free-streamline flow over a polygonal obstacle.
* l. n. trefethen, summer 1984
*
      implicit double precision (a-b,d-h,o-w,y)
      implicit complex*16 (c,x,z)
      common /dual/ sgx
      dimension z(20),x(20),betam(20),qwork(400)
      dimension x2(20),z2(20),betam2(20),qwork2(400)
      data zi,zoffst /(0.d0,1.d0),(2.d0,0.d0)/
      sgx = 1.
*
* input control parameters:
    1 print '(''n, nq?'')'
      read *,n,nq
      do 2 k = 1,n
        print '(''z('',i1,'') ?'')', k
        read *, aa,bb
    2   z(k) = dcmplx(aa,bb)
      tol = 10.**(-min(12,nq+2))
*
* solve mapping problem:
      call count0
      call qinitx(n,z,betam,nq,qwork)
      call ksolv(0,tol,n,c,x,z,betam,nq,qwork)
      aw = abs(c)
      call count
*
* print results:
      write (10,'(''w'',4f10.5)') -3.,-3.,3.,3.
      print '(''stagnation point:'',2f10.7)', z(n+1)
      write (10,'(''m -2.5 -2.3''/''t 17 stagnation point:'')')
      write (10,'(''m -1.  -2.3''/''t 20 '',2f10.6)') z(n+1)
      print '(''w:'',f10.7)', aw
      write (10,'(''m -2.5 -2.6''/''t 17                w:'')')
      write (10,'(''m -1.  -2.6''/''t 10 '',f10.6)') aw
*
* force computation 1 -- along obstacle:
      sgx = -1.
      call qinitx(n,z,betam2,nq,qwork2)
      x2(1) = x(1)
      z2(1) = z(1)
      write (10,101) z2(1) + zoffst
      kmq = 1
      do 12 k = 2,n
        km = k-1
        if (betam2(k).ge.-.9999) then
*         --- non-sharp edge ---
          x2(k) = x(k)
          z2(k) = zx(x2(k),k,x2(km),z2(km),kmq,n,c,x,betam2,nq,qwork2)
          write (10,102) z2(k) + zoffst
          kmq = k
        else
*         --- sharp edge ---
          scale = min( abs(x(k)-x(km)) , abs(x(k+1)-x(k)) )
          xx = x(k) - 1.d-2*scale
          zz = zx(xx,0,x2(km),z2(km),kmq,n,c,x,betam2,nq,qwork2)
          write (10,102) zz + zoffst
          xx2 = x(k) - .1*zi*scale
          zz = zx(xx2,0,xx,zz,0,n,c,x,betam2,nq,qwork2)
          x2(k) = x(k) + 1.d-2*scale
          z2(k) = zx(x2(k),0,xx2,zz,0,n,c,x,betam2,nq,qwork2)
          write (10,101) z2(k) + zoffst
          kmq = 0
          endif
   12   continue
      zf1 = zi * (z2(n)-z(n))
*
* force computation 2 -- via point x = -i :
      ztmp = zx(-zi,0,x(1),z(1),1,n,c,x,betam2,nq,qwork2)
      ztmp = zx(x(n),n,-zi,ztmp,0,n,c,x,betam2,nq,qwork2)
      zf2 = zi * (ztmp-z(n))
*
* force computation 3 -- via point x = 1.1-i :
      zz = (1.1d0,-1.d0)
      ztmp = zx(zz,0,x(1),z(1),1,n,c,x,betam2,nq,qwork2)
      ztmp = zx(x(n),n,zz,ztmp,0,n,c,x,betam2,nq,qwork2)
      zf3 = zi * (ztmp-z(n))
*
* print forces; all three sets of numbers should agree:
      print '(''  drag ='',3f12.6)', dreal(zf1),dreal(zf2),dreal(zf3)
      print '(''  lift ='',3f12.6)', dimag(zf1),dimag(zf2),dimag(zf3)
      print '('' total ='',3f12.6)', abs(zf1),abs(zf2),abs(zf3)
      write (10,'(''m -2.5 -2.''/''t 24 drag, lift, total force:'')')
      write (10,'(''m 0. -2.''/''t 30 '',3f10.6)') dreal(zf2),
     &  dimag(zf2), abs(zf2)
      sgx = 1.
*
* plot free streamlines:
   23 write (10,101) z(1)
      write (10,102) (z(k),k=2,n)
      print '('' no. of streamlines ?'')'
      read *, nstr
      if (nstr.eq.0) stop
      space = .1d0*(2.d0/aw)
      npts = 30
      do 4 nl = 0,1
        write (10,101) z(n-nl*(n-1))
        do 3 k = 0,npts
          xx = (-1.)**nl * (1. + k**3*sqrt(32.*space)/npts**3)
          call nearx(xx,x0,z0,k0,n,x,z,betam)
          zz = zx(xx,0,x0,z0,k0,n,c,x,betam,nq,qwork)
    3     write (10,102) zz
    4   continue
      print '(''finished free streamlines'')'
*
* plot internal streamlines:
      do 6 nl = -nstr,nstr
        do 5 k = -npts,npts
          if (nl.eq.0.and.k.gt.0) goto 6
          xx = sqrt(dcmplx(space*k**3*54./npts**3,nl*space))
          if (dimag(xx).lt.0.d0) xx = -xx
          xx = xx + x(n+1)
          call nearx(xx,x0,z0,k0,n,x,z,betam)
          zz = zx(xx,0,x0,z0,k0,n,c,x,betam,nq,qwork)
          if (k.eq.-npts)  write (10,101) zz
          write (10,102) zz
    5     continue
    6   continue
      print '(''finished internal streamlines'')'
*
* plot upstream equipotential lines:
      do 8 nl = -16,-1
        klim = 2*nstr
        if (nl.lt.-7) klim = nstr
        do 7 k = -klim,klim
          xx = sqrt(dcmplx(2.*nl*space,k*nstr*space/klim))
          if (dimag(xx).lt.0.d0) xx = -xx
          xx = xx + x(n+1)
          call nearx(xx,x0,z0,k0,n,x,z,betam)
          zz = zx(xx,0,x0,z0,k0,n,c,x,betam,nq,qwork)
          if (k.eq.-klim)  write (10,101) zz
          write (10,102) zz
    7     continue
    8   continue
      print '(''finished upstream equipotential lines'')'
*
* plot downstream equipotential lines:
      do 11 nn = 1,2
        do 10 nl = 0,27
          klim = 2*nstr
          if (nl.eq.0) klim = 10*klim
          if (nl.gt.5) klim = nstr
          do 9 k = 0,klim
            xx = sqrt(dcmplx(2.*nl*space,k*(-1)**nn*nstr*space/klim))
            if (dimag(xx).lt.0.d0) xx = -xx
            if (nn.eq.1.and.k.eq.0) xx = -xx
            xx = xx + x(n+1)
            call nearx(xx,x0,z0,k0,n,x,z,betam)
            zz = zx(xx,0,x0,z0,k0,n,c,x,betam,nq,qwork)
            if (k.eq.0)  write (10,101) zz
            write (10,102) zz
    9       continue
   10     continue
   11   continue
      print '(''finished downstream equipotential lines'')'
      call count
  101 format ('m ',2f9.4)
  102 format ('c ',2f9.4)
      end
*     kirch1
*
* computes the conformal map of the upper half plane (x) onto an
* asymmetrical free-streamline flow domain (z) of kirchhoff type.
* x is the square root of the velocity potential.
* the obstacle is composed of n-1 straight line segments with
* vertices z(k) indexed from 1 to n.  z(n+1) is the stagnation point.
* see monakhov, "boundary-value problems with free boundaries...", iv.1.
*
* adapted from scpack.  l. n. trefethen, june 1984.


      subroutine qinitx(n,z,betam,nq,qwork)
*
* computes nodes and weights for gauss-jacobi quadrature
*
      implicit double precision (a-b,d-h,o-w,y)
      implicit complex*16(c,x,z)
      dimension qwork(1),z(1),betam(1)
      call anglex(n,z,betam)
      iscr = nq*(2*n+2) + 1
      do 1 k = 1,n+1
        inodes = nq*(k-1) + 1
        iwts = nq*(n+k) + 1
    1   if (betam(k).ge.-.9999) call gaussj(nq,0.d0,
     &    betam(k),qwork(iscr),qwork(inodes),qwork(iwts))
      return
      end


      subroutine ksolv(iprint,tol,n,c,x,z,betam,nq,qwork)
*
* directs the solution of the modified schwarz-christoffel parameter
* problem for kirchhoff flow over a polygonal obstacle
*
      implicit double precision (a-b,d-h,o-w,y)
      implicit complex*16(c,x,z)
      common /param1/ nq2,c2,x2(20),z2(20),qwork2(460),betam2(20)
      dimension x(1),z(1),betam(1),qwork(1)
      dimension ajinv(20,20),scr(900),fval(19),y(19)
      external kfun
      np = n+1
      nm = n-1
      nm2 = n-2
      x(1) = -1.
      x(n) =  1.
*
* initial guess and ns01a control parameters:
      do 1 k = 2,nm
    1   y(k-1) = 0.
      dstep = 1.d-6
      dmax = 300.
      maxfun = 30*nm
*
* copy input data to /param1/ and solve nonlinear system:
      nq2 = nq
      do 2 k = 1,np
        x2(k) = x(k)
        z2(k) = z(k)
    2   betam2(k) = betam(k)
      nwdim = nq*(2*n+3)
      do 3 i = 1,nwdim
    3   qwork2(i) = qwork(i)
      call ns01a(nm2,y,fval,ajinv,dstep,dmax,tol,maxfun,
     &  iprint,scr,kfun)
*
* copy output data from /param1/ and print results:
      c = c2
      do 4 k = 2,np
    4   x(k) = x2(k)
      call nearx(x(np),x0,z0,k0,n,x,z,betam)
      z(np) = zx(x(np),0,x0,z0,k0,n,c,x,betam,nq,qwork)
      if (iprint.ge.0) call koutp(n,c,x,z,betam,nq)
      return
      end


      function zx(xx,kxx,x0,z0,k0,n,c,x,betam,nq,qwork)
*
* computes map z(xx)
*
      implicit double precision (a-b,d-h,o-w,y)
      implicit complex*16(c,x,z)
      dimension x(1),betam(1),qwork(1)
      zx = z0 + c * xquad(x0,k0,xx,kxx,n,x,betam,nq,qwork)
      return
      end


      subroutine yxtran(n,y,x,z,betam)
*
* transforms y(k) to x(k)
*
      implicit double precision (a-b,d-h,o-w,y)
      implicit complex*16(c,x,z)
      dimension y(1),x(1),z(1),betam(1)
      nm = n-1
      dx = 1.
      sumx = dx
      do 1 k = 2,nm
        dx = dx / exp(y(k-1))
    1   sumx = sumx + dx
      dx = 2./sumx
      x(2) = -1. + dx
      do 2 k = 3,nm
        dx = dx / exp(y(k-2))
    2   x(k) = x(k-1) + dx
      asum = 0.
      do 3 k = 2,nm
    3   asum = asum + betam(k)*acos(dreal(-x(k)))
      x(n+1) = -cos( dimag(log(z(n)-z(nm))) + asum )
      return
      end


      subroutine kfun(nm2,y,fval)
*
* this is the function whose zero must be found in ksolv
*
      implicit double precision (a-b,d-h,o-w,y)
      implicit complex*16(c,x,z)
      dimension fval(nm2),y(nm2)
      common /param1/ nq,c,x(20),z(20),qwork(460),betam(20)
      nm = nm2+1
      n = nm2+2
      call yxtran(n,y,x,z,betam)
      zdenom = xquad(x(1),1,x(2),2,n,x,betam,nq,qwork)
      c = (z(2)-z(1))/zdenom
      do 1 k = 2,nm
        kp = k+1
        xint = xquad(x(k),k,x(kp),kp,n,x,betam,nq,qwork)
        fval(k-1) = abs(z(kp)-z(k)) - abs(c*xint)
    1 continue
      return
      end


      subroutine koutp(n,c,x,z,betam,nq)
*
* prints a table of k, z(k), betam(k), and x(k)
*
      implicit double precision (a-b,d-h,o-w,y)
      implicit complex*16(c,x,z)
      dimension x(1),z(1),betam(1)
      print 101, n, nq
      do 1 k = 1,n
    1   print 102, k,z(k),betam(k),x(k)
      print 103, z(n+1), x(n+1)
      print '(/''  c = ('',e15.8,'','',e15.8,'')''/)', c
      return
  101 format (/' parameters defining map:',15x,'(n =',
     &  i3,')',6x,'(nq =',i3,')'//
     &  '  k', 9x,'z(k)',11x,'betam(k)',12x,'x(k)'/
     &  ' ---',8x,'----',11x,'--------',12x,'----'/)
  102 format (i3,'    (',f6.3,',',f6.3,')',f12.5,
     &  3x,'(',f11.8,',',f11.8,')')
  103 format ('stagn. (',f6.3,',',f6.3,')',15x,'(',f11.8,',',f11.8,')')
      end


      function xquad(xa,ka,xb,kb,n,x,betam,nq,qwork)
*
* computes the complex line integral of xprod from xa to xb along a
* straight line segment.  function xquad1 is called twice.
*
      implicit double precision (a-b,d-h,o-w,y)
      implicit complex*16(c,x,z)
      dimension x(1),betam(1),qwork(1)
      xmid = (xa + xb) / 2.
      xquad = xquad1(xa,xmid,ka,n,x,betam,nq,qwork)
     &       -xquad1(xb,xmid,kb,n,x,betam,nq,qwork)
      return
      end


      function xquad1(xa,xb,ka,n,x,betam,nq,qwork)
*
* computes the complex line integral of xprod from xa to xb along a
* straight line segment by compound one-sided gauss-jacobi quadrature.
*
      implicit double precision (a-b,d-h,o-w,y)
      implicit complex*16(c,x,z)
      dimension x(1),betam(1),qwork(1)
      data resprm /1.d0/
*
* step 1: one-sided gauss-jacobi quadrature for left endpoint:
      if (abs(xa-xb).gt.0.d0) goto 1
      xquad1 = (0.,0.)
      return
    1 r = min(1.d0,distx(xa,ka,n,x)*resprm/abs(xa-xb))
      xaa = xa + r*(xb-xa)
      xquad1 = xqsum(xa,xaa,ka,n,x,betam,nq,qwork)
*
* step 2: adjoin intervals of pure gaussian quadrature if necessary:
    2 if (r.eq.1.d0) return
      r = min(1.d0,distx(xaa,0,n,x)*resprm/abs(xaa-xb))
      xbb = xaa + r*(xb-xaa)
      xquad1 = xquad1 + xqsum(xaa,xbb,0,n,x,betam,nq,qwork)
      xaa = xbb
      goto 2
      end


      function distx(xx,ks,n,x)
*
* determines the distance from xx to the nearest singularity x(k)
* other than x(ks)
*
      implicit double precision (a-b,d-h,o-w,y)
      implicit complex*16(c,x,z)
      dimension x(1)
      distx = 9999.
      do 1 k = 1,n
        if (k.eq.ks) goto 1
        distx = min(distx,abs(xx-x(k)))
    1   continue
      return
      end


      function xqsum(xa,xb,ka,n,x,betam,nq,qwork)
*
* computes the integral of xprod from xa to xb by applying a
* one-sided gauss-jacobi formula with possible singularity at xa.
*        ka=0: gauss-legendre
*   ka=1 or n: gauss-legendre after change of variables (separation pt.)
*      1<ka<n: gauss-jacobi with exponent betam(ka)
*
      implicit double precision (a-b,d-h,o-w,y)
      implicit complex*16(c,x,z)
      dimension x(1),betam(1),qwork(1)
      xs = (0.,0.)
      xh = (xb-xa) / 2.
      xc = (xa+xb) / 2.
      k = ka
      if (k.eq.0) k = n+1
      iwt1 = nq*(k-1) + 1
      iwt2 = iwt1 + nq - 1
      ioffst = nq*(n+1)
      if (k.eq.1.or.k.eq.n) goto 2
      do 1 i = iwt1,iwt2
    1   xs = xs + qwork(ioffst+i)*xprod(xc+xh*qwork(i),ka,n,x,betam)
      xqsum = xs*xh
      if (abs(xh).ne.0.d0) xqsum = xqsum*abs(xh)**betam(k)
      return
    2 continue
      do 3 i = iwt1,iwt2
        t = (qwork(i)**2 + 2.*qwork(i) - 1.) / 2.
    3   xs = xs + qwork(ioffst+i)*
     &    (qwork(i)+1.)*xprod(xc+xh*t,ka,n,x,betam)
      xqsum = xs*xh
      return
      end


      function xprod(xxx,ks,n,x,betam)
*
* evaluates the modified schwarz-christoffel integrand
*
      implicit double precision (a-b,d-h,o-w,y)
      implicit complex*16(c,x,z)
      common /dual/ sgx
      dimension x(1),betam(1)
      common /logcnt/ ncount,ncnt1
      nm = n-1
      np = n+1
      xx = xxx
      if (sgx.lt.0.d0) xx = dconjg(xx)
      xx = xx + (0.d0,1.d-12)
      xsum = sgx * log( (1.-xx*x(np)+sqrt((1.-x(np)**2)*(1.-xx**2))) /
     &  (xx-x(np)) )
      do 1 k = 2,nm
        xtmp = xx - x(k)
        if (k.eq.ks) xtmp = xtmp / abs(xtmp)
        xtmp = xtmp / (1.-xx*x(k)+sqrt((1.-x(k)**2)*(1.-xx**2)))
    1   xsum = xsum + betam(k)*log(xtmp)
      xprod = exp(xsum)
      if (sgx.lt.0.d0) xprod = dconjg(xprod)
      xprod = (xxx-x(np)) * xprod
      ncount = ncount + nm
      return
      end


      subroutine nearx(xx,x0,z0,k0,n,x,z,betam)
*
* returns information associated with the nearest prevertex x(k) to xx
*
      implicit double precision (a-b,d-h,o-w,y)
      implicit complex*16(c,x,z)
      dimension x(1),z(1),betam(1)
      dist = 9999.
      do 1 k = 1,n
        if (betam(k).lt.-.9999) goto 1
        distxk = abs(xx-x(k))
        if (distxk.ge.dist) goto 1
        dist = distxk
        k0 = k
    1   continue
      x0 = x(k0)
      z0 = z(k0)
      return
      end


      subroutine anglex(n,z,betam)
*
* given vertices z(k), computes angle parameters betam(k)
*
      implicit double precision (a-b,d-h,o-w,y)
      implicit complex*16(c,x,z)
      common /dual/ sgx
      dimension z(1),betam(1)
      pi = acos(-1.d0)
      betam(1) = 0.
      do 1 k = 2,n-1
        betam(k) = dimag(log((z(k-1)-z(k))/(z(k+1)-z(k))))/pi - 1.d0
        if (betam(k).le.-1.d0) betam(k) = betam(k) + 2.d0
    1   betam(k) = sgx*betam(k)
      betam(n) = 0.
      betam(n+1) = 0.
      return
      end


      subroutine count0
      common /logcnt/ ncount,ncnt1
      ncount = 0
      ncnt1 = 0
      print '('' ------- log counter set to zero'')'
      return
      end


      subroutine count
      common /logcnt/ ncount,ncnt1
      ncdiff = ncount - ncnt1
      print '('' ------- no. logs: since last count'',i7,
     &  '',   total'',i8)', ncdiff, ncount
      ncnt1 = ncount
      return
      end