#include <complex>
#include <fstream>
 
using namespace std;

void SetHSV(double h, double s, double v, unsigned char color[3]) {
  double r=0, g=0, b=0;
  if(s==0)
    r = g = b = v;
  else {
    h /= 60;
    int i = floor(h);
    double f = h - i;
    double p = v*(1-s);
    double q = v*(1-s*f);
    double t = v*(1-s*(1-f));
    
    switch(i){
      case 0: r=v; g=t; b=p; break;
      case 1: r=q; g=v; b=p; break;
      case 2: r=p; g=v; b=t; break;
      case 3: r=p; g=q; b=v; break;
      case 4: r=t; g=p; b=v; break;
      case 5: r=v; g=p; b=q; break;
    }
  }
  int c;
  c = int(256*r); if(c>255) c = 255; color[0] = c;
  c = int(256*g); if(c>255) c = 255; color[1] = c;
  c = int(256*b); if(c>255) c = 255; color[2] = c;
}

complex<double> fun(complex<double>& c ){
  const complex<double> i(0., 1.);
  return (pow(c,2) - 1.) * pow(c -2. -i, 2) / (pow(c, 2) + 2. + 2. * i);
}
 
int main(){
  const int dimx = 800; const int dimy = 800;
  const double rmi = -3; const double rma =  3;
  const double imi = -3; const double ima =  3;

  ofstream f("complex.ppm", ios::binary);
  f << "P6" << endl
    << dimx << " " << dimy << endl
    << "255" << endl;

  for(int j=0; j < dimy; ++j){
    double im = ima - (ima -imi) *j /(dimy -1);
    for(int i=0; i < dimx; ++i){    
      double re = rmi +(rma -rmi) *i /(dimx -1);
      complex<double> c(re, im);
      complex<double> v = fun(c); 
      double a = arg(v)*180/M_PI;
      if(a < 0) a += 360;
      
      double m = abs(v);
      double ranges;
      double rangee = 1;

      while(m>rangee){
          ranges = rangee;
          rangee *= M_E;
      }

      double k   = (m-ranges)/(rangee-ranges);
      double sat = k < 0.5 ? 1-2.8*k : 2.8*k-1.1;
      if(sat < 0.3) sat = 0.3;
      else if(sat > 1) sat = 1;

      double val = k < 0.5 ? 1.4-1.6*k : 1.6*k-0.6;
      if(val < 0.6) val = 0.6;
      else if(val > 1) val = 1;

      unsigned char color[3];
      SetHSV(a,sat,val,color);
      f.write(color,3);
    }
  }
  return 0;
}