!      Demo of multiple stream/window capability
!
!      Copyright (C) 2004  Arjen Markus
!      Copyright (C) 2004-2016  Alan W. Irwin
!
!      This file is part of PLplot.
!
!      PLplot is free software; you can redistribute it and/or modify
!      it under the terms of the GNU Library General Public License as
!      published by the Free Software Foundation; either version 2 of the
!      License, or (at your option) any later version.
!
!      PLplot is distributed in the hope that it will be useful,
!      but WITHOUT ANY WARRANTY; without even the implied warranty of
!      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!      GNU Library General Public License for more details.
!
!      You should have received a copy of the GNU Library General Public
!      License along with PLplot; if not, write to the Free Software
!      Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA

!     N.B. the pl_test_flt parameter used in this code is only
!     provided by the plplot module to allow convenient developer
!     testing of either kind(1.0) or kind(1.0d0) floating-point
!     precision regardless of the floating-point precision of the
!     PLplot C libraries.  We do not guarantee the value of this test
!     parameter so it should not be used by users, and instead user
!     code should replace the pl_test_flt parameter by whatever
!     kind(1.0) or kind(1.0d0) precision is most convenient for them.
!     For further details on floating-point precision issues please
!     consult README_precision in this directory.
!
! Plots several simple functions from other example programs.
!
! This version sends the output of the first 4 plots (one page) to two
! independent streams.

program x14f
    use plplot
    implicit none

    integer digmax

    character(len=80) :: driver
    character(len=15) :: geometry_master
    character(len=15) :: geometry_slave

    integer :: fam, num, bmax

    real(kind=pl_test_flt), dimension(101) ::  x, y
    real(kind=pl_test_flt), dimension(6) ::  xs, ys
    real(kind=pl_test_flt)  :: xscale, yscale, xoff, yoff
    common /plotdat/ x, y, xs, ys, xscale, yscale, xoff, yoff
    real(kind=pl_test_flt) :: xp0, yp0
    integer :: xleng0, yleng0, xoff0, yoff0
    integer :: plparseopts_rc, plsetopt_rc
    logical :: valid_geometry

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

    !      Process command-line arguments
    plparseopts_rc = plparseopts(PL_PARSE_FULL)
    if(plparseopts_rc .ne. 0) stop "plparseopts error"

    !      If valid geometry specified on command line, use it for both streams.
    call plgpage(xp0, yp0, xleng0, yleng0, xoff0, yoff0)
    valid_geometry = (xleng0 .gt. 0 .and. yleng0 .gt. 0)

    !      Set up first stream

    if(valid_geometry) then
        call plspage(xp0, yp0, xleng0, yleng0, xoff0, yoff0)
    else
        plsetopt_rc = plsetopt( 'geometry', geometry_master)
        if(plsetopt_rc .ne. 0) stop "plsetopt error"
    endif

    call plssub(2, 2)
    call plinit()

    call plgdev(driver)
    call plgfam(fam,num,bmax)

    write(6,'(3A)') 'Demo of multiple output streams via the ', &
           trim(driver), ' driver.'
    write(6,'(A)') 'Running with the second stream as slave '// &
           'to the first.'
    write(6,*)
    !      flush unit 6 so this part of stdout is guaranteed to be written prior
    !      to stdout generated by second plinit below.
    flush(6)

    !      Start next stream

    call plsstrm(1)

    if(valid_geometry) then
        call plspage(xp0, yp0, xleng0, yleng0, xoff0, yoff0)
    else
        plsetopt_rc = plsetopt( 'geometry', geometry_slave)
    endif

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

    call plspause(.false.)
    call plsdev(driver)
    call plsfam(fam,num,bmax)
    plsetopt_rc = plsetopt('fflen','2')
    call plinit()

    !      Set up the data & plot
    !      Original case

    call plsstrm(0)

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

    !      Set up the data & plot

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

    !      Set up the data & plot

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

    !      Set up the data & plot

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

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

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

    !      Back to master

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

    !      To slave

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

    !      Back to master to wait for user to advance

    call plsstrm(0)
    call pleop()

    !      Call plend to finish off.

    call plend()

