C      $Id: x16af.fm4,v 1.6 2005/04/27 06:43:54 rlaboiss Exp $
C      Demonstration of plshade plotting
C      Reduce colors in cmap 0 so that cmap 1 is useful on a 16-color display
C
C      Copyright (C) 2004  Alan W. Irwin
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

      implicit none

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

      call plscmap0n(3)

C      Initialize plplot

      call plinit()

C      Rectangular coordinate plot

      call rect()

C      Polar coordinate plot

      call polar()

      call plend
      end

C      Plot function using the identity transform

      subroutine rect()

      implicit none
      integer   xdim, ydim, NX, NY, NCONTR
C      xdim and ydim are the static dimensions of the 2D arrays while
C      NX and NY are the defined area.
      parameter (xdim = 99, NX = 35, ydim = 100, NY = 46, NCONTR = 14)

      real*8    z(xdim, ydim), w(xdim, ydim), clevel(NCONTR)
      real*8    xmin, xmax, ymin, ymax, zmin, zmax, x, y
      real*8    shade_min, shade_max, sh_color
      integer   i, j, sh_cmap, sh_width
      integer   min_color, min_width, max_color, max_width

      xmin = -1.d0
      ymin = -1.d0
      xmax =  1.d0
      ymax =  1.d0

C      Set up for plshade call

      sh_cmap   = 1
      min_color = 1
      min_width = 0
      max_color = 0
      max_width = 0

C      Set up data arrays

      do i = 1, NX
        x = (i - 1 - (NX/2)) / dble (NX/2)
        do j = 1, NY
          y = (j - 1 - (NY/2)) / dble (NY/2) - 1.0d0
          z(i,j) = x*x - y*y + (x - y) / (x*x + y*y + 0.1d0)
          w(i,j) = 2*x*y
        enddo
      enddo

      call a2mnmx(z, NX, NY, zmin, zmax, xdim)
      do  i = 1, NCONTR
        clevel(i) = zmin + (zmax - zmin) * (i + 0.5d0) /
     &    dble (NCONTR)
      enddo

C      Plot using identity transform

      call pladv(0)
      call plvpor(0.1d0, 0.9d0, 0.1d0, 0.9d0)
      call plwind(-1.0d0, 1.0d0, -1.0d0, 1.0d0)

      do i = 1, NCONTR
        shade_min = zmin + (zmax - zmin) * dble (i - 1) /
     &    dble (NCONTR)
        shade_max = zmin + (zmax - zmin) * dble (i)     /
     &    dble (NCONTR)
        sh_color = dble (i - 1) / dble (NCONTR - 1)
        sh_width = 2
        call plpsty(0)
        call plshade0(z, NX, NY, ' ',
     &    dble ( -1.d0), 1.0d0, -1.0d0, 1.0d0,
     &    shade_min, shade_max,
     &    sh_cmap, sh_color, sh_width,
     &    min_color, min_width, max_color, max_width, xdim)
      enddo

      call plcol(1)
      call plbox('bcnst', 0.0d0, 0, 'bcnstv', 0.0d0, 0)
      call plcol(2)
      call pllab('distance', 'altitude', 'Bogon flux')

      end

C      Routine for demonstrating use of transformation arrays in contour plots.

      subroutine polar()

      implicit none
      integer   xdim, ydim, NX, NY, NCONTR, NBDRY
      real*8    TWOPI
C      xdim and ydim are the static dimensions of the 2D arrays while
C      NX and NY are the defined area.
      parameter (xdim = 99, NX = 40, ydim = 100, NY = 64)
      parameter (NCONTR = 14, NBDRY=200)
      parameter (TWOPI=6.2831853071795864768d0)

      real*8    z(xdim, ydim), ztmp(xdim, ydim+1)
      real*8    xg(xdim, ydim+1), yg(xdim, ydim+1),
     &   xtm(NBDRY), ytm(NBDRY)
      real*8    clevel(NCONTR)
      real*8    xmin, xmax, ymin, ymax, zmin, zmax
      real*8    xpmin, xpmax, ypmin, ypmax
      real*8    r, theta, rmax, x0, y0
      real*8    eps, q1, d1, q1i, d1i, q2, d2, q2i, d2i
      real*8    div1, div1i, div2, div2i

      real*8    shade_min, shade_max, sh_color
      real*8    xtick, ytick
      integer   nxsub, nysub
      integer   ncolbox, ncollab
      integer   i, j, kx, lx, ky, ly
      integer   sh_cmap, sh_width
      integer   min_color, min_width, max_color, max_width
      character*8 xopt, yopt

