/* mandelbrot
* Copyright (C) 2008 Wikimedia Foundation
*
* 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 2 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, write to the Free Software
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/
#include <stdio.h>
#include <iostream>
#include <cmath>
using namespace std;
/********************************** settings **********************************/
const int w = 4000; // image width
const int h = 1500; // image height
const double mx = 0.7, my = 0.0; // relative position of the origin
const int u = 0.3 * w; // unit
const int iters = 2000; // maximum number of iterations
const bool smooth = true; // flattens the color steps
const int subpix = 16; // must be >= 1; decrease for much faster evaluation!
/******************************************************************************/
double fc(double x)
{ // makes the color slightly smoother
if (x <= 0.0) return(0.0);
if (x >= 1.0) return(1.0);
return(2.0 * (1.0 + x * sqrt(x) - (1 - x) * sqrt(1 - x)) - 3.0 * x);
}
void color(double x, double rgb[])
{
x += 0.52; // arbitrary phase shift
x -= (int)x;
if (x < 0)
x += 1.0;
x *= 6;
int xx = x;
x -= xx;
switch (xx)
{
case 0: rgb[2] = 0; rgb[0] = 1.0; rgb[1] = fc(1.0 - x);
break;
case 1: rgb[1] = 0; rgb[0] = 1.0; rgb[2] = fc(x);
break;
case 2: rgb[1] = 0; rgb[2] = 1.0; rgb[0] = fc(1.0 - x);
break;
case 3: rgb[0] = 0; rgb[2] = 1.0; rgb[1] = fc(x);
break;
case 4: rgb[0] = 0; rgb[1] = 1.0; rgb[2] = fc(1.0 - x);
break;
case 5: rgb[2] = 0; rgb[1] = 1.0; rgb[0] = fc(x);
break;
}
}
double flatten(double x0, double y0, double x, double y)
{ // creates a continuous gradient instead of color steps
int steps = 0;
double sqrxy = x * x + y * y;
if (sqrxy >= 4.0)
{
while (sqrxy < 65536) // high value >> 2
{
double x_old = x;
x = x * x - y * y;
y = 2 * x_old * y;
x += x0;
y += y0;
steps++;
sqrxy = x * x + y * y;
}
double m = steps + 1.5 - log2(log2(sqrxy));
return(m);
}
return(0.0);
}
bool outcircle(double cx, double cy, double r, double x, double y)
{ // checks if (x,y) is outside the circle around (x0,y0) with radius r
x -= cx;
y -= cy;
if (x * x + y * y > r * r)
return(true);
return(false);
}
int main()
{
// open image file
FILE *image_file = fopen("mandelbrot.ppm", "wb");
fprintf(image_file, "P6\n%d %d\n255\n", w, h);
// iterate pixels
for(int row = 0; row < h; row++) {
for(int column = 0; column < w; column++)
{
double sum[3] = { 0.0, 0.0, 0.0 };
// iterate subpixels
for (int i = 0; i < subpix; i++) {
for (int j = 0; j < subpix; j++)
{
// pixel -> coordinates
double x0 = (-w * mx + column
+ (j + 0.5) / subpix) / u;
double y0 = (h * (1.0 - my)
- row - (i + 0.5) / subpix) / u;
double x = x0, y = y0;
int n = -1;
// skip values we know they are inside
if ((outcircle(-0.11, 0.0, 0.63, x0, y0) || x0 > 0.1)
&& outcircle(-1.0, 0.0, 0.25, x0, y0)
&& outcircle(-0.125, 0.744, 0.092, x0, y0)
&& outcircle(-1.308, 0.0, 0.058, x0, y0)
&& outcircle(0.0, 0.25, 0.35, x0, y0))
{
for (int it = 0; it <= iters; it++)
{ // iteration steps
/********* fractal generation ********/
double x_old = x;
x = x * x - y * y;
y = 2 * x_old * y;
x += x0;
y += y0;
/*************************************/
if (x * x + y * y > 4)
{ // value is outside
n = it;
break;
}
}
}
double subpixel[3];
if (n < 0)
{ // black color inside
for (int c = 0; c < 3; c++)
subpixel[c] = 0;
}
else
{
double m = n;
if (smooth)
m += flatten(x0, y0, x, y);
if (m <= -1.8)
m = -0.8;
m = 0.23 * log(1.8 + m);
color(m, subpixel);
}
for (int c = 0; c < 3; c++)
sum[c] += subpixel[c];
}}
// subpixels finished -> make arithmetic mean
char pixel[3];
for (int c = 0; c < 3; c++)
pixel[c] = (int)(255.0 * sum[c] / (subpix * subpix) + 0.5);
fwrite(pixel, 1, 3, image_file);
//pixel finished
}}
fclose(image_file);
return(0);
}