C      $Id: x14f.fm4,v 1.5 2006/05/27 22:02:53 hbabcock Exp $
C      Demo of multiple stream/window capability
C
C      Copyright (C) 2004  Arjen Markus
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

C Plots several simple functions from other example programs.
C
C This version sends the output of the first 4 plots (one page) to two
C independent streams.

      implicit none

      integer digmax
      integer xleng0 , yleng0 , xoff0 , yoff0
      integer xleng1 , xoff1 , yoff1

      character*80 driver
      character*15 geometry_master
      character*15 geometry_slave

      real*8  x(101), y(101)
      real*8  xs(6), ys(6)
      real*8  xscale, yscale, xoff, yoff
      common /plotdat/ x, y, xs, ys, xscale, yscale, xoff, yoff
C      some fortran compilers demand typing of intrinsic lnblnk, and
C      although this is not demanded on g77 it also works there.
      integer lnblnk
      integer PL_PARSE_FULL
      parameter(PL_PARSE_FULL = 1)

      xleng0 = 400
      yleng0 = 300
      xoff0  = 200
      yoff0  = 200
      xleng1 = 400
      xleng1 = 300
      xoff1  = 500
      yoff1  = 500

      geometry_master = '500x410+100+200'
      geometry_slave = '500x410+650+200'

C      Process command-line arguments
      call plparseopts(PL_PARSE_FULL)

      call plgdev(driver)

      write(*,*) 'Demo of multiple output streams via the ',
     &  driver(:lnblnk(driver)), ' driver.'
      write(*,*) 'Running with the second stream as slave ',
     &  'to the first.'
      write(*,*)

C      Set up first stream

      call plsetopt( 'geometry', geometry_master)

      call plsdev(driver)
      call plssub(2, 2)
      call plinit()

C      Start next stream

      call plsstrm(1)

C      Turn off pause to make this a slave (must follow master)

      call plsetopt( 'geometry', geometry_slave)
      call plspause(0)
      call plsdev(driver)
      call plinit()

C      Set up the data & plot
C      Original case

      call plsstrm(0)

      xscale = 6.d0
      yscale = 1.d0
      xoff = 0.d0
      yoff = 0.d0
      call plot1()

C      Set up the data & plot

      xscale = 1.d0
      yscale = 1.d+6
      call plot1()

C      Set up the data & plot

      xscale = 1.d0
      yscale = 1.d-6
      digmax = 2
      call plsyax(digmax, 0)
      call plot1()

C      Set up the data & plot

      xscale = 1.d0
      yscale = 0.0014d0
      yoff = 0.0185d0
      digmax = 5
      call plsyax(digmax, 0)
      call plot1()

C      To slave
C      The pleop() ensures the eop indicator gets lit.

      call plsstrm(1)
      call plot4()
      call pleop()

C      Back to master

      call plsstrm(0)
      call plot2()
      call plot3()

C      To slave

      call plsstrm(1)
      call plot5()
      call pleop()

C      Back to master to wait for user to advance

      call plsstrm(0)
      call pleop()

C      Call plend to finish off.

      call plend()
      end

C======================================================================

      subroutine plot1()
      implicit none

      real*8 x(101), y(101)
      real*8 xs(6), ys(6)
      real*8  xscale, yscale, xoff, yoff,
     &  xmin, xmax, ymin, ymax
      integer i
      common /plotdat/ x, y, xs, ys, xscale, yscale, xoff, yoff

      do i = 1, 60
        x(i) = xoff + xscale * dble(i)/60.0d0
        y(i) = yoff + yscale * x(i)**2
      enddo

      xmin = x(1)
      xmax = x(60)
      ymin = y(1)
      ymax = y(60)

      do i = 1, 6
        xs(i) = x((i-1)*10+4)
        ys(i) = y((i-1)*10+4)
      enddo

C      Set up the viewport and window using PLENV. The range in X is
C      0.0 to 6.0, and the range in Y is 0.0 to 30.0. The axes are
C      scaled separately (just = 0), and we just draw a labelled
C      box (axis = 0).

      call plcol0(1)
      call plenv( xmin, xmax, ymin, ymax, 0, 0 )
      call plcol0(6)
      call pllab( '(x)', '(y)', '#frPLplot Example 1 - y=x#u2' )

