program test1
      implicit complex(c,w,z)
      dimension z(20),w(20),betam(20),qwork(344)
      zero = (0.,0.)
      zi = (0.,1.)
c
c set up problem:
      n = 4
      wc = cmplx(0.,sqrt(2.))
      w(1) = zi
      w(2) = zero
      w(3) = (1.e20,1.e20)
      w(4) = zero
      betam(1) = 1.
      betam(2) = -.5
      betam(3) = -2.
      betam(4) = -.5
c
c compute nodes and weights for parameter problem:
      nptsq = 5
      call qinit(n,betam,nptsq,qwork)
c
c solve parameter problem:
c   (initial guess must be given to avoid accidental exact solution)
      iprint = 0
      iguess = 1
      do 1 k = 1,4
    1 z(k) = exp(cmplx(0.,k-4.))
      tol = 1.e-6
      call scsolv(iprint,iguess,tol,errest,n,c,z,
     &  wc,w,betam,nptsq,qwork)
c
c compare wsc(z) to exact values for various z:
      do 10 i = 1,4
        zz = (.3,0.) * cmplx(i-2.,.2*i+.5)
        ww = wsc(zz,0,zero,wc,0,n,c,z,betam,nptsq,qwork)
        ztmp = -zi * (zz-zi) / (zz+zi)
        wwex = zi * sqrt(-ztmp**2 + (1.,0.))
        err = abs(ww-wwex)
        write (6,201) zz,ww,wwex,err
   10   continue
      write (6,200)
c
c compare zsc(w) to exact values for various w:
      do 20 i = 1,6
        ww = cmplx(i-2.,sqrt(i+1.))
        ier = 0
        zz = zsc(ww,0,zero,zero,wc,0,tol,ier,
     &    n,c,z,wc,w,betam,nptsq,qwork)
        wtmp = zi * sqrt((-1.,0.)-ww**2)
        zzex = -zi * (wtmp-zi) / (wtmp+zi)
        err = abs(zz-zzex)
        write (6,202) ww,zz,zzex,err
   20   continue
c
      stop
  200 format (1x)
  201 format (' z,w,wex,err: ',3('(',f6.3,',',f6.3,') '),d11.4)
  202 format (' w,z,zex,err: ',3('(',f6.3,',',f6.3,') '),d11.4)
      end
      program test2
      implicit complex (c,w,z)
      dimension w(10),ibrk(4),qwork(300)
      n = 6
      w(1) = (0.,0.)
      w(2) = (2.,0.)
      w(3) = (2.,1.)
      w(4) = (1.,1.)
      w(5) = (1.,2.)
      w(6) = (0.,2.)
      wc = (.5,.5)
      ibrk(1) = 1
      ibrk(2) = 2
      ibrk(3) = 5
      ibrk(4) = 6
c
c main loop: different accuracy specifications:
      do 10 ndig = 1,5,2
        write (6,202) ndig
        r = resist(n,w,wc,ibrk,ndig,errest,qwork)
        write (6,201) r,errest
   10   continue
      stop
c
  201 format ('   r,errest:',2d23.15)
  202 format (/' ndig =',i3,':')
      end
      function resist(n,w,wc,ibrk,ndig,errest,qwork)
