#include <iostream>
#include <math.h>
#include <stdlib.h>
#include <fstream>
#include <map>
using namespace std;

struct Colour
{
  int r, g, b;
  Colour();
  Colour(int red, int green, int blue);
};
Colour::Colour() { r = 0; g = 0; b = 0; }
Colour::Colour(int red, int green, int blue)
{
  r = red;
  if (r > 255) r = 255;  if (r < 0) r = 0;
  g = green;
  if (g > 255) g = 255;  if (g < 0) g = 0;
  b = blue;
  if (b > 255) b = 255;  if (b < 0) b = 0;
}
class ColourMap
{
  private:
    map<double,Colour> mapData;
  public:
    ColourMap();
    bool addEntry(double mapIndex, Colour mapColour);
    Colour getColour(double mapIndex);
};
ColourMap::ColourMap()
{
  mapData[0] = Colour(0, 0, 0);
  mapData[1] = Colour(255, 255, 255);
}
bool ColourMap::addEntry(double mapIndex, Colour mapColour)
{
  bool result = false;
  if (mapIndex >= 0 && mapIndex <= 1) {
    mapData[mapIndex] = mapColour;
	result = true; }
  return result;
}
Colour ColourMap::getColour(double mapIndex)
{
  Colour result = Colour(0, 0, 0);
  if (mapIndex == 0) { result = mapData.begin()->second; }
  if (mapIndex > 0 && mapIndex <= 1)
  {
    map<double,Colour>::iterator p = mapData.begin();
	while (p->first < mapIndex) p++;
	double index2 = p->first;
	Colour mapColour2 = p->second;
	p--;
	double index1 = p->first;
	Colour mapColour1 = p->second;
	double midIndex = (mapIndex-index1)/(index2-index1);
	result.r = int(mapColour1.r + (mapColour2.r-mapColour1.r)*midIndex);
	result.g = int(mapColour1.g + (mapColour2.g-mapColour1.g)*midIndex);
	result.b = int(mapColour1.b + (mapColour2.b-mapColour1.b)*midIndex);
  }
  return result;
}

const int mn = 5;
int t;
double x, y, z, newx, newy, newz;
double at[18];

void calculatePoint();
void error(const char* p);

void error(const char* p)
{
  cerr << p << '\n';
  exit(1);
}

void calculatePoint()
{
  switch(t)
  {
    case 0:
      newx = x + at[0]*at[3]*(y - x);
      newy = y + at[3]*(at[1]*x - y - z*x);
      newz = z + at[3]*(x*y - at[2]*z);
      break;
    case 1:
      newx = x + at[4]*(-at[0]*x - y*y - z*z + at[0]*at[2]);
      newy = y + at[4]*(-y + x*y - at[1]*x*z + at[3]);
      newz = z + at[4]*(-z + at[1]*x*y + x*z);
      break;
    case 2:
      newx = sin(at[0]*y) - z*cos(at[1]*x);
      newy = z*sin(at[2]*x) - cos(at[3]*y);
      newz = sin(x);
      break;
    case 3:
      newx = at[0] + y - z*y;
      newy = at[1] + z - x*z;
      newz = at[2] + x - y*x;
      break;
    case 4:
      newx = at[0] + y - z*(at[1] + y);
      newy = at[2] + z - x*(at[3] + z);
      newz = at[4] + x - y*(at[5] + x);
      break;
    case 5:
      newx = at[0] + x*(at[1]+at[2]*x+at[3]*y) + y*(at[4]+at[5]*y);
      newy = at[6] + y*(at[7]+at[8]*y+at[9]*z) + z*(at[10]+at[11]*z);
      newz = at[12] + z*(at[13]+at[14]*z+at[15]*x) + x*(at[16]+at[17]*x);
      break;
    default:
	  break;
  }
}

