/* 
 * MDist - find out if a point is in the Mandelbrot set, and estimate
 *         the distance to the border.  That is, compute the radius of
 *         a disk that is entirely inside (or outside) the set.
 *
 * Copyright Niels Möller <nisse@lysator.liu.se>
 *
 * Released under the GNU General Public Licence
 *
 * $VER: MDist.c 1.2
 */

#include "MDist.h"
#include <assert.h>

#ifdef DBUG
#define bug printf
#define D(x) x 
#else !DBUG
#define D(x) 
#endif !DBUG

enum Distance {Small, Medium, Large};

MDist_t MDist(double point[2], double *radius, double *potential, 
	      struct MDistParameters *parameters)
{
  complex
    c = point[0] + 1i*point[1],
    z = c,
    /* partial derivative (d/dc) Z[n] */
    dc = 1;
  int
    /* iteration counter */
    n = 0,
    /* counter to bound the number of iterations when Rz>4 */
    nOutside = 0;
  enum Distance 
    /* 
     * Set to Small if the abs(Z[n])^2 <= 4.
     * Medium if 4 < abs(Z[n])^2 < MaxZAbs
     * Large if MaxZAbs < abs(Z[n]).
     */
    zSize = Small;
  MDist_t
    result;
  double
    /* Rz holds abs(z)^2 */
    Rz;

  assert(point);
  assert(radius);
  assert(potential);
  assert(parameters);

  D(bug("Enter MDist, c = %.15lf + i %.15lf.\n", RE(c), IM(c) ));

  *potential = *radius = 0.0; 
  /* 
   * Start iterating Zn and dZn/dc.
   * Z[n+1] = Z[n]*Z[n] + c
   * (d/dc) Z[n+1] = 2 * Z * (d/dc)Z[n] + 1
   */
  do
    {
      dc = 2*z*dc+1;
      z = SQR(z) + c;
      Rz = ABS2(z);

      if (Rz > parameters->MaxZAbs)
	zSize = Large;
      else
	if (Rz > 4)
	  zSize = Medium;
      n++;
      D(bug("Iteration %d: |z|^2 = %lf\n", n, Rz));
    }
  while ( ((zSize == Small) && (n < parameters->MaxIter)) ||
	 ((zSize == Medium) && (nOutside++ < parameters->MaxIterOutside) ));
				
  if (zSize == Small)
    {
      /* Find period, and compute radius of disk inside M */
      complex
	z0 = z,
	/* the derivative dz0/dc */
	dc0 = dc,
	/* more partial derivatives with respect to c and z0 */
	dz = 1.0,
	/* second order partial derivatives */
	dzz = 0.0,
	dcz = 0.0;
      double 
	/* squared distance from z[n] to z0 */
	diff;

      /*reset counter */
      n = 0;
      dc = 0.0;
      /* 
       * Iterate over one period, to compute length and attractivity of
       * a the periodic cycle.
       *
       * Z[n+1]           = Z[n]*Z[n] + c
       * (d/dc) Z[n+1]    = 2 * Z[n] * (d/dc)Z[n] + 1
       * (d/dz) Z[n+1]    = 2 * Z[n] * (d/dz)Z[n]
       * (d2/dcdz) Z[n+1] = 2 * (d/dz)Z[n] * (d/dc)Z[n] +
       *                    2 * Z[n] * (d2/dcdz)Z[n]
       * (d2/dz2) Z[n+1]  = 2 * SQR( (d/dz)Z[n] ) +
       *                    2 * (d2/dz2)Z[n]
       */

      do
	{
	  dcz = 2*(dc*dz+z*dcz);
	  dzz = 2*(SQR(dz)+z*dzz);
	  dc = 2*z*dc+1;
	  dz = 2*z*dz;
	  z = SQR(z) + c;
	  Rz = ABS2(z);
	  diff = ABS2(z-z0);
	  if (Rz > 4)
	    zSize = Medium;

	  D(bug("Z0= %lf + %lf i, Z%d = %lf + %lf i\n",
		RE(z0), IM(z0), n, RE(z), IM(z) ));

	}
      while ( (++n < parameters->MaxPeriod) && (zSize == Small) &&
	     (diff > (parameters->ERel * Rz + parameters->EAbs)));
      if ( (n == parameters->MaxPeriod) || (zSize > Small))
	{
	  result = ( (zSize == Small) ? InsideUnknown : OutsideUnknown) ;
	}
      else
	{
	  /* Estimate distance as
	   *                       | dZ[n] |2
	   *                   1 - |-------|
	   *   1                   | dZ[0] |
	   *  --- * ------------------------------------
	   *   4     |  d2Z[n]      d2Z[n]     dZ[0]  |
	   *         | --------- + -------- * ------- |
	   *         |  dcdZ[0]     dZ[n]2       dc   |
	   */

	  double 
	    A = 1 - ABS2(dz),
	    B = ABS2( dcz + dzz * dc0 );
	
	  D(bug("Period: %d\n",n));

	  if (A>0.0)
	    {
	      *radius = .25 * A / sqrt(B);
	      result = Inside;
	    }
	  else
	    {
	      result = InsideUnknown; /* Close to border, iteratives stayed small,
				       * but no attractive cycle were found */  
	    }
	}
    }
  else
    {
      /* 
       * Point was outside M. If Z[N] is large, compute potential:
       *
       *       -(N-1)
       * Gm = 2       * ln Z[N]
       *
       * and estimate distance to M:
       *
       *                      2(N-1)    Gm
       *   1     2 sinh Gm * 2       * e
       *  --- * ---------------------------- 
       *   4             | dZ[N] |
       *                 |-------|
       *                 |   dc  |
       *
       *       1     ln(Rz) * sqrt(Rz)
       *   ~= --- * ------------------
       *       4          |dZ/dc|
       *
       */

      if (zSize < Large)
	{
	  result = OutsideUnknown;
	}
      else
	{
	  double A = log(Rz);
	  D(bug("MDist: outside, |z|^2 = %lf, A = %lf\n", Rz, A));
	  *radius = .25 * A * sqrt(Rz / ABS2(dc));
	  *potential = ldexp(A,-n);
	  result = Outside;
	}
    }
  D(bug("Exit MDist.\n"));
  
  assert(*radius >= 0.0);
  assert(*potential >= 0.0);
  return result ;
}