C      Plot the data points

      call plcol0(9)
      call plpoin( 6, xs, ys, 9 )

C      Draw the line through the data

      call plcol0(4)
      call plline( 60, x, y )
      call plflush
      end

C======================================================================

      subroutine plot2()
      implicit none
      real*8  x(101), y(101)
      real*8  xs(6), ys(6)
      real*8  xscale, yscale, xoff, yoff
      integer i
      common /plotdat/ x, y, xs, ys, xscale, yscale, xoff, yoff

C======================================================================
C
C      Set up the viewport and window using PLENV. The range in X is
C      -2.0 to 10.0, and the range in Y is -0.4 to 2.0. The axes are
C      scaled separately (just = 0), and we draw a box with axes
C      (axis = 1).

      call plcol0(1)
      call plenv(-2.0d0, 10.0d0, -0.4d0, 1.2d0, 0, 1 )
      call plcol0(2)
      call pllab( '(x)', 'sin(x)/x',
     &            '#frPLplot Example 1 - Sinc Function' )

C      Fill up the arrays

      do i = 1, 100
        x(i) = (i-20.0d0)/6.0d0
        y(i) = 1.0d0
        if (x(i) .ne. 0.0d0) y(i) = sin(x(i)) / x(i)
      enddo
C      Draw the line

      call plcol0(3)
      call plline( 100, x, y )
      call plflush
      end

C======================================================================

      subroutine plot3()
C
C      For the final graph we wish to override the default tick intervals,
C      and so do not use PLENV
      implicit none

      real*8 PI
      parameter (PI = 3.1415926535897932384d0)
      real*8  x(101), y(101)
      real*8  xs(6), ys(6)
      real*8  xscale, yscale, xoff, yoff
      integer i
      common /plotdat/ x, y, xs, ys, xscale, yscale, xoff, yoff
      call pladv(0)

C      Use standard viewport, and define X range from 0 to 360 degrees,
C      Y range from -1.2 to 1.2.

      call plvsta()
      call plwind( 0.0d0, 360.0d0, -1.2d0, 1.2d0 )

C      Draw a box with ticks spaced 60 degrees apart in X, and 0.2 in Y.

      call plcol0(1)
      call plbox( 'bcnst', 60.0d0, 2, 'bcnstv', 0.2d0, 2 )

C      Superimpose a dashed line grid, with 1.5 mm marks and spaces. With
C      only a single mark and space element, we do not need arrays

      call plstyl( 1, 1500, 1500 )
      call plcol0(2)
      call plbox( 'g', 30.0d0, 0, 'g', 0.2d0, 0 )
      call plstyl( 0, 0, 0 )

      call plcol0(3)
      call pllab( 'Angle (degrees)', 'sine',
     &            '#frPLplot Example 1 - Sine function' )

      do i = 1, 101
        x(i) = 3.6d0 * (i-1)
        y(i) = sin( x(i) * PI/180.0d0 )
      enddo

      call plcol0(4)
      call plline( 101, x, y )
      call plflush
      end

C======================================================================

      subroutine plot4()
      implicit none

      character*3 text
      real*8 x0(0:360), y0(0:360)
      real*8 x(0:360), y(0:360), dtr, theta, dx, dy, r
      integer i, j, nsp
      real*8 PI
      parameter (PI = 3.1415926535897932384d0)

      dtr = PI/180.0d0
      do i=0,360
        x0(i) = cos(dtr * dble (i))
        y0(i) = sin(dtr * dble (i))
      enddo

C      Set up viewport and window, but do not draw box

      call plenv(-1.3d0, 1.3d0, -1.3d0, 1.3d0, 1, -2)
      do i = 1,10
        do j = 0,360
          x(j) = 0.1d0*i*x0(j)
          y(j) = 0.1d0*i*y0(j)
        enddo

C        Draw circles for polar grid

        call plline(361,x,y)
      enddo
      call plcol0(2)
      do i = 0,11
        theta = 30.0d0*i
        dx = cos(dtr*theta)
        dy = sin(dtr*theta)

C        Draw radial spokes for polar grid

        call pljoin(0.0d0, 0.0d0, dx, dy)
        write (text,'(i3)') nint(theta)

C        Write labels for angle

        text = text(nsp(text):)