contains

    !======================================================================

    subroutine plot1()
        use plplot
        implicit none

        real(kind=pl_test_flt), dimension(101) :: x, y
        real(kind=pl_test_flt), dimension(6) :: xs, ys
        real(kind=pl_test_flt) :: 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 * real(i,kind=pl_test_flt)/60.0_pl_test_flt
            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

        !      Set up the viewport and window using PLENV. The range in X is
        !      0.0 to 6.0, and the range in Y is 0.0 to 30.0. The axes are
        !      scaled separately (just = 0), and we just draw a labelled
        !      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' )

        !      Plot the data points

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

        !      Draw the line through the data

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

    !======================================================================

    subroutine plot2()
        use plplot
        implicit none
        real(kind=pl_test_flt), dimension(101) :: x, y
        real(kind=pl_test_flt), dimension(6) :: xs, ys
        real(kind=pl_test_flt) :: xscale, yscale, xoff, yoff
        integer :: i
        common /plotdat/ x, y, xs, ys, xscale, yscale, xoff, yoff

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

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

        !      Fill up the arrays

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

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

    end subroutine plot2
    !======================================================================

    subroutine plot3()
        !
        !      For the final graph we wish to override the default tick intervals,
        !      and so do not use_ PLENV

        use plplot, double_PI => PL_PI
        implicit none
        real(kind=pl_test_flt), parameter :: PI = double_PI
        real(kind=pl_test_flt), dimension(101) ::  x, y
        real(kind=pl_test_flt), dimension(6) :: xs, ys
        real(kind=pl_test_flt) :: xscale, yscale, xoff, yoff
        integer :: i
        common /plotdat/ x, y, xs, ys, xscale, yscale, xoff, yoff
        call pladv(0)

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

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

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

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

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

        call plstyl( (/1500/), (/1500/) )
        call plcol0(2)
        call plbox( 'g', 30.0_pl_test_flt, 0, 'g', 0.2_pl_test_flt, 0 )
        call plstyl( (/integer ::/), (/integer ::/) )

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

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

        call plcol0(4)
        call plline(x, y)
        call plflush
    end subroutine plot3

    !======================================================================

    subroutine plot4()

        use plplot, double_PI => PL_PI
        implicit none
        real(kind=pl_test_flt), parameter :: PI = double_PI
        character(len=3) :: text
        real(kind=pl_test_flt), dimension(0:360) :: x0, y0, x, y
        real(kind=pl_test_flt) :: dtr, theta, dx, dy, r
        integer :: i, j

        dtr = PI/180.0_pl_test_flt
        do i=0,360
            x0(i) = cos(dtr * real(i,kind=pl_test_flt))
            y0(i) = sin(dtr * real(i,kind=pl_test_flt))
        enddo

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

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

            !        Draw circles for polar grid

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

            !        Draw radial spokes for polar grid

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

            !        Write labels for angle

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

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

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

        !      Flush the plot at end

        call plflush
    end subroutine plot4

    !======================================================================

    integer function nsp(text)
        !      ==================

        !      Find first non-space character
        use plplot
        implicit none

        character*(*) text
        integer l

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

    !======================================================================

    subroutine plot5()

        use plplot
        implicit none
        integer :: i, j, nptsx, nptsy, xdim, ydim
        !      xdim and ydim are the absolute static dimensions.
        !      nptsx, and nptsy are the (potentially dynamic) defined area of the 2D
        !      arrays that is actually used.
        parameter (xdim=99, ydim=100, nptsx=35,nptsy=46)

        real(kind=pl_test_flt), dimension(xdim,ydim) ::  z, w
        real(kind=pl_test_flt), dimension(11) :: clevel
        real(kind=pl_test_flt) :: xx, yy
        real(kind=pl_test_flt), dimension(6) :: tr

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

        tr(1) = 2._pl_test_flt/real(nptsx-1,kind=pl_test_flt)
        tr(2) = 0.0_pl_test_flt
        tr(3) = -1.0_pl_test_flt
        tr(4) = 0.0_pl_test_flt
        tr(5) = 2._pl_test_flt/real(nptsy-1,kind=pl_test_flt)
        tr(6) = -1.0_pl_test_flt

        !      Calculate the data matrices.
        do i=1,nptsx
            xx = real(i-1-(nptsx/2),kind=pl_test_flt)/real(nptsx/2,kind=pl_test_flt)
            do j=1,nptsy
                yy = real(j-1-(nptsy/2),kind=pl_test_flt)/real(nptsy/2,kind=pl_test_flt) - 1.0_pl_test_flt
                z(i,j) = xx*xx - yy*yy
                w(i,j) = 2._pl_test_flt*xx*yy
            enddo
        enddo

        !      Plot using identity transform
        call plenv(-1.0_pl_test_flt, 1.0_pl_test_flt, -1.0_pl_test_flt, 1.0_pl_test_flt, 0, 0)
        call plcol0(2)
        call plcont(z,1,nptsx,1,nptsy,clevel,tr)
        call plstyl( (/1500/), (/1500/) )
        call plcol0(3)
        call plcont(w,1,nptsx,1,nptsy,clevel,tr)
        call plstyl( (/integer ::/), (/integer ::/) )
        call plcol0(1)
        call pllab('X Coordinate', 'Y Coordinate', &
               'Streamlines of flow')

        call plflush
    end subroutine plot5
end program x14f