C<PRE>
C     This program displays an animated color contour plot of the function
C
C	  f(x,y) = cos(a*x**2 + y**2)/sqrt(a*x**2 + y**2 + 1)
C
C     for various values of a, using the Pittsburgh Supercomputing
C     Center's drawcgm library to create a cgm graphics metafile.
C     The cgm file can then be displayed to an X windows display, a
C     Macintosh, or a PC using the Pittsburgh Supercomputing Center's 
C     gplotm, MacGplot, or PCPlot programs, respectively.  These are 
C     available by anonymous ftp from ftp.psc.edu.
C
C     ***************************   variable key  ****************************
C     length,paramlen
C		the length of each side of the square membrane.
C     n,ndim	the number of interior grid points in each of the x and y
C  		directions.
C     dx	the grid spacing.
C     u		the array of values of f.
C     color	an integer representing a color from the color table.  0 is
C		the background color, usually white, and 1 is the default
C		label color, usually black.  The other colors range from
C		2 to 255 and depend on which color table is chosen.
C     ***********************************************************************
C
      parameter (ndim=100, paramlen=0.5, tmax=10.0)
      real length,dx,u(0:ndim+1,0:ndim+1),charsize
      integer n,i,j,k,frame,color
      character*45 title
C
      n = ndim
      length = paramlen
      dx = length/float(n+1)
C
C     ***********************************************************************
C     Specify the CGM graphics output file.
C     ***********************************************************************
C
      call outfil("raindrops.cgm")
C
C     ***********************************************************************
C     Set background color to blue and text color to white for the title page.
C     ***********************************************************************
C
      redvalue = 0.0
      greenvalue = 0.0
      bluevalue = 1.0
      call setclr(0,redvalue,greenvalue,bluevalue)
C
      redvalue = 1.0
      greenvalue = 1.0
      bluevalue = 1.0
      call setclr(1,redvalue,greenvalue,bluevalue)
C
C     ***********************************************************************
C     Initialize the graphics software
C     ***********************************************************************
C
      call grfini()
C
C     ***********************************************************************
C     Create a title frame.
C     ***********************************************************************
C
      title = 'Raindrops'
      xmin = 0.4
      ymin = 0.7
      charsize = 0.05
      color = 1
      call label(xmin,ymin,title,color,charsize)
C
      title = 'H. Edward Donley'
      xmin = 0.1
      ymin = 0.15
      charsize = 0.025
      color = 1
      call label(xmin,ymin,title,color,charsize)
C
      title = 'Indiana University of Pennsylvania'
      xmin = 0.1
      ymin = 0.1
      charsize = 0.025
      color = 1
      call label(xmin,ymin,title,color,charsize)
C
C     ***********************************************************************
C     Establish a color table
C     ***********************************************************************
C
C        ********************************************************************
C        Uncomment the next statement if you want to use a predefined color
C        table.  See /usr/local/doc/drawcgm.doc for a list of these tables.
C        ********************************************************************
C
C      call setctb(5)
C
C        ********************************************************************
C        Uncomment the next statement if you want to read a color table from
C        your own file.  You can create your own color table using the
C        program mkspectrum.f.
C        ********************************************************************
C
      call getctb(0,255,'spectrum.ctb',ierror)
C
C     ***********************************************************************
C     Calculation of 2d array of function values
C     ***********************************************************************
C
      do 30 frame = 1, 20
         a = float(frame) * 0.1
         do 20 i = 0, n+1
	    x2 = a * (20.0*float(i)/float(n+1) - 10.0)**2
	    do 10 j = 0, n+1 
	      x2y2 = x2 + (20.0*float(j)/float(n+1)-10.0)**2
	      u(i,j) = cos(x2y2)/sqrt(x2y2 + 1.0)
   10	    continue
   20    continue
C
C        ********************************************************************
C        Graph a color contour plot of the function
C        ********************************************************************
C
         call graph(u,length,n,dx)
C
   30 continue
C
C     ***********************************************************************
C     Close the graphics software
C     ***********************************************************************
C
      call grfcls
C
      stop
      end
C
C
      subroutine graph(u,length,n,dx)
C     This subroutine generates and displays one frame of the animation using
C     drawcgm.  The function is displayed as a color contour plot.
C
C     ****************************  variable key  ****************************
C     length	the length of each side of the square membrane.
C     n		the number of interior grid points in each of the x and y
C  		directions.
C     dx	the grid spacing.
C     u		the displacement from the xy plane at time k*dt.
C     ubig	the displacements interpolated onto a finer grid.
C     iubig	the conversion of ubig to colors--represented by integers from
C		2 to 255.
C     ibigscale	the size of the finer grid, used to give smoother transitions
C		between colors.
C		indicates whether or not the animation frame is the first one.
C     ***********************************************************************
C
      parameter (ibigscale = 200)
      real u(0:n+1,0:n+1),length,dx,ubig(ibigscale,ibigscale)
      integer n,k,iubig(ibigscale,ibigscale)
      character title*49
C
C
C     ***********************************************************************
C     Start a new frame in the animation
C     ***********************************************************************
C
      call newfrm()
C
C     ***********************************************************************
C     Use piecewise linear interpolation to enlarge the data set.  This gives
C     a smoother transition between colored regions.
C     ***********************************************************************
C
      nsmallx = n + 2
      nsmally = n + 2
      nbigx = ibigscale
      nbigy = ibigscale
      call strlin(u,nsmallx,nsmally,ubig,nbigx,nbigy)
C
C     ***********************************************************************
C     Scale the real values in u into a range of integers, 2 to 255, that
C     span the color table.
C     ***********************************************************************
C
      umax = ubig(1,1)
      umin = ubig(1,1)
      do 20 i = 1,ibigscale
	 do 10 j = 1,ibigscale
	    uij = ubig(i,j)
	    if(uij .gt. umax)then
	       umax = uij
	    elseif(uij .lt. umin)then
	       umin = uij
	    endif
   10	 continue
   20 continue
      if(umin .eq. umax)then
	 umin = umin - 1.0
	 umax = umax + 1.0
      endif
      call rtoint(ubig,iubig,nbigx,nbigy,umin,umax,2,255)
C
C     ***********************************************************************
C     Generate the color image from the integer array.  The image is positioned
C     on the display surface using drawcgm's coordinate system.  The display
C     surface is represented by x and y coordinates ranging from 0.0 to 1.0.
C     ***********************************************************************
C
      xmin = 0.1
      ymin = 0.25
      xmax = 0.6
      ymax = 0.75
      call drawit(iubig,nbigx,nbigy,xmin,ymin,xmax,ymax)
C
C     ***********************************************************************
C     Generate a vertical color bar.
C     ***********************************************************************
C
      xmin = 0.7
      ymin = 0.25
      xmax = 0.8
      ymax = 0.75
      charsize = 0.025
      call getctb(0,255,'ctb6.dat',ierror)
      call vrtcbr(xmin,ymin,xmax,ymax,2,255,'below','above',1,charsize)
C
C     ***********************************************************************
C     Draw a caption below the graph.
C     ***********************************************************************
C
      title = 'f(x,y) = cos(a*x**2 + y**2)/(a*x**2 + y**2 + 1.0)'
      charsize = 0.025
      xmin = 0.1
      ymin = 0.2
      call label(xmin,ymin,title,1,charsize)
C
      return
      end
C</PRE>