// A thing to generate a Feigenbaum diagram.
// by Geoff Richards.


#include <iostream>
#include <cmath>
#include <cassert>


const double X1 = 1, X2 = 4, Y1 = 0, Y2 = 1;
const int N_ITER = 600, N_HIDE = 300;
const int XS = 1024, YS = 800;

const int SHADES = 16;
const char *GREY = "XZYxzyt%$&+-;,. ";
const char *HEX = "0123456789ABCDEF";

char **img;


void
plot (int x, double y)
{
   assert (x >= 0);
   assert (y >= 0);
   assert (x < XS);
   assert (y < YS);

   int yi = int (y);
   double dy = y - yi;

   int level = int (fabs (0.5 - dy) * (SHADES * 2.0));
   assert (level >= 0 && level <= SHADES);
   img[x][YS - yi] += level;
   if (img[x][YS - yi] >= SHADES)
      img[x][YS - yi] = SHADES - 1;

}


void
write_xpm ()
{
   cout << "/* XPM */\n"
           "static char * feigenbaum_xpm[] = {\n"
           "\"" << XS << ' ' << YS << " " << SHADES << " 1\"\n";
   for (int c = 0; c < SHADES; ++c)
   {
      cout << '\"' << GREY[c] << "\tc #";
      int col = 255 - int (c / (SHADES - 1.0) * 255);
      for (int i = 0; i < 3; ++i)
         cout << HEX[col / 0x10] << HEX[col & 0xF];
      cout << "\"\n";
   }

   for (int y = 0; y < YS; ++y)
   {
      cout << '\"';

      for (int x = 0; x < XS; ++x)
      {
         assert (img[x][y] >= 0 && img[x][y] < 16);
         cout << GREY[img[x][y]];
      }

      if (y < YS - 1)
         cout << "\",\n";
      else
         cout << "\"};\n";
   }
}


int
main ()
{
   img = new char*[XS];
   for (int x = 0; x < XS; ++x)
   {
      img[x] = new char[YS];
      for (int y = 0; y < YS; ++y)
         img[x][y] = 0;
   }

   for (int xx = 0; xx < XS * 4; ++xx)
   {
      double x = xx / 4.0;
      const double a = X1 + (X2 - X1) / XS * x;
      double v = 0.5;

      for (int i = 0; i < N_ITER; ++i)
      {
         v = a * v * (1 - v);
         const double y = (v - Y1) / ((Y2 - Y1) / YS);

         if (i > N_HIDE && x >= 0 && x < XS && y >= 0 && y < YS)
            plot (int (x), y);
      }
   }

   write_xpm();
}