subroutine xstest(runid)

  !  test plots, output of xstraln library routines
  !  dmc apr 2011 initial version, only have test for xs_fxget routines which
  !  provide a fast (r,z) -> (rho,theta) "inverse" map.

  use xstraln_calls
  implicit none

  character*(*), intent(in) :: runid

  !-----------------------------------------------------------------
  logical :: lxout,lchar
  integer :: itty,lnout
  !-----------------------------------------------------------------

  itty = lnout(0)

  do
     if(lxout(0)) then
        write(itty,*) ' runid: '//trim(runid)
        write(itty,*) '  xstest (xstraln test) options:'
        write(itty,*) '   i -- (r,z) -> (rho,theta) fast inverse map test.'
        write(itty,*) '   q -- quit this menu, go back.'
     endif
     call uread('xstest option:$')

     if(lchar('q')) then
        exit

     else if(lchar('i')) then
        call xs_fxget_test(runid)

     else
        call uermsg('xstest: unknown option.',' ',0)
     endif
  enddo

end subroutine xstest

subroutine xs_fxget_test(runid)

  !  test plots:  xs_fxget routines which
  !  provide a fast (r,z) -> (rho,theta) "inverse" map.

  use xstraln_calls
  implicit none

  character*(*), intent(in) :: runid

  !-----------------------------------------------------------------
  logical :: lxout,lchar
  integer :: itty,lnout
  integer :: iwarn,ierr

  integer :: isize_grid = 201
  real*8, dimension(:), allocatable :: rgrid,zgrid
  real*8, dimension(:), allocatable :: rhowk,thwk,phwk,zwk
  real*8, dimension(:,:), allocatable :: rho2d,th2d,rho2d2,th2d2
  integer, dimension(:), allocatable :: iervec

  real*8 :: rmin,rmax,zmin,zmax

  real*8 :: rmin_loc,rmax_loc
  real*8 :: zmin_loc,zmax_loc
  real*8 :: arg1,arg2,argmin,argmax

  real*8, parameter :: zero = 0.0d0
  real*8, parameter :: one = 1.0d0
  !-----------------------------------------------------------------

  itty = lnout(0)

  call eq_rzlims(rmin,rmax,zmin,zmax,ierr)
  if(ierr.ne.0) return

  rmin_loc = rmin
  rmax_loc = rmax

  zmin_loc = zmin
  zmax_loc = zmax

  do 
     call alloc
     call gset

     if(lxout(0)) then
        write(itty,*) ' runid: '//trim(runid)
        write(itty,'(a,i4)') '   gridsize: ',isize_grid
        write(itty,'(a,2(1x,1pe11.4))') '   rmin,rmax (m): ',rmin_loc,rmax_loc
        write(itty,'(a,2(1x,1pe11.4))') '   zmin,zmax (m): ',zmin_loc,zmax_loc
        write(itty,*) '  xs_fxget_test options:'
        write(itty,*) '   n -- reset gridsize'
        write(itty,*) '   r -- reset rmin,rmax'
        write(itty,*) '   z -- reset zmin,zmax'
        write(itty,*) '   p -- **calculate map**, provide plotting options.'
        write(itty,*) '   q -- quit this menu, go back.'
     endif

     call uread('xs_fxget_test option:$')

     if(lchar('q')) then
        exit

     else if(lchar('n')) then
        if(lxout(0)) write(itty,*) ' enforced limits: min=11,max=9999'
        call ureadi('xs_fxget_test gridsize:$',isize_grid)
        isize_grid = min(9999,max(11,isize_grid))

     else if(lchar('r')) then
        call get_minmax('r',rmin,rmax,iwarn)
        if(iwarn.eq.0) then
           rmin_loc = argmin
           rmax_loc = argmax
        endif

     else if(lchar('z')) then
        call get_minmax('z',zmin,zmax,iwarn)
        if(iwarn.eq.0) then
           zmin_loc = argmin
           zmax_loc = argmax
        endif

     else if(lchar('p')) then
        call testplot

     else
        call uermsg('xs_fxget_test: unknown option.',' ',0)
     endif
  enddo
       
