#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

#define XMIN -1.4
#define YMIN -1.1
#define ZMIN -1.4
#define XMAX 1.4
#define YMAX 1.1
#define ZMAX 1.4
#define XGRAIN .02
#define YGRAIN .02
#define ZGRAIN .01

// globals, yeah
//double c[]={-0.213,-0.0410,-0.563,-0.560};
double c[]={-0.125,-0.256,0.847,0.0895};
//double c[]={0.7, 0.5, 0.2, 0.1};

double clast=0.; // this is the 5th loose variable: converting from 3d point to 4d quaternion.

double eye[]={0., 0., -3.};
double light[]={2., -2., 3.}; //-2, 3, -3

int maxcolor=511;
int maxiter=10;

double zarr[(int)((XMAX-XMIN)/XGRAIN+2)][3];
int zarri;
double max[]={XMAX, YMAX, ZMAX};
double min[]={XMIN, YMIN, ZMIN};
double grain[]={XGRAIN, YGRAIN, ZGRAIN};

// bigz// double biggestz=-10.0;

// fn prototypes
void colorrow(double y); //, double zarr[][3], int zarri
double dot(double *a, double *b);
void cross(double *a, double *b, double *product);
void processrow(double viewy);
double processvector(double x, double y);
int jpoint(double *zin);
double magnitude(double *vec, int elts);
double fake_mag(double *vec, int elts);
void mkunitvec(double *vec, int elts);


double stof(char *s) {
  

  return 0;
}

// here be dragans
int main(int argc, char *argv[]) {

  if (argc>1) {
    char *end;
    clast=strtod(argv[1], &end);
  }
  

  //PGM heading // IGNORED due to #
  printf("P2 %.0f %.0f %d\n", (((max[0]-min[0])/grain[0])) - 1.0, (((max[1]-min[1])/grain[1])) - 1.0, maxcolor);

  printf("#c: %f %f %f %f; clast=%f; light %f %f %f; eye %f %f %f; min %f %f %f; max %f %f %f; grain %f %f %f\n", c[0],c[1],c[2],c[3],clast,light[0],light[1],light[2],eye[0],eye[1],eye[2],min[0],min[1],min[2],max[0],max[1],max[1], grain[0],grain[1],grain[2]);

  double viewy=max[1];
  
  zarri=0;
  processrow(viewy);
  zarri++;  

  viewy -= grain[1];
  processrow(viewy);
  
  zarri++;
  while (viewy>min[1]) {
    viewy -=grain[1];
    processrow(viewy);
    //printf("#going to color: %f\n", zarr[0][0]);
    colorrow(viewy+grain[1]);
    if (zarri==2)
      zarri=0;
    else zarri++;
  }
  //printf("#biggestz: %f\n",biggestz);
  return 0;
}

void colorrow(double y) {
  //printf("#colorrow(%f)\n",y);
  int zi, zabove, zbelow;
  double x=min[0];
  int i=0;
  int col;

  double avec[3], bvec[3], nvec[3], lightvec[3];

  if (zarri==0) zi=2;
  else zi=zarri-1;

  if (zi==0) {zabove=2; zbelow=1;} 
  else if (zi==1) {zabove=0; zbelow=2;}
  else {zabove=1; zbelow=0;}
  while (x < max[0]-grain[0]) {

    if (zarr[i][zi] == max[2]) {
      printf("%d ", maxcolor);
    }
    else {
      if (x==min[0]){
	avec[0]=0-grain[0];
	avec[1]=0;
	avec[2]=zarr[i][zi] - zarr[i+1][zi];

	bvec[0]=0;
	bvec[1]=2*grain[1];
	bvec[2]=zarr[i][zabove] - zarr[i][zbelow];
      }
      else {
	avec[0]=2*grain[0];
	avec[1]=0;
	avec[2]=zarr[i+1][zi] - zarr[i-1][zi];

	bvec[0]=0;
	bvec[1]= 2*grain[1];
	bvec[2]=zarr[i][zabove] - zarr[i][zbelow];
      }

      lightvec[0]=x-light[0];
      lightvec[1]=y-light[1];
      lightvec[2]=zarr[i][zi] - light[2];
      
      //nvec=cross(avec, bvec)
      cross(avec, bvec, nvec);

      mkunitvec(lightvec, 3);
      mkunitvec(nvec, 3);
      
      col = (int)(maxcolor * (dot(lightvec, nvec)));
      printf("%d ", abs(col));
    }
    x+=grain[0];
    i++;
  }
  printf("\n");
}


double dot(double *a, double *b) {
  return (a[0]*b[0]+a[1]*b[1]+a[2]*b[2]);
}

void cross(double *a, double *b, double *product) {
  product[0] = a[1]*b[2] - a[2]*b[1];
  product[1] = a[2]*b[0] - a[0]*b[2];
  product[2] = a[0]*b[1] - a[1]*b[0];
}

void processrow(double viewy) {
  //printf("#processrow(%f), zarri: %d\n",viewy, zarri);
  double viewx=min[0];
  int pix=0;
  //double thisz;
  while (viewx < max[0]) {
    //thisz= processvector(viewx, viewy);
    zarr[pix][zarri] = processvector(viewx, viewy);;
    //if ((thisz > biggestz) && thisz!=max[2]) {biggestz=thisz;}
    viewx += grain[0];
    pix++;
  }
  //printf("#processrow leaving(%f), zarri: %d\n",viewy, zarri);  
}
  

double processvector(double x, double y) {
  //printf("#processvector(%f, %f)\n", x, y);
  double point[4]={x, y, min[2], clast};
  int j=0;
  while (point[2] < max[2]) {
    //printf("jpoint at %f %f %f\n",x,y,point[2]);
    j=jpoint(point);
    if (j==1) {
      //printf("#processvector exitwounds a: %f\n",point[2]);
      return point[2];
    }
    point[2] += grain[2];
  }
  //printf("#processvector exitwounds b : %f\n", max[2]);
  return max[2];
}

int jpoint(double *zin) {
  int i=0,j;
  double znew[4];
  double z[4]={zin[0], zin[1], zin[2], zin[3]};
  double ztwo, mag;
  while (i<maxiter) {
    
    i++;
    ztwo=2*z[0];
    znew[0]=( z[0]*z[0] - z[1]*z[1] - z[2]*z[2] - z[3]*z[3] + c[0] );
    znew[1]=( ztwo*z[1] + c[1] );
    znew[2]=( ztwo*z[2] + c[2] );
    znew[3]=( ztwo*z[3] + c[3] );

    for (j=0;j<4;j++) {z[j]=znew[j];}
    mag=fake_mag(z, 4);
    if (mag > 8.) {
      //printf("#jpoint leaving\n");
      return 0;
    }
  }
  return 1;
}

double magnitude(double *vec, int elts) {
  double sum=0;
  int i;
  for (i=0; i<elts; i++) {
    sum+=pow(vec[i], 2);
  }
  return sqrt(sum);
}

double fake_mag(double *vec, int elts) {
  double sum=0;
  int i;
  for (i=0; i<elts; i++) {
    sum+=pow(vec[i], 2);
  }
  return sum;
}

void mkunitvec(double *vec, int elts) {
  double d=magnitude(vec, elts);
  int i;
  for(i=0; i<elts; i++) {
    vec[i]/=d;
  }
}