! test_spiral.f90 find min and max for various ranges 

program test_spiral
  implicit none
  double precision :: xmin, xmax, ymin, ymax
  double precision :: zmin, zmax, atx, aty

  interface
    subroutine range(xmin, xmax, ymin, ymax, zmin, zmax, atx, aty)
      double precision, intent(in) :: xmin, xmax, ymin, ymax
      double precision, intent(out) :: zmin, zmax, atx, aty
    end subroutine range
  end interface

  print *, 'test_spiral.f90 running'
  xmin = -1.0
  xmax = 1.0
  ymin = -1.0
  ymax = 1.0
  call range(xmin, xmax, ymin, ymax, zmin, zmax, atx, aty)
  print *, '     xmin      xmax      ymin      ymax'
  print *, '     zmin      zmax       atx       aty'
  print '(2(4f10.6/))', xmin,xmax,ymin,ymax,zmin,zmax,atx,aty 

  xmin = -2.0
  xmax = 2.0
  ymin = -2.0
  ymax = 2.0
  call range(xmin, xmax, ymin, ymax, zmin, zmax, atx, aty)
  print '(2(4f10.6/),/)', xmin,xmax,ymin,ymax,zmin,zmax,atx,aty 

  xmin = -4.0
  xmax = 4.0
  ymin = -4.0
  ymax = 4.0
  call range(xmin, xmax, ymin, ymax, zmin, zmax, atx, aty)
  print '(2(4f10.6/),/)', xmin,xmax,ymin,ymax,zmin,zmax,atx,aty 

  xmin = 0.0
  xmax = 1.0
  ymin = 0.0
  ymax = 1.0
  call range(xmin, xmax, ymin, ymax, zmin, zmax, atx, aty)
  print '(2(4f10.6/),/)', xmin,xmax,ymin,ymax,zmin,zmax,atx,aty 

  xmin = -0.1
  xmax = 0.1
  ymin = -0.1
  ymax = 0.1
  call range(xmin, xmax, ymin, ymax, zmin, zmax, atx, aty)
  print '(2(4f10.6/),/)', xmin,xmax,ymin,ymax,zmin,zmax,atx,aty 

end program test_spiral

subroutine range(xmin, xmax, ymin, ymax, zmin, zmax, atx, aty)
  implicit none
  double precision, intent(in) :: xmin, xmax, ymin, ymax
  double precision, intent(out) :: zmin, zmax, atx, aty
  double precision :: tmin, tmax
  integer :: i, j
  double precision :: x, y, z, xh, yh
  integer :: n=100

  interface
    function spiral(x, y) result (z)
      double precision, intent(in) :: x, y
      double precision :: z
    end function spiral
  end interface

  xh = (xmax-xmin)/n
  yh = (ymax-ymin)/n
  z = spiral(xmin,ymin)
  tmin = z
  tmax = z
  atx = xmin
  aty = ymin

  do i=0,n
    x = xmin+i*xh
    do j=0,n
      y = ymin+j*yh
      z = spiral(x,y)
      if(z<tmin) then
        tmin = z
        atx = x
        aty = y
      end if
      if(z>tmax) tmax = z
    end do ! on j
  end do ! on i
  zmin = tmin
  zmax = tmax
end subroutine range