--[[--

Mandulia -- Mandelbrot/Julia explorer
Copyright (C) 2010  Claude Heiland-Allen <claudiusmaximus@goto10.org>

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program 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 General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.

--]]--

do

  --[[--

  SYNOPSIS

    s = exteriordistance(x, y, n, r)

  INPUTS

    x, y  :  point in the complex plane
    n     :  iteration limit
    r     :  escape radius

  OUTPUTS

    s     :  distance from Mandelbrot set (or nil if inside)

  EXAMPLE

    s = exteriordistance(x, y, 100, 100)
    if s ~= nil and s < 0.01 then ... else ... end

  --]]--
  function exteriordistance(x, y, n, r)

    local px = 0
    local py = 0
    local cx = x
    local cy = y
    local dx = 1
    local dy = 0
    local r2 = r * r
    local escape = false
    local px1, py1, dx1, dy1, px2, py2, pxy, d2, p, d

    for i = 1,n do
      -- p_n+1 := p_n * p_n + c
      -- d_n+1 := 2 * p_n * d_n + 1
      px2 = px * px
      py2 = py * py
      d2 = px2 + py2
      if d2 > r2 then escape = true ; break end
      pxy = px * py
      px1 = px2 - py2 + cx
      py1 = 2 * pxy + cy
      dx1 = 2 * (px * dx - py * dy) + 1
      dy1 = 2 * (dx * py + dy * px)
      px = px1 ; py = py1 ; dx = dx1 ; dy = dy1
    end

    if escape then
      p = math.sqrt(px * px + py * py)
      d = math.sqrt(dx * dx + dy * dy)
      return 2 * p * math.log(p) / d
    else
      return nil
    end

  end

end