C      $Id: x22f.fm4,v 1.4 2005/04/27 06:43:54 rlaboiss Exp $
C      Vector plot demo.
C
C      Copyright (C) 2004  Alan W. Irwin
C      Copyright (C) 2004  Andrew Ross
C
C      This file is part of PLplot.
C
C      PLplot is free software; you can redistribute it and/or modify
C      it under the terms of the GNU General Library Public License as
C      published by the Free Software Foundation; either version 2 of the
C      License, or (at your option) any later version.
C
C      PLplot is distributed in the hope that it will be useful,
C      but WITHOUT ANY WARRANTY; without even the implied warranty of
C      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C      GNU Library General Public License for more details.
C
C      You should have received a copy of the GNU Library General Public
C      License along with PLplot; if not, write to the Free Software
C      Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA

C      Does several contour plots using different coordinate mappings.
      implicit none

      integer narr, fill
      parameter (narr=6)
      real*8 arrow_x(narr),arrow_y(narr),
     &	arrow2_x(narr),arrow2_y(narr)

      data arrow_x/-0.5d0, 0.5d0, 0.3d0, 0.5d0, 0.3d0, 0.5d0/
      data arrow_y/0.d0, 0.d0, 0.2d0, 0.d0, -0.2d0, 0.d0/
      data arrow2_x/-0.5d0, 0.3d0, 0.3d0, 0.5d0, 0.3d0, 0.3d0/
      data arrow2_y/0.d0, 0.d0, 0.2d0, 0.d0, -0.2d0, 0.d0/

      integer PL_PARSE_FULL
      parameter(PL_PARSE_FULL = 1)
C      Process command-line arguments
      call plparseopts(PL_PARSE_FULL)

      call plinit


      call circulation

      fill = 0

C     Set arrow style using arrow_x and arrow_y the
C     plot using these arrows
      call plsvect(arrow_x, arrow_y, narr, fill)
      call constriction

C     Set arrow style using arrow_x and arrow_y the
C     plot using these arrows
      fill = 1
      call plsvect(arrow2_x, arrow2_y, narr, fill)
      call constriction

      call potential

      call plend
      
      end

C     vector plot of the circulation around the origin
      subroutine circulation()
      implicit none

      integer i, j, nx, ny
      parameter (nx=20, ny=20)

      real*8 u(nx, ny), v(nx, ny), xg(nx,ny), yg(nx,ny)

      real*8 dx, dy, xmin, xmax, ymin, ymax
      real*8 xx, yy, scaling

      dx = 1.0d0
      dy = 1.0d0

      xmin = -dble(nx)/2.0d0*dx
      xmax = dble(nx)/2.0d0*dx
      ymin = -dble(ny)/2.0d0*dy
      ymax = dble(ny)/2.0d0*dy

      do i=1,nx
        xx = (dble(i)-nx/2.0d0-0.5d0)*dx
        do j=1,ny
          yy = (dble(j)-ny/2.0d0-0.5d0)*dy
	  xg(i,j) = xx
	  yg(i,j) = yy
          u(i,j) = yy
          v(i,j) = -xx
        enddo
      enddo

      call plenv(xmin, xmax, ymin, ymax, 0, 0)
      call pllab('(x)', '(y)', 
     &	'#frPLplot Example 22 - circulation')
      call plcol0(2)
      scaling = 0.0d0
      call plvec2(u,v,nx,ny,scaling,xg,yg)
      call plcol0(1)

      end

