/*
    Code to calculate Mandelbrot fractal.
    https://en.wikipedia.org/wiki/Mandelbrot_set

 Algorithm:

    1)  Let c=x+iy be a coordinate in the plane [-2:2]x[-2:2].
        and initialize z[0]=0

    2)  Iterate z[n+1] = z[n]^2 + c  while |z|<=2 and max iteration not reached.

        sqrt(..) is expensive and not necessary since
        Let z=a+ib then norm(z) = |z| = sqrt( a^2+b^2)
        |z|<=2 corresponds to |z|^2 = a^2+b^2 <= 4
        Hence no need waste time on sqrt's.

    3)  If maximal iteration reached output 0
                                else output iteration/maxiter

 Experiment 1:

   Try changing z[n+1] = z[n]^3 + c  what happens?

   These are called Multibrot fractals

 Experiment 2:

  Change code to:
  Initialize z[0]=x+iy  and c=-0.79 +  0.15i be a constant.

  Run Mandelbrot as before, the resulting set is a Julia set.

  https://en.wikipedia.org/wiki/Julia_set
  https://www.karlsims.com/julia.html

  Matlab to convert to image:
  data=importdata("imagedata.dat");
  imagesc(data);
  axis equal

  Choose e.g. colormap hot
*/


#include <fstream>
#include <iostream>
#include <complex>

using namespace std;

int main()
{
   const int width=1024;    // This is expensive!
   const int height=1024;
   const int itermax=1000;

   const double xmin=-2;
   const double xmax=2;
   const double ymin=-2;
   const double ymax=2;

   ofstream fo("imagedata.dat");

   // Weird order, but Matlab is column-row order
   for (int j=0; j<height; j++)
      {
        for (int i=0; i<width; i++)
           {
              if (j%128==0) cout << "Done " << 100.0*j/height << "%\n";

              // Calculate coordinate.
              double x= (xmax-xmin)*i/width  + xmin;
              double y= (ymax-ymin)*j/height + ymin;

              // Define complex numbers
              complex<double> c( x,y );
              complex<double> z( 0,0 );

              int iter=0;
              while (norm(z)<=4 && iter<itermax)
                 {
                     z=z*z+c;
                     iter++;
                 }

              // choose output
              if (iter==itermax) fo << "0 ";
                            else fo << iter << " ";

           }
         fo << endl;      // shift line after each column.
      }
   fo.close();
}