contains

  subroutine testplot

    ! call xs_fxget variants; plotting options

    integer :: jj,ii
    real*8 :: diff_rho,diff_th

    do jj=1,isize_grid
       zwk=zgrid(jj)
       iervec = 0

       call xs_fxget(isize_grid,one,phwk,rgrid,zwk,rho2d(:,jj),th2d(:,jj), &
            iervec)

       ierr=maxval(iervec)
       if(ierr.ne.0) exit

       do ii=1,isize_grid

          call xs_fxget(one,phwk(1),rgrid(ii),zgrid(jj), &
               rho2d2(ii,jj),th2d2(ii,jj),ierr)
          if(ierr.ne.0) exit

          diff_rho = abs(rho2d2(ii,jj)-rho2d(ii,jj))
          diff_th  = abs(th2d2(ii,jj)-th2d(ii,jj))

          if(max(diff_rho,diff_th).gt.1.0d-10) then
             write(itty,'(a,2(1x,1pe11.4),a,2(1x,1pe11.4))') ' %at (r,z)=', &
                  rgrid(ii),zgrid(jj),' found (diff_rho,diff_th) = ', &
                  diff_rho,diff_th
          endif
       enddo
       if(ierr.ne.0) exit
    enddo

    call r8_grf3f1(rgrid,rho2d,zgrid,zwk,isize_grid,isize_grid,isize_grid, &
         'r','rho','z','m',' ','m','xs_fxget results',runid,0)

    call r8_grf3f1(rgrid,th2d,zgrid,zwk,isize_grid,isize_grid,isize_grid, &
         'r','theta','z','m','rad','m','xs_fxget results',runid,0)

  end subroutine testplot

  subroutine get_minmax(var,minmin,maxmax,iwarn)
    character*(*), intent(in) :: var
    real*8, intent(in) :: minmin,maxmax
    integer, intent(out) :: iwarn

    !  get minmax pair, enforcing passed limits

    iwarn = 0

    if(lxout(0)) then
       write(itty,*) ' this routine enforces min < max and: '
       write(itty,'(a,a,a,1pe11.4)') '...minumum usable ',var,'min: ',minmin
       write(itty,'(a,a,a,1pe11.4)') '...maximum usable ',var,'max: ',maxmax
    endif

    call ureadr8('new '//var//'min (m):$',arg1)
    call ureadr8('new '//var//'max (m):$',arg2)

    if((arg1.lt.minmin).or.(arg1.gt.maxmax)) then
       write(itty,*) ' min/max limit enforced on "min" input.'
       arg1 = max(minmin, min(maxmax, arg1))
    endif

    if((arg2.lt.minmin).or.(arg2.gt.maxmax)) then
       write(itty,*) ' min/max limit enforced on "max" input.'
       arg2 = max(minmin, min(maxmax, arg2))
    endif

    argmin = min(arg1,arg2)
    argmax = max(arg1,arg2)

    if(argmin.eq.argmax) then
       write(itty,*) ' warning, min = max not supported...'
       iwarn=1
    endif
  end subroutine get_minmax

  subroutine alloc

    integer :: isize_now

    if(.not.allocated(rgrid)) then
       isize_now = 0
    else
       isize_now = size(rgrid)
    endif

    if(isize_now.ne.isize_grid) then
       if(allocated(rgrid)) then
          deallocate(rgrid,zgrid)
          deallocate(rhowk,thwk,phwk,zwk)
          deallocate(rho2d,th2d,rho2d2,th2d2)
          deallocate(iervec)
       endif

       allocate(rgrid(isize_grid),zgrid(isize_grid))
       allocate(rhowk(isize_grid),thwk(isize_grid))
       allocate(phwk(isize_grid),zwk(isize_grid))
       allocate(rho2d(isize_grid,isize_grid))
       allocate(th2d(isize_grid,isize_grid))
       allocate(rho2d2(isize_grid,isize_grid))
       allocate(th2d2(isize_grid,isize_grid))
       allocate(iervec(isize_grid))
    endif

  end subroutine alloc

  subroutine gset

    integer :: ii

    do ii=1,isize_grid
       rgrid(ii) = ((ii-1)*rmax_loc + (isize_grid-ii)*rmin_loc)/ &
            (isize_grid - 1)
       zgrid(ii) = ((ii-1)*zmax_loc + (isize_grid-ii)*zmin_loc)/ &
            (isize_grid - 1)

       rhowk(ii) = zero
       thwk(ii) = zero
       phwk(ii) = zero
       zwk(ii) = zero

       rho2d(:,ii) = zero
       th2d(:,ii) = zero
       rho2d2(:,ii) = zero
       th2d2(:,ii) = zero

       iervec(ii) = 0
    enddo

  end subroutine gset

end subroutine xs_fxget_test
