! spiral.f90 spiral trough 

function spiral(x, y) result (z)
  implicit none
  double precision, intent(in) :: x, y
  double precision :: z
  double precision :: halfpi = 1.57079
  double precision :: r, theta, tz

  r = sqrt(x*x+y*y)
  if( r < 0.001) then
    theta = 0.0
  else
    theta = atan2(y,x)
  end if
  tz = theta/halfpi
  if(tz < 0.0) tz = tz + 4.0
  z = 0.75*r + sin((r+tz)*halfpi)
end function spiral