--[[-- Mandulia -- Mandelbrot/Julia explorer Copyright (C) 2010 Claude Heiland-Allen 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 . --]]-- 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