C      Set up for plshade call

      sh_cmap = 1
      min_color = 1
      min_width = 0
      max_color = 0
      max_width = 0

      kx = 1
      lx = NX
      ky = 1
      ly = NY

C      Set up r-theta grids
C      Tack on extra cell in theta to handle periodicity.

      do i = 1, NX
        r = i - 0.5d0
        do j = 1, NY
          theta = TWOPI/dble(NY) * (j-0.5d0)
          xg(i,j) = r * cos(theta)
          yg(i,j) = r * sin(theta)
        enddo
        xg(i, NY+1) = xg(i, 1)
        yg(i, NY+1) = yg(i, 1)
      enddo
      call a2mnmx(xg, NX, NY, xmin, xmax, xdim)
      call a2mnmx(yg, NX, NY, ymin, ymax, xdim)

      rmax = r
      x0 = (xmin + xmax)/2.d0
      y0 = (ymin + ymax)/2.d0

C      Potential inside a conducting cylinder (or sphere) by method of images.
C      Charge 1 is placed at (d1, d1), with image charge at (d2, d2).
C      Charge 2 is placed at (d1, -d1), with image charge at (d2, -d2).
C      Also put in smoothing term at small distances.

      eps = 2.d0

      q1 = 1.d0
      d1 = r/4.d0

      q1i = - q1*r/d1
      d1i = r**2/d1

      q2 = -1.d0
      d2 = r/4.d0

      q2i = - q2*r/d2
      d2i = r**2/d2

      do i = 1, NX
        do j = 1, NY
          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
        enddo
      enddo

C      Tack on extra cell in theta to handle periodicity.

      do i = 1, NX
        do j = 1, NY
          ztmp(i,j) = z(i,j)
        enddo
        ztmp(i, NY+1) = z(i, 1)
      enddo
      call a2mnmx(z, NX, NY, zmin, zmax, xdim)

C      Set up contour levels.

      do i = 1, NCONTR
        clevel(i) = zmin + (i-0.5d0)*abs(zmax - zmin)/dble(NCONTR)
      enddo

C      Advance graphics frame and get ready to plot.

      ncolbox = 1
      ncollab = 2

      call pladv(0)
      call plcol(ncolbox)

C      Scale window to user coordinates.
C      Make a bit larger so the boundary does not get clipped.

      eps = 0.05d0
      xpmin = xmin - abs(xmin)*eps
      xpmax = xmax + abs(xmax)*eps
      ypmin = ymin - abs(ymin)*eps
      ypmax = ymax + abs(ymax)*eps

      call plvpas(0.1d0, 0.9d0, 0.1d0, 0.9d0, 1.0d0)
      call plwind(xpmin, xpmax, ypmin, ypmax)

      xopt = ' '
      yopt = ' '
      xtick = 0.d0
      nxsub = 0
      ytick = 0.d0
      nysub = 0

      call plbox(xopt, xtick, nxsub, yopt, ytick, nysub)

C      Call plotter once for z < 0 (dashed), once for z > 0 (solid lines).

      do i = 1, NCONTR
        shade_min = zmin + (zmax - zmin) * dble (i - 1) /
     &    dble (NCONTR)
        shade_max = zmin + (zmax - zmin) * dble (i)     /
     &    dble (NCONTR)
        sh_color = dble (i - 1) / dble (NCONTR - 1)
        sh_width = 2
        call plpsty(0)

        call plshade2(z, NX, NY, ' ',
     &    -1.0d0, 1.0d0, -1.0d0, 1.0d0,
     &    shade_min, shade_max,
     &    sh_cmap, sh_color, sh_width,
     &    min_color, min_width, max_color, max_width, xg, yg, xdim)
      enddo

C      Draw boundary.

      do i = 1, NBDRY
        theta = (TWOPI)/(NBDRY-1) * dble(i-1)
        xtm(i) = x0 + rmax * cos(theta)
        ytm(i) = y0 + rmax * sin(theta)
      enddo
      call plcol(ncolbox)
      call plline(NBDRY, xtm, ytm)

      call plcol(ncollab)
      call pllab(' ', ' ',
     &'Shielded potential of charges in a conducting sphere')

      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
