/******************************************************************************/
/*                                                                            */
/*  MSETDIST - The msetdist algorithm for Mandelbrot and Julia sets           */
/*                                                                            */
/******************************************************************************/

#define STRICT
#include <windows.h>
#include <commctrl.h>
#include <assert.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <ctype.h>
#include <malloc.h>
#include <new.h>
#include <float.h>

#include "const.h"
#include "classes.h"
#include "extern.h"
#include "funcdefs.h"
#include "julia.rh"

double msetdist (
   COMPLEX *param ,   // c 
   COMPLEX *start ,   // Z0
   int max_dwell ,    // Dwell limit, should be as large as practical
   int julia ,        // Julia? (vs. Mandelbrot)
   double *orbit_r ,  // Work area max_dwell long
   double *orbit_i    // Ditto
   )
{
   int i, iter ;
   double real2, imag2, real, imag ;
   double rd, id, rtemp, itemp ;

   iter = 0 ;
   real = start->real ;
   imag = start->imag ;

   do {
      orbit_r[iter] = real ;
      orbit_i[iter] = imag ;
      real2 = real * real ;
      imag2 = imag * imag ;
      if (real2 + imag2 > 100000.0)
         break ;
      rtemp = 2. * real ;
      real = real2 - imag2 + param->real ;
      imag = rtemp * imag + param->imag ;
      } while (++iter < max_dwell) ;

   if (real2 + imag2 > 100000.0) {
      rd = 1.0 ;  // Recall when called for Mandelbrot, starts at c not 0
      id = 0.0 ;
      for (i=0 ; i<iter ; i++) {   /// Must not be <=, even though orbit[iter] is defined
         // (rtemp,itemp) is gradient at iter i+1
         rtemp = 2.0 * (orbit_r[i] * rd  -  orbit_i[i] * id) ;
         itemp = 2.0 * (orbit_r[i] * id  +  orbit_i[i] * rd) ;
         if (! julia)
            rtemp += 1.0 ;
         if ((fabs ( rtemp ) > 1.e30)  ||  (fabs ( itemp ) > 1.e30))
            return 0.0 ;
         rd = rtemp ;
         id = itemp ;
         }
      return log  ( real2 + imag2 ) *
             sqrt ( real2 + imag2 ) /
             sqrt ( rd * rd  +  id * id ) ;
      }

   return 0.0 ;
}