C        Slightly off zero to avoid floating point logic flips at
C        90 and 270 deg.
        if (dx.ge.-0.00001d0) then
          call plptex(dx, dy, dx, dy, -0.15d0, text)
        else
          call plptex(dx, dy, -dx, -dy, 1.15d0, text)
        end if
      enddo
C      Draw the graph

      do i=0,360
        r = sin(dtr*dble (5*i))
        x(i) = x0(i) * r
        y(i) = y0(i) * r
      enddo
      call plcol0(3)
      call plline(361,x,y)

      call plcol0(4)
      call plmtex('t', 2.0d0, 0.5d0, 0.5d0,
     &  '#frPLplot Example 3 - r(#gh)=sin 5#gh')

C      Flush the plot at end

      call plflush
      end

C======================================================================

      integer function nsp(text)
C      ==================

C      Find first non-space character
      implicit none

      character*(*) text
      integer l, len

      l = len(text)
      nsp = 1
      do while(text(nsp:nsp).eq.' ' .and. nsp.lt.l)
        nsp = nsp+1
      enddo
      end

C======================================================================

      subroutine plot5()

      implicit none
      real*8 PI
      parameter (PI = 3.1415926535897932384d0)
      integer i, j, nptsx, nptsy, xdim, ydim
C      xdim and ydim are the absolute static dimensions.
C      nptsx, and nptsy are the (potentially dynamic) defined area of the 2D
C      arrays that is actually used.
      parameter (xdim=99, ydim=100, nptsx=35,nptsy=46)

      real*8 z(xdim, ydim), w(xdim, ydim), clevel(11),
     &  xg1(xdim), yg1(ydim),
     &  xg2(xdim, ydim), yg2(xdim, ydim)
      real*8 tr, xx, yy, argx, argy, distort

      common /plplot/ tr(6)

      data clevel /-1.d0, -0.8d0, -0.6d0, -0.4d0, -0.2d0,
     &  0.d0, 0.2d0, 0.4d0, 0.6d0 ,0.8d0, 1.d0/

      tr(1) = 2.d0/dble(nptsx-1)
      tr(2) = 0.0d0
      tr(3) = -1.0d0
      tr(4) = 0.0d0
      tr(5) = 2.d0/dble(nptsy-1)
      tr(6) = -1.0d0

C      Calculate the data matrices.
4     do i=1,nptsx
        xx = dble(i-1-(nptsx/2))/dble (nptsx/2)
        do j=1,nptsy
          yy = dble(j-1-(nptsy/2))/dble (nptsy/2) - 1.0d0
          z(i,j) = xx*xx - yy*yy
          w(i,j) = 2.d0*xx*yy
        enddo
      enddo

C      Build the 1-d coord arrays.
      distort = 0.4d0
      do i=1,nptsx
        xx = -1.d0 + dble(i-1)*2.d0/dble(nptsx-1)
        xg1(i) = xx + distort*cos(0.5d0*PI*xx)
      enddo

      do j=1,nptsy
        yy = -1.d0 + dble(j-1)*2.d0/dble(nptsy-1)
        yg1(j) = yy - distort*cos(0.5d0*PI*yy)
      enddo

C      Build the 2-d coord arrays.
      do i=1,nptsx
        xx = -1.d0 + dble(i-1)*2.d0/dble(nptsx-1)
        argx = 0.5d0*PI*xx
        do j=1,nptsy
          yy = -1.d0 + dble(j-1)*2.d0/dble(nptsy-1)
          argy = 0.5d0*PI*yy
          xg2(i,j) = xx + distort*cos(argx)*cos(argy)
          yg2(i,j) = yy - distort*cos(argx)*cos(argy)
        enddo
      enddo

C      Plot using identity transform
      call plenv(-1.0d0, 1.0d0, -1.0d0, 1.0d0, 0, 0)
      call plcol0(2)
      call plcont(z,xdim,ydim,1,nptsx,1,nptsy,clevel,11)
      call plstyl(1,1500,1500)
      call plcol0(3)
      call plcont(w,xdim,ydim,1,nptsx,1,nptsy,clevel,11)
      call plstyl(0,1500,1500)
      call plcol0(1)
      call pllab('X Coordinate', 'Y Coordinate',
     &  'Streamlines of flow')

      call plflush
      end