c
c this function returns the resistance of a polygonally
c shaped resistor.  computations are based on the schwarz-
c christoffel transformation.  we normalize by assuming that
c a square has resistance 1.
c
c input parameters:
c
c   n      number of vertices of the polygon.  n must satisfy
c          4 .le. n .le. 20.
c
c   w      complex array of dimension at least n containing
c          the positions of the vertices, viewed as complex numbers.
c          the vertices must be listed in counterclockwise order
c          around the polygon.  it is a good idea to keep the
c          w(k) roughly on the scale of unity.
c
c   wc     a point in the interior of the polygon.  try to
c          pick wc as central as possible in the sense that
c          as few parts of the polygon as possible are shielded
c          from it by reentrant edges.
c
c   ibrk   array of dimension at least 4 containing indices of
c          the vertices which define the breaks
c          between constant-voltage and insulated portions
c          of the boundary.  the program will assume that
c          the voltage is applied between side w(ibrk(1))-w(ibrk(2))
c          and side w(ibrk(3))-w(ibrk(4)), with the other two
c          sides insulated.  the break vertices must be numbered
c          in counterclockwise order; thus the integers ibrk(i)
c          must increase with i, except that they may wrap
c          once around n.
c
c   ndig   input integer giving the desired number of digits
c          of accuracy in the result.  ndig must be at least 1.
c          it should be no more than a couple of digits less
c          than full-word precision.
c
c output parameter:
c
c   errest rough but conservative estimate of the size of
c          the error in the value returned.
c
c work space parameter:
c
c   qwork  real work array.  dimension at least (ndig+1)*(2n+3),
c          but no more than 460.
c
c
c some advice:
c
c   the program is not infallible.  if it yields strange
c   messages about not converging, it's probably best to
c   try a simpler geometry.
c
c   the amount of time required is roughly proportional to ndig
c   and also roughly proportional to n**3.
c   thus problems with many corners can
c   take quite a while -- several minutes of cpu time on
c   the ibm 370-168 for a problem with n = 20.
c
c lloyd n. trefethen
c courant institute of mathematical sciences
c new york university
c november 1979, revised july 1983
c
      implicit complex(c,w,z)
c     map from disk to resistor:
      dimension z(20),w(1),ibrk(1),betam(20),qwork(1)
c     map from disk to rectangle with equal resistance:
      dimension z2(4),w2(4),betam2(4)
c
      zero = (0.,0.)
      pi = acos(-1.)
c
c compute angles and check input parameters:
      call angles(n,w,betam)
      if (ndig.lt.1) write (6,101)
      if (ndig.lt.1) stop
      if (n.lt.4) write (6,102)
      if (n.lt.4) stop
      do 1 k = 1,4
        if (ibrk(k).gt.n .or. ibrk(k).lt.1) goto 2
    1   continue
      goto 3
    2 write (6,103)
      stop
    3 continue
c
c set up first map:
      nptsq = ndig + 1
      tol = max(1.e-5, 10.**(-nptsq))
      call qinit(n,betam,nptsq,qwork)
      call scsolv(-2,0,tol,eest,n,c,z,wc,w,betam,nptsq,qwork)
c
c set up map to rectangle:
c (qwork is overwritten to save space)
      n2 = 4
      c2 = (1.,0.)
      do 9 k = 1,4
        betam2(k) = -.5
    9   z2(k) = z(ibrk(k))
      nptsq2 = ndig + 1
      call qinit(n2,betam2,nptsq2,qwork)
      do 12 k = 1,4
   12   w2(k)=wsc(z2(k),k,zero,zero,0,n2,c2,z2,betam2,nptsq2,qwork)
c
c compute length, width, resistance, and error estimate:
      xlen1 = abs(w2(2)-w2(3))
      xlen2 = abs(w2(4)-w2(1))
      errl = max(abs(xlen1-xlen2),2.*eest) / xlen1
      xwid1 = abs(w2(2)-w2(1))
      xwid2 = abs(w2(4)-w2(3))
      errw = max(abs(xwid1-xwid2),2.*eest) / xwid1
      resist = (xlen1+xlen2) / (xwid1+xwid2)
      errest = resist * (errl + errw)
c
      return
c
  101 format (' *** error in resist *** ndig should be at',
     &  ' least 1.'/)
  102 format (' *** error in resist *** n must be no',
     &  ' greater than 20'/' and no smaller than 4')
  103 format (' *** error in resist *** each ibrk(i) must',
     &  ' be in the',/,' range from 1 to n')
      end
      program test3
