/* * brat.c * * FUNCTION: * Explore Hausdorf measure of mandelbrot set. * * HISTORY: * quick hack -- Linas Vepstas October 1989 * modernize -- Linas Vepstas March 1996 */ #include #include #include #include #include /*-------------------------------------------------------------------*/ /* This routine computes a convergent for the Mandelbrot set iterator. */ void mandelbrot_cutoff ( float *glob, unsigned int sizex, unsigned int sizey, double re_center, double im_center, double width, int itermax) { unsigned int i,j, globlen; double re_start, im_start, delta; double re_position, im_position; double re, im, tmp, mod, modulus=0.0; double dre, dim, dmod; double ddre, ddim, ddmod; double zpre, zpim, zppre, zppim; double mp, mpp; double phi, phip, phipp; int loop; double omod=0.0, frac, theta; double escape_radius = 1.0e30; double ren, tl; double tau; double *regulator, *rp, *rpp, *rppp, *rpppp; double sum_n, sum_np, sum_npp, sum_nppp, sum_npppp; double sum_re, sum_im, sum_mod; double sum_rep, sum_imp, sum_modp; double sum_repp, sum_impp, sum_modpp; double sum_dre, sum_dim, sum_dmod; double sum_ddre, sum_ddim, sum_ddmod; double sum_ddrep, sum_ddimp, sum_ddmodp; double sum_ddrepp, sum_ddimpp, sum_ddmodpp; double sum_zpre, sum_zpim, sum_zpmod; double sum_zppre, sum_zppim, sum_zppmod; /* first, compute the regulator, so that the itermax'th iteration makes * a negligable contribution (about 1e-30) */ tau = 16.0 / ((double) itermax); /* set up smooth ramp * regulator is exponential, and rp is derivative w.r.t. tau, * rpp is second deriv. w.r.t. tau */ sum_n = sum_np = sum_npp = sum_nppp = sum_npppp = 0.0; regulator = (double *) malloc ((itermax+1)*sizeof (double)); rp = (double *) malloc ((itermax+1)*sizeof (double)); rpp = (double *) malloc ((itermax+1)*sizeof (double)); rppp = (double *) malloc ((itermax+1)*sizeof (double)); rpppp = (double *) malloc ((itermax+1)*sizeof (double)); for (i=0; i escape_radius*escape_radius) break; } #if 0 modulus = sqrt (modulus); modulus = sqrt (sum_re*sum_re + sum_im*sum_im); frac = log (log (modulus)) * tl; glob [i*sizex +j] = modulus / sum_n; glob [i*sizex +j] = sum_mod / sum_n; glob [i*sizex +j] = sum_mod - modulus; /* the following computes an almost-flat, divergence free thing */ glob [i*sizex +j] = modulus / sum_n; glob [i*sizex +j] /= (sqrt (re_position*re_position+im_position*im_position)); /* the interesting one is the z-prime-prime */ modulus = sqrt (sum_ddre*sum_ddre + sum_ddim*sum_ddim); glob [i*sizex +j] = modulus / sum_n; // glob [i*sizex +j] -= 0.25 * exp( -0.75 * log((re_position-0.25)*(re_position-0.25)+im_position*im_position)); /* --------------------------------------------------------- */ /* the interesting one is the z-prime-prime */ /* here we use a taylor expansion to extrapolate to tau=0 */ /* first, we need the derivatives of modulus w.r.t tau */ modulus = sqrt (sum_ddre*sum_ddre + sum_ddim*sum_ddim); mp = (sum_ddrep * sum_ddre + sum_ddimp * sum_ddim) / modulus; mpp = sum_ddrep * sum_ddrep + sum_ddre * sum_ddrepp; mpp += sum_ddimp * sum_ddimp + sum_ddim * sum_ddimpp - mp*mp; mpp /= modulus; /* next, we need derivatives of m/n w.r.t tau */ mod = modulus / sum_n; dmod = (mp - mod * sum_np) / sum_n; ddmod = (mpp - 2.0 * dmod * sum_np - mod * sum_npp) / sum_n; /* finally the taylor expansion */ glob [i*sizex +j] = mod - tau* (dmod - 0.5 * tau * ddmod); // glob [i*sizex +j] -= 0.25 * exp( -0.75 * log((re_position-0.25)*(re_position-0.25)+im_position*im_position)); /* --------------------------------------------------------- */ /* OK, lets do the taylor expansion for just-plain Z */ modulus = sqrt (sum_re*sum_re + sum_im*sum_im); mp = (sum_rep * sum_re + sum_imp * sum_im) / modulus; mpp = sum_rep * sum_rep + sum_re * sum_repp; mpp += sum_imp * sum_imp + sum_im * sum_impp - mp*mp; mpp /= modulus; /* next, we need derivatives of m/n w.r.t tau */ mod = modulus / sum_n; dmod = (mp - mod * sum_np) / sum_n; ddmod = (mpp - 2.0 * dmod * sum_np - mod * sum_npp) / sum_n; /* finally the taylor expansion */ glob [i*sizex +j] = mod - tau* (dmod - 0.5 * tau * ddmod); // printf ("%9.6g %9.6g %9.6g\n", re_position, glob[i*sizex+j], mod); /* ok, this part should be the divergent part ... */ theta = 0.5 * atan2 (-im_position, 0.25-re_position); mod = (re_position-0.25)*(re_position-0.25)+im_position*im_position; mod = pow (mod, 0.25); re = - mod * cos(theta); im = - mod * sin(theta); re += 0.5; tmp = sqrt(re*re+im*im); if (0.5 < tmp) tmp = 0.5; glob [i*sizex +j] -= tmp; /* --------------------------------------------------------- */ /* the interesting one is the z-prime-prime */ /* here we use a taylor expansion to extrapolate to tau=0 */ /* first, we need the drerivatives of modulus w.r.t tau */ modulus = sqrt (sum_ddre*sum_ddre + sum_ddim*sum_ddim); mp = (sum_ddrep * sum_ddre + sum_ddimp * sum_ddim) / modulus; mpp = sum_ddrep * sum_ddrep + sum_ddre * sum_ddrepp; mpp += sum_ddimp * sum_ddimp + sum_ddim * sum_ddimpp - mp*mp; mpp /= modulus; /* finally the taylor expansion */ /* subtract the main-body divergence */ tmp = 0.25 * exp( -0.75 * log((re_position-0.25)*(re_position-0.25)+im_position*im_position)); glob [i*sizex +j] = (modulus-tmp*sum_n); glob [i*sizex +j] -= tau* ((mp-tmp*sum_np) - 0.5 * tau * (mpp-tmp*sum_npp)); #endif /* --------------------------------------------------------- */ /* OK, lets do the taylor expansion for just-plain Z */ /* here, we subtract the leading divergence */ modulus = sqrt (sum_re*sum_re + sum_im*sum_im); mp = (sum_rep * sum_re + sum_imp * sum_im) / modulus; mpp = sum_rep * sum_rep + sum_re * sum_repp; mpp += sum_imp * sum_imp + sum_im * sum_impp - mp*mp; mpp /= modulus; /* finally the taylor expansion */ glob [i*sizex +j] = modulus - tau* (mp - 0.5 * tau * mpp); theta = 0.5 * atan2 (-im_position, 0.25-re_position); mod = (re_position-0.25)*(re_position-0.25)+im_position*im_position; mod = pow (mod, 0.25); re = - mod * cos(theta); im = - mod * sin(theta); re += 0.5; tmp = sqrt(re*re+im*im); if (0.5 < tmp) tmp = 0.5; glob [i*sizex +j] -= tmp * (sum_n - tau* (sum_np - 0.5 * tau * sum_npp)); #if 0 /* --------------------------------------------------------- */ /* the interesting one is the z-prime-prime */ /* here we use a taylor expansion to extrapolate to tau=0 */ /* first, we need the derivatives of modulus w.r.t tau */ /* we compute the phase */ phi = atan2 (sum_ddim, sum_ddre); modulus = 1.0 / sqrt (sum_ddre*sum_ddre + sum_ddim*sum_ddim); phip = (sum_ddre * sum_ddimp - sum_ddrep * sum_ddim) * modulus; phipp = (sum_ddre * sum_ddimp - sum_ddrep * sum_ddim) * modulus; phipp *= -2.0* (sum_ddre * sum_ddrep + sum_ddim * sum_ddimp) * modulus; phipp += (sum_ddre * sum_ddimpp - sum_ddrepp * sum_ddim) * modulus; glob [i*sizex +j] = (phi + M_PI)/(2.0*M_PI); // glob [i*sizex +j] = (phi - tau* (phip - 0.5 * tau * phipp) +M_PI)/(2.0*M_PI); theta = -1.5 * atan2 (-im_position, 0.25-re_position); mod = (re_position-0.25)*(re_position-0.25)+im_position*im_position; mod = pow (mod, -0.75); re = 0.25 * mod * cos(theta); im = 0.25 * mod * sin(theta); glob [i*sizex +j] = (sum_ddre/sum_n-re)*(sum_ddre/sum_n-re); glob [i*sizex +j] += (sum_ddim/sum_n-im)*(sum_ddim/sum_n-im); glob [i*sizex +j] = sqrt (glob[i*sizex+j]); /* --------------------------------------------------------- */ theta = 0.5 * atan2 (-im_position, 0.25-re_position); mod = (re_position-0.25)*(re_position-0.25)+im_position*im_position; mod = pow (mod, 0.25); re = - mod * cos(theta); im = - mod * sin(theta); re += re_position; im += im_position; glob [i*sizex +j] = sqrt (re*re +im*im); #endif /* --------------------------------------------------------- */ re_position += delta; } im_position -= delta; /*top to bottom, not bottom to top */ } } /*-------------------------------------------------------------------*/ /* this routine fills in the exterior of the mandelbrot set using */ /* the classic algorithm */ void mandelbrot_out ( float *glob, unsigned int sizex, unsigned int sizey, double re_center, double im_center, double width, int itermax) { unsigned int i,j, globlen; double re_start, im_start, delta; double re_position, im_position; double re, im, tmp; int loop; double modulus=0.0, frac; double escape_radius = 3.1; double ren, tl; ren = log( log (escape_radius)) / log(2.0); tl = 1.0/ log(2.0); delta = width / (double) sizex; re_start = re_center - width / 2.0; im_start = im_center + width * ((double) sizey) / (2.0 * (double) sizex); globlen = sizex*sizey; for (i=0; i escape_radius*escape_radius) break; } modulus = sqrt (modulus); frac = log (log (modulus)) *tl; /* frac = Re (c/z*z) */ tmp = (re*re-im*im); im = 2.0*re*im; re = tmp; frac = re*re_position + im*im_position; frac /= re*re+im*im; if (0.0 > frac) frac = -frac; #ifdef JUNK /* second order correction */ frac = frac * ( 2.0 - frac + 2.0*ren); nprev = loop - nprev; prev = frac - prev - (double)nprev; /* printf ("ij %d %d deltan %d delta frac %f \n", i, j, nprev, prev); */ nprev = loop; prev = frac; glob [i*sizex +j] = ((double) loop) - frac; #endif glob [i*sizex +j] = frac; glob [i*sizex +j] = ((double) loop); /* glob [i*sizex +j] = ((float) (loop%10)) / 10.0; if (loop == itermax) { glob[i*sizex+j] = 0.0; } else {glob[i*sizex+j]=0.9999;} */ re_position += delta; } im_position -= delta; /*top to bottom, not bottom to top */ } } /*-------------------------------------------------------------------*/ /* this routine fills in the interior and exterior of the mandelbrot set using * the classic algorithm. The derivative (infinitessimal flow) * is used to play games. */ void dmandelbrot_out ( float *glob, unsigned int sizex, unsigned int sizey, double re_center, double im_center, double width, int itermax) { unsigned int i,j, globlen; double re_start, im_start, delta; double re_position, im_position; double re, im, tmp; double dre, dim; double ddre, ddim; double d3re, d3im; double d4re, d4im; int loop; double modulus, phi, phip, phi3, frac; double escape_radius = 3450.0; double ren, tl; ren = log( log (escape_radius)) / log(2.0); tl = 1.0/ log(2.0); delta = width / (double) sizex; re_start = re_center - width / 2.0; im_start = im_center + width * ((double) sizey) / (2.0 * (double) sizex); globlen = sizex*sizey; for (i=0; i escape_radius*escape_radius) break; } modulus = (re*re + im*im); modulus = sqrt (modulus); frac = log (log (modulus)) *tl; frac = ((double) loop) - frac + 1.0; /* compute zprime/z */ tmp = re*dre + im*dim; /* divergence */ dim = re*dim - im*dre; /* curl */ dre = tmp; dre /= (re*re + im*im); dim /= (re*re + im*im); /* compute zprimeprime/z */ tmp = re*ddre + im*ddim; /* divergence */ ddim = re*ddim - im*ddre; /* curl */ ddre = tmp; ddre /= (re*re + im*im); ddim /= (re*re + im*im); /* compute z'''/z */ tmp = re*d3re + im*d3im; /* divergence */ d3im = re*d3im - im*d3re; /* curl */ d3re = tmp; d3re /= (re*re + im*im); d3im /= (re*re + im*im); /* compute z'(4)/z */ tmp = re*d4re + im*d4im; /* divergence */ d4im = re*d4im - im*d4re; /* curl */ d4re = tmp; d4re /= (re*re + im*im); d4im /= (re*re + im*im); /* phase */ phi = atan2 (dim, dre); phi += M_PI; phi /= 2.0*M_PI; phip = atan2 (ddim, ddre); phip += M_PI; phip /= 2.0*M_PI; phi3 = atan2 (d3im, d3re); phi3 += M_PI; phi3 /= 2.0*M_PI; modulus = sqrt (dre*dre+dim*dim); modulus = sqrt (ddre*ddre+ddim*ddim); // modulus = sqrt (d4re*d4re+d4im*d4im); modulus /= (double) loop; modulus *= log((double) loop); modulus *= log((double) loop); modulus = sqrt (dre*dre+dim*dim); glob [i*sizex +j] = modulus; re_position += delta; } im_position -= delta; /*top to bottom, not bottom to top */ } } /*-------------------------------------------------------------------*/ /* this routine attempts to determine if * a given point is inside or outside the mandelbrot set. */ void mandelbrot_decide ( float *glob, unsigned int sizex, unsigned int sizey, double re_center, double im_center, double width, int itermax) { unsigned int i,j, globlen; double re_start, im_start, delta; double re_position, im_position; double re, im, tmp; double dre, dim; double ddre, ddim; double zppre, zppim; int loop; double modulus, limit; double escape_radius = 50.0; int state; double *limits; clock_t start, stop; int hunds; #define OUTSIDE 1 #define UNDECIDED 6 #define INSIDE 4 #define ERROR 8 delta = width / (double) sizex; re_start = re_center - width / 2.0; im_start = im_center + width * ((double) sizey) / (2.0 * (double) sizex); globlen = sizex*sizey; for (i=0; i escape_radius*escape_radius) { if ((UNDECIDED == state) || (OUTSIDE == state)) { state = OUTSIDE; } else { state = ERROR; if (ERROR < loop) state = loop; } break; } /* compute zprimeprime/z */ zppre = re*ddre + im*ddim; /* divergence */ zppim = re*ddim - im*ddre; /* curl */ zppre /= (re*re + im*im); zppim /= (re*re + im*im); modulus = sqrt (zppre*zppre+zppim*zppim); // modulus /= (double) loop; limit = limits[loop]; /* check to see if we can identify this as an interior point */ if ((10 < loop) && (limit > modulus) && (UNDECIDED == state)) { state = INSIDE; // record time saved ... // state = loop; break; } } if (ERROR > state) { glob [i*sizex +j] = ((double) state) / 8.1; } else { glob [i*sizex +j] = ((double) state) / ((double)itermax); } re_position += delta; /* Measure cpu time spent. On Linux, there are a million * clocks per sec, which will overflow in about an hour. * so we have to measure time frequently! */ if (j%100==0) { stop = clock(); hunds += (stop-start) / (CLOCKS_PER_SEC/100); start = stop; } } im_position -= delta; /*top to bottom, not bottom to top */ } stop = clock(); hunds += (stop-start) / (CLOCKS_PER_SEC/100); start = stop; free (limits); /* do some counting. */ { int inside=0, uncertain=0; double rem; for (i=0; i 8.1*glob[i]) && (((double) UNDECIDED)-0.1 < 8.1*glob[i])) { uncertain ++; } if ((((double) INSIDE)+0.1 > 8.1*glob[i]) && (((double) INSIDE)-0.1 < 8.1*glob[i])) {inside ++; } } rem = ((double) uncertain) / ((double)(uncertain+inside)); printf (">> %d %d %d %d %f %d\n", itermax, uncertain, inside, uncertain+inside, rem, hunds); } } /*-------------------------------------------------------------------*/ /* this routine fills in the exterior of the mandelbrot set using */ /* the classic algorithm */ void mandelbrot_inverse ( float *glob, unsigned int sizex, unsigned int sizey, double re_center, double im_center, double width, int itermax) { unsigned int i,j, globlen; double re_start, im_start, delta; double re_position, im_position; double re, im, tmp; int loop; int ire, jim; double modulus=0.0, phase, tp, frac; double escape_radius = 2000.1; double ren, tl; ren = log( log (escape_radius)) / log(2.0); tl = 1.0/ log(2.0); delta = width / (double) sizex; delta *= 0.25; re_start = re_center - width / 2.0; im_start = im_center + width * ((double) sizey) / (2.0 * (double) sizex); globlen = sizex*sizey; for (i=0; i escape_radius*escape_radius) break; } if (loop == itermax) continue; modulus = sqrt (modulus); frac = ren - log (log (modulus)) *tl; frac = ((double) loop) - ren + frac; modulus = log (modulus); phase = atan2 (im, re); tp = 1.0 / pow (2.0, loop); modulus *= tp; phase *= tp; modulus = exp (modulus); re = cos (phase) * modulus; im = sin (phase) * modulus; re -= 1.0; im /= re; re *= 24.0; im *= 1.0; ire = sizex * (re + 0.5); if (0 > ire) ire = 0; if (sizex <= ire) ire = sizex-1; im *= ((double) sizey)/ ((double) sizex); jim = sizey * (im + 0.5); if (0 > jim) jim = 0; if (sizey <= jim) jim = sizey-1; glob [jim*sizex +ire] ++; re_position += delta; } im_position -= delta; /*top to bottom, not bottom to top */ } } /*-------------------------------------------------------------------*/ /* this routine fills in the interior of the mandelbrot set using */ /* the classic algorithm */ void mandelbrot_stop ( float *glob, unsigned int sizex, unsigned int sizey, double re_center, double im_center, double width, int itermax) { unsigned int i,j, globlen; double re_start, im_start, delta; double re_position, im_position; double re, im, tmp; int loop; delta = width / (double) sizex; re_start = re_center - width / 2.0; im_start = im_center + width * ((double) sizey) / (2.0 * (double) sizex); globlen = sizex*sizey; for (i=0; i 154.0) break; } /* glob [i*sizex +j] = (2.0+ re)/2.5; */ /* glob [i*sizex +j] = re; */ glob [i*sizex +j] = im; re_position += delta; } im_position -= delta; /*top to bottom, not bottom to top */ } } /*-------------------------------------------------------------------*/ /* this routine fills in the interior of the mandelbrot set using */ /* the Cliff Pickover's stalks algorithm */ void mandelbrot_stalk ( float *glob, unsigned int sizex, unsigned int sizey, double re_center, double im_center, double width, int itermax, double stalkx, double stalky) { unsigned int i,j, globlen; double re_start, im_start, delta; double re_position, im_position; double re, im, tmp, tmpx, tmpy; double visited_x, visited_y; int loop; delta = width / (double) sizex; re_start = re_center - width / 2.0; im_start = im_center + width * ((double) sizey) / (2.0 * (double) sizex); globlen = sizex*sizey; for (i=0; i 154.0) break; } visited_x = re; visited_y = im; tmp = re*re - im*im + re_position; im = 2.0*re*im + im_position; re = tmp; #endif MAXDIST for (loop=1; loop 154.0) break; #ifdef HAND_BUILT_COLORMAP Hand-built colormap -- visually bad idea if (((-0.5 < re) && (re < 0.5)) ||((-0.5 < im) && (im < 0.5))) if (0.1 > glob [i*sizex +j]) glob [i*sizex +j] = 0.1; if (((-0.2 < re) && (re < 0.2)) ||((-0.2 < im) && (im < 0.2))) if (0.2 > glob [i*sizex +j]) glob [i*sizex +j] = 0.2; if (((-0.1 < re) && (re < 0.1)) ||((-0.1 < im) && (im < 0.1))) if (0.3 > glob [i*sizex +j]) glob [i*sizex +j] = 0.3; if (((-0.05 < re) && (re < 0.05)) ||((-0.05 < im) && (im < 0.05))) if (0.5 > glob [i*sizex +j]) glob [i*sizex +j] = 0.5; if (((-0.02 < re) && (re < 0.02)) ||((-0.02 < im) && (im < 0.02))) if (0.6 > glob [i*sizex +j]) glob [i*sizex +j] = 0.6; if (((-0.01 < re) && (re < 0.01)) ||((-0.01 < im) && (im < 0.01))) if (0.7 > glob [i*sizex +j]) glob [i*sizex +j] = 0.7; if (((-0.005 < re) && (re < 0.005)) ||((-0.005 < im) && (im < 0.005))) if (0.8 > glob [i*sizex +j]) glob [i*sizex +j] = 0.8; if (((-0.002 < re) && (re < 0.002)) ||((-0.002 < im) && (im < 0.002))) if (0.9 > glob [i*sizex +j]) glob [i*sizex +j] = 0.9; if (((-0.001 < re) && (re < 0.001)) ||((-0.001 < im) && (im < 0.001))) if (1.0 > glob [i*sizex +j]) glob [i*sizex +j] = 1.0; #endif /* HAND_BUILT_COLORMAP */ #ifdef STALK tmpx = re-stalkx; if (0.0 > tmpx) tmpx = -tmpx; tmpy = im-stalky; if (0.0 > tmpy) tmpy = -tmpy; ltmpx = -log (tmpx); if (ltmpx > glob [i*sizex +j]) glob [i*sizex +j] = ltmpx; ltmpy = -log (tmpy); if (ltmpy > glob [i*sizex +j]) glob [i*sizex +j] = ltmpy; #endif /* STALK */ #ifdef CROSS tmpx = re-stalkx; if (0.0 > tmpx) tmpx = -tmpx; tmpy = im-stalky; if (0.0 > tmpy) tmpy = -tmpy; ltmpx = -log (tmpx); if (0.2 > tmpy) { if (ltmpx > glob [i*sizex +j]) glob [i*sizex +j] = ltmpx; } ltmpy = -log (tmpy); if (0.2 > tmpx) { if (ltmpy > glob [i*sizex +j]) glob [i*sizex +j] = ltmpy; } #endif /* CROSS */ #define CIRCSTALK #ifdef CIRCSTALK tmpx = re - stalkx; tmpy = im - stalky; tmp = tmpx*tmpx + tmpy*tmpy -0.04; tmp = tmpx*tmpx + tmpy*tmpy; if (0.0 > tmp) tmp = -tmp; tmp = -log (tmp); if (tmp > glob [i*sizex +j]) glob [i*sizex +j] = tmp; #endif CIRCSTALK #ifdef MAXDIST tmpx = re-stalkx - visited_x; tmpy = im-stalky - visited_y; tmpx = re - stalkx; tmpy = im - stalky; tmp = tmpx*tmpx + tmpy*tmpy; tmp = -log (tmp); tmpx = re - stalkx-0.001; tmpy = im - stalky-0.001; tmp += log (tmpx*tmpx + tmpy*tmpy); if (tmp > glob [i*sizex +j]) glob [i*sizex +j] = tmp; #endif MAXDIST } re_position += delta; } im_position -= delta; /*top to bottom, not bottom to top */ } } /*-------------------------------------------------------------------*/ /* this routine fills in the exterior of the circle map set using */ /* the classic algorithm */ void circle_out ( float *glob, unsigned int sizex, unsigned int sizey, double re_center, double im_center, double width, int itermax, double time) { unsigned int i,j, globlen; double re_start, im_start, delta; double re_position, im_position; double re_omega, im_omega; double re_K, im_K; double re, im, tmp; int n, loop; double ep, em, es, ec; delta = width / (double) sizex; re_start = re_center - width / 2.0; im_start = im_center + width * ((double) sizey) / (2.0 * (double) sizex); globlen = sizex*sizey; for (i=0; i 1.0e6) break; n = (int) re; if (0.0 > re) n--; re -= (double) n; if (im*im > 2284.0) break; } glob [i*sizex +j] = loop; re_position += delta; } im_position -= delta; /*top to bottom, not bottom to top */ } tmp = 1.0 / ((double) itermax); for (i=0; i= isamp) isamp = 1; isamp *= globlen*itermax; irow = isamp / sizey; for (i=0; i 1.0e6) break; n = (int) re; if (0.0 > re) n--; re -= (double) n; if (im*im > 2284.0) break; horiz_pix = (int) (x_slope * re - x_off); vert_pix = (int) (y_slope * im - y_off); if ( (horiz_pix >= 0) && (horiz_pix < sizex) && (vert_pix >= 0) && (vert_pix < sizey) && (SETTLE_COUNT < loop)) { glob [vert_pix*sizex + horiz_pix] ++; } } } /* renormalize */ tmp = 1.0 / ((double) itermax*(LOOP_COUNT-SETTLE_COUNT)) ; for (i=0; i= isamp) isamp = 1; isamp *= globlen*itermax; irow = isamp / sizey; for (i=0; i 7.0) break; } for (loop=0; loop < LOOP_COUNT; loop++) { tmp = re*re - im*im + re_position; im = 2.0*re*im + im_position; re = tmp; if ((re*re + im*im) > 7.0) break; horiz_pix = (int) (x_slope * re - x_off); vert_pix = (int) (y_slope * im - y_off); if ( (horiz_pix >= 0) && (horiz_pix < sizex) && (vert_pix >= 0) && (vert_pix < sizey)) { glob [vert_pix*sizex + horiz_pix] ++; } } } /* renormalize */ for (i=0; i= isamp) isamp = 1; isamp *= globlen*itermax; irow = isamp / sizey; for (i=0; i 7.0) break; } for (loop=0; loop < LOOP_COUNT; loop++) { tmp = re*re - im*im + re_position; im = 2.0*re*im + im_position; re = tmp; if ((re*re + im*im) > 7.0) break; horiz_pix = (int) (x_slope * (re-re_position) - x_off); vert_pix = (int) (y_slope * (im-im_position) - y_off); if ( (horiz_pix >= 0) && (horiz_pix < sizex) && (vert_pix >= 0) && (vert_pix < sizey)) { glob [vert_pix*sizex + horiz_pix] ++; } } } /* renormalize */ for (i=0; i= isamp) isamp = 1; isamp *= globlen*itermax; irow = isamp / sizey; for (i=0; i 7.0) break; horiz_pix = (int) (x_slope * re - x_off); vert_pix = (int) (y_slope * im - y_off); if ( (horiz_pix >= 0) && (horiz_pix < sizex) && (vert_pix >= 0) && (vert_pix < sizey)) { glob [vert_pix*sizex + horiz_pix] += loop; square [vert_pix*sizex + horiz_pix] += loop*loop; cube [vert_pix*sizex + horiz_pix] += loop*loop*loop; quart [vert_pix*sizex + horiz_pix] += loop*loop*loop*loop; quint [vert_pix*sizex + horiz_pix] += loop*loop*loop*loop*loop; count [vert_pix*sizex + horiz_pix] ++; } } } /* renormalize */ ollie = 1.0 / ((double) LOOP_COUNT); for (i=0; i= isamp) isamp = 1; isamp *= globlen*itermax; irow = isamp / sizey; for (i=0; i 7.0) break; if ((re[1]*re[1] + im[1]*im[1]) > 7.0) break; if ((re[2]*re[2] + im[2]*im[2]) > 7.0) break; if ((re[3]*re[3] + im[3]*im[3]) > 7.0) break; if ((re[4]*re[4] + im[4]*im[4]) > 7.0) break; #ifdef STRAIGHT horiz_pix = (int) (x_slope * re[0] - x_off); vert_pix = (int) (y_slope * im[0] - y_off); #else horiz_pix = (int) (x_slope * (re[0]-re_position[0]) - x_off); vert_pix = (int) (y_slope * (im[0]-im_position[0]) - y_off); #endif if ( (horiz_pix >= 0) && (horiz_pix < sizex) && (vert_pix >= 0) && (vert_pix < sizey)) { dist = 0.0; for (j=1; j<5; j++) { tmp = im[j] - im[0]; dist += tmp*tmp; tmp = re[j] - re[0]; dist += tmp*tmp; } dist = sqrt (dist); dist = log (dist) - logdelta; dist /= (double) loop; glob [vert_pix*sizex + horiz_pix] += dist; count [vert_pix*sizex + horiz_pix] ++; } } } /* renormalize */ for (i=0; i= isamp) isamp = 1; isamp *= globlen*itermax; irow = isamp / sizey; for (i=0; i 7.0) break; if ((re[1]*re[1] + im[1]*im[1]) > 7.0) break; if ((re[2]*re[2] + im[2]*im[2]) > 7.0) break; if ((re[3]*re[3] + im[3]*im[3]) > 7.0) break; if ((re[4]*re[4] + im[4]*im[4]) > 7.0) break; horiz_pix = (int) (x_slope * re[0] - x_off); vert_pix = (int) (y_slope * im[0] - y_off); if ( (horiz_pix >= 0) && (horiz_pix < sizex) && (vert_pix >= 0) && (vert_pix < sizey)) { for (j=1; j<5; j++) { tmp = (re[j] - re[0]) * (re[j] - re[0]); tmp += (im[j] - im[0]) * (im[j] - im[0]); tmp = 1.0 / sqrt (tmp); vecx [vert_pix*sizex + horiz_pix] += (re[j] - re[0]) *tmp; vecy [vert_pix*sizex + horiz_pix] += (im[j] - im[0]) *tmp; } count [vert_pix*sizex + horiz_pix] ++; } } } /* renormalize */ for (i=0; i= isamp) isamp = 1; isamp *= globlen*itermax; irow = isamp / sizey; for (i=0; i 7.0) break; horiz_pix = (int) (x_slope * re - x_off); vert_pix = (int) (y_slope * im - y_off); if ( (horiz_pix >= 0) && (horiz_pix < sizex) && (vert_pix >= 0) && (vert_pix < sizey)) { glob [vert_pix*sizex + horiz_pix] += last; last = ++ count [vert_pix*sizex + horiz_pix]; } } } /* renormalize */ for (i=0; i= isamp) isamp = 1; isamp *= globlen*itermax; irow = isamp / sizey; for (i=0; i 7.0) break; horiz_pix = (int) (x_slope * re - x_off); vert_pix = (int) (y_slope * im - y_off); if ( (horiz_pix >= 0) && (horiz_pix < sizex) && (vert_pix >= 0) && (vert_pix < sizey) && (2 < loop)) { /* dist = (renew-re)*(renew-re); dist += (imnew-im)*(imnew-im); dist = atan2 (imnew-im, renew-re); */ dist = (renew-re)*(re-reold) + (imnew-im)*(im-imold); norm = (renew-re)*(renew-re) + (imnew-im)*(imnew-im); norm *= (re-reold)*(re-reold) + (im-imold)*(im-imold); if (1.0e-30 < norm) { dist = dist / sqrt (norm); } else { dist = 0.0; } glob [vert_pix*sizex + horiz_pix] += dist; count [vert_pix*sizex + horiz_pix] ++; } reold = re; imold = im; re = renew; im = imnew; } } /* renormalize */ for (i=0; i= isamp) isamp = 1; isamp *= globlen*itermax; irow = isamp / sizey; for (i=0; i 7.0) break; horiz_pix = (int) (x_slope * re - x_off); vert_pix = (int) (y_slope * im - y_off); if ( (horiz_pix >= 0) && (horiz_pix < sizex) && (vert_pix >= 0) && (vert_pix < sizey)) { /* glob [vert_pix*sizex + horiz_pix] += ys; */ glob [vert_pix*sizex + horiz_pix] = xs; count [vert_pix*sizex + horiz_pix] ++; } re = renew; im = imnew; } } /* renormalize */ /* for (i=0; i= isamp) isamp = 1; isamp *= globlen*itermax; irow = isamp / sizey; for (i=0; i 7.0) break; #define OM 10.0 horiz_pix = (int) (x_slope * re - x_off); vert_pix = (int) (y_slope * im - y_off); if ( (horiz_pix >= 0) && (horiz_pix < sizex) && (vert_pix >= 0) && (vert_pix < sizey)) { xs = sin (OM * renew); ys = sin (OM * imnew); glob [vert_pix*sizex + horiz_pix] += xs*xs*ys*ys; count [vert_pix*sizex + horiz_pix] ++; } re = renew; im = imnew; } } /* renormalize */ for (i=0; i= isamp) isamp = 1; isamp *= globlen*itermax; irow = isamp / sizey; for (i=0; i= sizex) hopix = sizex-1; if (vopix >= sizey) vopix = sizey-1; loco = vopix*sizex + hopix; if (i%(irow)==0) printf(" start row %ld\n", i/irow); re = re_position; im = im_position; for (loop=0; loop < LOOP_COUNT; loop++) { tmp = re*re - im*im + re_position; im = 2.0*re*im + im_position; re = tmp; if ((re*re + im*im) > 7.0) break; horiz_pix = (int) (x_slope * re - x_off); vert_pix = (int) (y_slope * im - y_off); if ( (horiz_pix >= 0) && (horiz_pix < sizex) && (vert_pix >= 0) && (vert_pix < sizey)) { ico = count [vert_pix*sizex + horiz_pix] ++; glob [loco] += ico; } } } /* renormalize */ for (i=0; i 7.0) break; } horiz_pix = (int) (x_slope * re_position - x_off); vert_pix = (int) (y_slope * im_position - y_off); glob [vert_pix*sizex + horiz_pix] = ((float) (loop%10)) / 10.0; } } /*-------------------------------------------------------------------*/ extern FILE *Fopen(); int main (int argc, char *argv[]) { float *data; /* my data array */ unsigned int data_width, data_height;/* data array dimensions */ double re_center, im_center, width; int itermax; double renorm, tmp; FILE *fp; int i; char buff [80]; if (5 > argc) { fprintf (stderr, "Usage: %s [ ]\n", argv[0]); exit (1); } itermax = 1; data_width = 200; data_height = 200; data_width = atoi (argv[2]); data_height = atoi (argv[3]); itermax = atoi (argv[4]); data = (float *) malloc (data_width*data_height*sizeof (float)); re_center = -0.6; im_center = 0.0; width = 2.8; if (8 == argc) { re_center = atof (argv[5]); im_center = atof (argv[6]); width = atof (argv[7]); } /* Do the interior now */ renorm = 1.0; if (!strcmp(argv[0], "brat")) mandelbrot_out (data, data_width, data_height, re_center, im_center, width, itermax); if (!strcmp(argv[0], "cutoff")) mandelbrot_cutoff (data, data_width, data_height, re_center, im_center, width, itermax); if (!strcmp(argv[0], "inf")) dmandelbrot_out (data, data_width, data_height, re_center, im_center, width, itermax); if (!strcmp(argv[0], "decide")) mandelbrot_decide (data, data_width, data_height, re_center, im_center, width, itermax); if (!strcmp(argv[0], "manvert")) mandelbrot_inverse (data, data_width, data_height, re_center, im_center, width, itermax); if (!strcmp(argv[0], "measure")) mandelbrot_measure (data, data_width, data_height, re_center, im_center, width, itermax, renorm); if (!strcmp(argv[0], "offset")) mandelbrot_offset (data, data_width, data_height, 0.5, 0.0, 3.0, itermax, renorm); if (!strcmp(argv[0], "mstop")) mandelbrot_stop (data, data_width, data_height, re_center, im_center, width, itermax); if (!strcmp(argv[0], "stalk")) mandelbrot_stalk (data, data_width, data_height, re_center, im_center, width, itermax, 0.0, 0.0); if (!strcmp(argv[0], "orig")) mandelbrot_orig (data, data_width, data_height, re_center, im_center, width, itermax, renorm); if (!strcmp(argv[0], "phase")) mandelbrot_phase (data, data_width, data_height, re_center, im_center, width, itermax, renorm); if (!strcmp(argv[0], "age")) mandelbrot_age (data, data_width, data_height, re_center, im_center, width, itermax, renorm); #ifdef STRAIGHT if (!strcmp(argv[0], "lyapunov")) mandelbrot_lyapunov (data, data_width, data_height, re_center, im_center, width, itermax, renorm); #else if (!strcmp(argv[0], "lyapunov")) mandelbrot_lyapunov (data, data_width, data_height, 0.5, 0.0, 3.0, itermax, renorm); #endif if (!strcmp(argv[0], "migrate")) mandelbrot_migrate (data, data_width, data_height, re_center, im_center, width, itermax, renorm); if (!strcmp(argv[0], "squige")) mandelbrot_squige (data, data_width, data_height, re_center, im_center, width, itermax, renorm); if (!strcmp(argv[0], "next")) mandelbrot_next (data, data_width, data_height, re_center, im_center, width, itermax, renorm); if (!strcmp(argv[0], "whack")) whack (data, data_width, data_height, re_center, im_center, width, itermax, renorm); if (!strcmp(argv[0], "circout")) { re_center = 0.0; im_center = 0.0; width = 1.0; circle_out (data, data_width, data_height, re_center, im_center, width, itermax, 0.33); } if (!strcmp(argv[0], "circin")) { re_center = 0.5; im_center = 0.0; width = 1.0; circle_in (data, data_width, data_height, re_center, im_center, width, itermax, 0.05); } /* ---------------------------------------------------- */ /* make a movie */ if (!strcmp(argv[0], "stalkmov")) { #define NFRAMES 60 for (i=0; i