C     vector plot of the flow through a constricted pipe
      subroutine constriction()
      implicit none

      real*8 PI
      parameter (PI = 3.1415926535897932384d0)

      integer i, j, nx, ny
      parameter (nx=20, ny=20)

      real*8 u(nx, ny), v(nx, ny), xg(nx,ny), yg(nx,ny)

      real*8 dx, dy, xmin, xmax, ymin, ymax
      real*8 xx, yy, Q, b, dbdx, scaling

      dx = 1.0d0
      dy = 1.0d0

      xmin = -dble(nx)/2.0d0*dx
      xmax = dble(nx)/2.0d0*dx
      ymin = -dble(ny)/2.0d0*dy
      ymax = dble(ny)/2.0d0*dy

      Q = 2.0d0
      do i=1,nx
        xx = (dble(i)-dble(nx)/2.0d0-0.5d0)*dx
        do j=1,ny
          yy = (dble(j)-dble(ny)/2.0d0-0.5d0)*dy
	  xg(i,j) = xx
	  yg(i,j) = yy
          b = ymax/4.0d0*(3.0d0-cos(PI*xx/xmax))
          if (abs(yy).lt.b) then
             dbdx = ymax/4.0d0*sin(PI*xx/xmax)*yy/b
             u(i,j) = Q*ymax/b
             v(i,j) = u(i,j)*dbdx
          else
             u(i,j) = 0.0d0
             v(i,j) = 0.0d0
          endif
        enddo
      enddo

      call plenv(xmin, xmax, ymin, ymax, 0, 0)
      call pllab('(x)', '(y)', 
     &	'#frPLplot Example 22 - constriction')
      call plcol0(2)
      scaling = -0.5d0
      call plvec2(u,v,nx,ny,scaling,xg,yg)
      call plcol0(1)

      end


      subroutine potential()
      implicit none

      real*8 PI
      parameter (PI = 3.1415926535897932384d0)

      integer i, j, nr, ntheta, nper, nlevel
      parameter (nr=20, ntheta=20, nper=100, nlevel=10)

      real*8 u(nr, ntheta), v(nr, ntheta), z(nr, ntheta)
      real*8 xg(nr,ntheta), yg(nr,ntheta)
      real*8 clevel(nlevel), px(nper), py(nper)

      real*8 xmin, xmax, ymin, ymax, zmin, zmax, rmax
      real*8 xx, yy, r, theta, scaling, dz

      real*8 eps, q1, d1, q1i, d1i, q2, d2, q2i, d2i
      real*8 div1, div1i, div2, div2i

      rmax = dble(nr)

      eps = 2.0d0

      q1 = 1.0d0;
      d1 = rmax/4.0d0;
      
      q1i = - q1*rmax/d1;
      d1i = rmax**2.0d0/d1;
      
      q2 = -1.0d0;
      d2 = rmax/4.0d0;
      
      q2i = - q2*rmax/d2;
      d2i = rmax**2.0d0/d2;
           
      do i = 1, nr
         r = 0.5 + dble(i-1)
        do j = 1, ntheta
          theta = 2.*PI/dble(ntheta-1)*(dble(j)-0.5)
          xx = r*cos(theta)
          yy = r*sin(theta) 
          xg(i,j) = xx
          yg(i,j) = yy
          div1 = sqrt((xg(i,j)-d1)**2 + (yg(i,j)-d1)**2 + eps**2)
          div1i = sqrt((xg(i,j)-d1i)**2 + (yg(i,j)-d1i)**2 + eps**2)

          div2 = sqrt((xg(i,j)-d2)**2 + (yg(i,j)+d2)**2 + eps**2)
          div2i = sqrt((xg(i,j)-d2i)**2 + (yg(i,j)+d2i)**2 + eps**2)

          z(i,j) = q1/div1 + q1i/div1i + q2/div2 + q2i/div2i
          u(i,j) = -q1*(xx-d1)/div1**3 - q1i*(xx-d1i)/div1i**3 -
     1         q2*(xx-d2)/div2**3 - q2i*(xx-d2i)/div2i**3
          v(i,j) = -q1*(yy-d1)/div1**3 - q1i*(yy-d1i)/div1i**3 -
     1         q2*(yy+d2)/div2**3 - q2i*(yy+d2i)/div2i**3
        enddo
      enddo

      call a2mnmx(xg, nr, ntheta, xmin, xmax, nr)
      call a2mnmx(yg, nr, ntheta, ymin, ymax, nr)
      call a2mnmx(z, nr, ntheta, zmin, zmax, nr)

      call plenv(xmin, xmax, ymin, ymax, 0, 0)
      call pllab('(x)', '(y)', 
     &	'#frPLplot Example 22 - potential gradient vector plot')

C     plot contours of the potential
      dz = abs(zmax - zmin)/dble (nlevel)
      do i = 1, nlevel
         clevel(i) = zmin + (i-0.5d0)*dz
      enddo
      call plcol0(3)
      call pllsty(2)
      call plcon2(z,nr,ntheta,1,nr,1,ntheta,clevel,nlevel,xg,yg)
      call pllsty(1)
      call plcol0(1)

      call plcol0(2)
      scaling = 25.0d0
      call plvec2(u,v,nr,ntheta,scaling,xg,yg)
      call plcol0(1)

      do i=1,nper
         theta = 2.0d0*PI/dble(nper-1)*dble(i)
         px(i) = rmax*cos(theta)
         py(i) = rmax*sin(theta)
      enddo

      call plline(nper,px,py)

      end

C----------------------------------------------------------------------------
C      Subroutine a2mnmx
C      Minimum and the maximum elements of a 2-d array.

      subroutine a2mnmx(f, nx, ny, fmin, fmax, xdim)
      implicit none

      integer   i, j, nx, ny, xdim
      real*8    f(xdim, ny), fmin, fmax

      fmax = f(1, 1)
      fmin = fmax
      do j = 1, ny
        do  i = 1, nx
          fmax = max(fmax, f(i, j))
          fmin = min(fmin, f(i, j))
        enddo
      enddo
      end