c
c computes the conformal map from a polygon (1) onto a rectangle (2)
c and plots a corresponding grid of size nx by ny.
c the corners of the rectangle are (0,i), (0,0), (h,0), (h,i).
c
c nick trefethen, icase, july 1983
c
      implicit complex(c,w,z)
      dimension wgrid(0:41,0:41)
      dimension z1(12),w1(12),betam1(12),qwork1(270),ibrk(4)
      dimension z2(4),w2(4),betam2(4),qwork2(110)
      data zero /(0.,0.)/
c
c input data
      print 90
   90 format (' vertices? (terminate with re w(k) > 100 )')
      k = 0
   91 k = k+1
      read *,x,y
      w1(k) = cmplx(x,y)
      if (x.lt.100.) goto 91
      n1 = k-1
      print 92
   92 format (' image of zero in the polygon?')
      read *,x,y
      wc1 = cmplx(x,y)
      print 93
   93 format (' four distinguished vertices?')
      read *,(ibrk(k),k=1,4)
      nq1 = 3
      nq2 = nq1
      tol = 10.**(-(nq1+1))
      print 94
   94 format (' no. of streamlines, equipotential lines?')
      read *,ny,nx
      nxp = nx + 1
      nyp = ny + 1
c
c compute parameters for map from disk to polygon
      call angles(n1,w1,betam1)
      call qinit(n1,betam1,nq1,qwork1)
      call scsolv(0,0,tol,eest,n1,c1,z1,wc1,w1,betam1,nq1,qwork1)
c
c compute parameters for map from disk to rectangle
      n2 = 4
      c2 = (1.,0.)
      do 2 k = 1,4
        betam2(k) = -.5
    2   z2(k) = z1(ibrk(k))
      call qinit(n2,betam2,nq2,qwork2)
      do 3 k = 1,4
    3   w2(k) = wsc(z2(k),k,zero,zero,0,n2,c2,z2,betam2,nq2,qwork2)
      c2 = (0.,1.) / (w2(1)-w2(2))
      wc2 = -c2*w2(2)
      do 4 k = 1,4
    4   w2(k) = c2*w2(k) + wc2
      write (6,102) (w2(k),k=1,4)
  102 format (' vertices of image rectangle: ',
     &  4(/'   (',e13.6,',',e13.6,')')/)
      h = abs(w2(3)-w2(2))
c
c compute grid points
      do 12 j = 0,nyp
        i1 = 0
        i2 = nxp
        if (j.eq.0.or.j.eq.nyp) i1 = 1
        if (j.eq.0.or.j.eq.nyp) i2 = nx
        do 11 i = i1,i2
          ww2 = cmplx((i*h)/nxp,float(j)/nyp)
          call nearw(ww2,zn2,wn2,kn2,n2,z2,wc2,w2,betam2)
          ier = 0
          iguess = 1
          if (i.eq.i1) iguess = 0
          zz = zsc(ww2,iguess,zz,zn2,wn2,kn2,1.e-3,ier,n2,c2,
     &      z2,wc2,w2,betam2,nq2,qwork2)
          call nearz(zz,zn1,wn1,kn1,n1,z1,wc1,w1,betam1)
          wgrid(i,j) = wsc(zz,0,zn1,wn1,kn1,n1,c1,z1,betam1,nq1,qwork1)
   11     continue
   12   write (6,105) j,nyp
  105   format (' finished row',i3,'/',i2,' of grid points')
c
c draw plot
      do 5 k = 1,n1
    5   write (10,103) w1(k)
      write (10,104) w1(1)
  103 format (2f10.5)
  104 format (2f10.5,'" "')
      do 6 j = 1,ny
        write (10,103) (wgrid(i,j),i=0,nx)
    6   write (10,104) wgrid(nxp,j)
      do 7 i = 1,nx
        write (10,103) (wgrid(i,j),j=0,ny)
    7   write (10,104) wgrid(i,nyp)
c
      stop
      end

.