int main (int argc, char * const argv[]) {
  int imageSize, pn, iterations, mapno, brightfunc;
  unsigned** image;
  unsigned** zimage;
  double xmax, xmin, ymax, ymin, zmax, zmin;
  // set defaults
  imageSize = 800;
  iterations = int(1e6);
  mapno = 1;
  brightfunc = 0;
  ColourMap cmp[mn];
  cmp[0].addEntry(0, Colour(255,255,255));
  cmp[0].addEntry(1, Colour(255,255,255));
  cmp[1].addEntry(0, Colour(255,128,128));
  cmp[1].addEntry(0.5, Colour(128,255,128));
  cmp[1].addEntry(1, Colour(128,128,255));
  cmp[2].addEntry(0, Colour(255,64,0));
  cmp[2].addEntry(0.5, Colour(255,255,128));
  cmp[2].addEntry(1, Colour(255,255,255));
  cmp[3].addEntry(0, Colour(255,255,64));
  cmp[3].addEntry(0.8, Colour(0,128,255));
  cmp[3].addEntry(1, Colour(255,255,255));
  cmp[4].addEntry(0, Colour(192,64,255));
  cmp[4].addEntry(0.5, Colour(255,255,64));
  cmp[4].addEntry(1, Colour(255,128,128));
  // read arguments
  if (argc < 2)
	error("not enough arguments");
  const char* infile = argv[1];
  string s = infile;
  if (s.find(".") < s.size())
	s = s.substr(0,s.find("."));
  s += ".tga";
  const char* outfile = s.data();
  for (int i = 2; i < argc; i++)
  {
    string argstr = argv[i];
	if (argstr[0] == '-')
	{
	  char opt = argstr[1];
	  string arg = argstr.substr(2,s.size()-2);
	  switch (opt)
	  {
	  case 's':
        imageSize = atoi(arg.c_str());
        if (imageSize < 32 || imageSize > 4000) {
		  cout << "Warning: image size outside normal bounds; using default.";
          imageSize = 800; }
		break;
	  case 'i':
	    iterations = int(atof(arg.c_str()));
		if (iterations < 0) {
		  cout << "Warning: negative iterations; using default.";
		  iterations = int(1e6); }
		if (iterations > int(1e10)) cout << "Warning: iteration count is high.";
		break;
	  case 'c':
	    mapno = atoi(arg.c_str());
		if (mapno > mn-1) mapno = mn-1;
		if (mapno < 0) mapno = 0;
		break;
	  case 'b':
	    brightfunc = atoi(arg.c_str());
		if (brightfunc > 2) brightfunc = 2;
		if (brightfunc < 0) brightfunc = 0;
		break;
	  default:
	    break;
	  }
	}
  }
  ifstream in(infile);
  in >> t;
  in >> pn;
  for (int i = 0; i < pn; i++)
	in >> at[i];

  // set up
  x = 0;  y = 0;  z = 0;
  if (t == 0 || t == 1)
  {
    x = (rand()/RAND_MAX)*0.1 - 0.05;
	y = (rand()/RAND_MAX)*0.1 - 0.05;
	z = (rand()/RAND_MAX)*0.1 - 0.05;
  }
  image = new unsigned*[imageSize];
  zimage = new unsigned*[imageSize];
  for (int i = 0; i < imageSize; i++) {
    image[i] = new unsigned[imageSize];
	zimage[i] = new unsigned[imageSize]; }
  for (int i = 0; i < imageSize; i++)
    for (int j = 0; j < imageSize; j++) {
	  image[i][j] = 0;
	  zimage[i][j] = 0; }

  // iterate test to find centre, scaling etc
  xmax = 0;  ymax = 0;  zmax = 0;
  xmin = 0;  ymin = 0;  zmin = 0;
  for (int iter = 0; iter < 10000; iter++)
  {
    calculatePoint();
    if (newx > xmax) xmax = newx;
    if (newx < xmin) xmin = newx;
    if (newy > ymax) ymax = newy;
    if (newy < ymin) ymin = newy;
	if (newz > zmax) zmax = newz;
	if (newz < zmin) zmin = newz;
    x = newx;  y = newy;  z = newz;
  }
  double xwidth = xmax - xmin;
  double ywidth = ymax - ymin;
  double zwidth = zmax - zmin;
  double maxWidth = max(xwidth, max(ywidth, zwidth));
  double scale = 0.95*(imageSize/maxWidth);
  double cx = xmin + xwidth/2;
  double cy = ymin + ywidth/2;
  double cz = zmin + zwidth/2;

  // iterate final image
  for (int iter = 0; iter < iterations; iter++)
  {
    calculatePoint();
	int xind = floor((newx-cx)*scale) + imageSize/2;
	int yind = floor((newy-cy)*scale) + imageSize/2;
	int zind = floor((newz-cz)*scale) + imageSize/2;
	if (yind > zimage[xind][zind])
      zimage[xind][zind] = yind;
    image[xind][zind] += 1;
	x = newx;
	y = newy;
	z = newz;
  }

  // calculate max brightness
  int pmax = 0;
  for (int i = 0; i < imageSize; i++)
    for (int j = 0; j < imageSize; j++)
	  if (image[i][j] > pmax) pmax = image[i][j];

  // save .tga file
  ofstream tga(outfile, ios_base::binary);
  // no of chars in id field
  tga.put(char(0));
  // color map type
  tga.put(char(0));
  // image type code
  tga.put(char(2));
  // color map spec; leave blank
  for (int i = 0; i < 5; i++)
	tga.put(char(0));
  // image spec
  // x origin
  tga.put(char(0));  tga.put(char(0));
  // y origin
  tga.put(char(0));  tga.put(char(0));
  // width of image
  unsigned short si = imageSize;
  tga.put(char(si));
  tga.put(char(si>>8));
  // height of image
  tga.put(char(si));
  tga.put(char(si>>8));
  // image pixel size
  tga.put(char(24));
  // image descriptor byte
  tga.put(char(0));
  // image id field - absent
  // color map data - absent
  // image data
  double brightness, xval;
  char* p = new char[3];
  Colour c;
  for (int i = 0; i < imageSize; i++)
    for (int j = 0; j < imageSize; j++)
	{
	  xval = double(image[i][j])/pmax;
	  switch (brightfunc)
	  {
	    case 0:
          brightness = pow(xval, 0.25);
		  break;
		case 1:
		  brightness = 1-exp(-xval*32);
		  break;
		case 2:
		  brightness = sqrt(2*xval-xval*xval);
		  break;
	  }
	  c = cmp[mapno].getColour(double(zimage[i][j])/imageSize);
      p[0] = (char)int(c.b*brightness);
	  p[1] = (char)int(c.g*brightness);
	  p[2] = (char)int(c.r*brightness);
	  tga.write(p, 3);
    }
  delete [] p;
  delete [] image;
  delete [] zimage;
  return 0;
}