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


#define SAMPLESAMOUNT       (32000)
#define PI                  (3.141592653589793238)

typedef struct sector SECTOR;
struct sector {
	char   type;
	// geeft aan of dit een blad of een node is
	// type = 0  ==> leeg (dus nog geen punt of subsectors)
	// type = 1  ==> blad (bevat punt)
	// type = 2  ==> node (bevat subsectors)
	SECTOR* child1;
	SECTOR* child2;
	SECTOR* child3;
	SECTOR* child4;
	double  x;
	double  y;
	double  xmin;
	double  xmax;
	double  ymin;
	double  ymax;
};

SECTOR* SECTOR_createnew(double x0, double x1, double y0, double y1);
void    SECTOR_destruct(SECTOR* sect);
void    SECTOR_addpoint(SECTOR* sect, double x, double y);
void    SECTOR_showall(SECTOR* sect);
void    SECTOR_findnearestpoint(SECTOR* sect, double* x, double* y, double* bestfound);


int main(void) {
  SECTOR*   sect           = SECTOR_createnew(-1.0001, 1.0001, -1.0001, 1.0001);
  FILE*     filehandle1;
  FILE*     filehandle2;
  FILE*     filehandle3;
  double*   x              = (double*) malloc(SAMPLESAMOUNT * sizeof(double));
  double*   y              = (double*) malloc(SAMPLESAMOUNT * sizeof(double));
  int       i              = 0;
  double	goldenratio    = .5*(1+sqrt(5));
  double*   goldenratio_x  = (double*) malloc(12 * SAMPLESAMOUNT * sizeof(double));
  double*   goldenratio_y  = (double*) malloc(12 * SAMPLESAMOUNT * sizeof(double));
  double*   outmindist     = (double*) malloc(SAMPLESAMOUNT * sizeof(double));


  /******* CREATE GOLDENRATIO SET *******/
  for (i = 0; i < SAMPLESAMOUNT * 12; i++) {
	printf("\rINITIALIZING...      : %d of %d", i, SAMPLESAMOUNT*12);
    double phi   = 2*PI*goldenratio*(double)i;
	double sqrtr = sqrt((double)i)/sqrt((double)SAMPLESAMOUNT*12.);
    goldenratio_x[i] = cos(phi)*sqrtr;
    goldenratio_y[i] = sin(phi)*sqrtr;
  }


  /******* CREATE DISTRIBUTION *******/

  x[0] = 0.0;
  y[0] = 0.0;
  SECTOR_addpoint(sect, x[0], y[0]);

  for (i = 1; i < SAMPLESAMOUNT; i++) {
    int    i2;
	int    groffset = 0;
    double olddist  = 0.0;

    printf("\rSETUP DISTRIBUTION   : %d of %d   ", i, SAMPLESAMOUNT);
	if (!(i % 1000)) printf("\007");

	for (i2 = SAMPLESAMOUNT*12 - i; i2 >= 0; i2--) {
	  double mindist = 5.;
	  SECTOR_findnearestpoint(sect, goldenratio_x + i2, goldenratio_y + i2, &mindist);
	  if (mindist > olddist) {
		groffset = i2;
		olddist = mindist;
		outmindist[i] = mindist;
		x[i] = goldenratio_x[i2];
		y[i] = goldenratio_y[i2];
	  }
	}

	goldenratio_x[groffset] = goldenratio_x[SAMPLESAMOUNT*12-i];
	goldenratio_y[groffset] = goldenratio_y[SAMPLESAMOUNT*12-i];

    SECTOR_addpoint(sect, x[i], y[i]);
  }
  printf("\rsetup distribution   : done                 \r\n");


  /******* WRITE SAMPLES TO DISK *******/

  printf("SAVE SAMPLES TO DISK...");
  filehandle1 = fopen("newraddata.inc", "w");
  filehandle2 = fopen("newraddata.txt", "w");
  filehandle3 = fopen("newraddata_mindist.txt", "w");

  fprintf(filehandle1, "#declare ThingAnimated = union {\r\n");
  for (i = 0; i < SAMPLESAMOUNT; i++) {
    double z = 1.0-x[i]*x[i]-y[i]*y[i];
	if (z >= 0.) {
	  z = sqrt(z);
	} else {
	  z = 0.;
	}
	fprintf(filehandle1, "  #if (FRAMENUMBER > %d)\r\n    sphere {<%2.12f, %2.12f, %2.12f>, PointSize}\r\n  #end\r\n", i, x[i], y[i], z);
	fprintf(filehandle3, "#if (FRAMENUMBER > %d)\r\n  #declare MINDIST = %2.12f\r\n#end\r\n", outmindist[i]);
    fprintf(filehandle2, "\t{%2.12f, %2.12f, %2.12f}", x[i], y[i], z);
	if (i < SAMPLESAMOUNT-1) fprintf(filehandle2, ",\n");
  }
  fprintf(filehandle1, "}\r\n");

  fclose(filehandle1);
  fclose(filehandle2);
  fclose(filehandle3);
  printf("\rsave samples to disk : done                \r\n\007\007\007");



  /******* FREE ALLOCATED MEMORY *******/
  free(x);
  free(y);
  free(goldenratio_x);
  free(goldenratio_y);
  free(outmindist);
  SECTOR_destruct(sect);
  return (0);
}
























SECTOR* SECTOR_createnew(double x0, double x1, double y0, double y1) {
  SECTOR* newsect = (SECTOR*) malloc(sizeof(SECTOR));

  if (newsect) {
    newsect->child1 = 0;
    newsect->child2 = 0;
    newsect->child3 = 0;
    newsect->child4 = 0;
	newsect->xmin   = x0;
	newsect->xmax   = x1;
	newsect->ymin   = y0;
	newsect->ymax   = y1;
    newsect->type   = 0;
  }
  return newsect;
}

void SECTOR_destruct(SECTOR* sect) {
  if (sect) {
    SECTOR_destruct(sect->child1);
	SECTOR_destruct(sect->child2);
	SECTOR_destruct(sect->child3);
	SECTOR_destruct(sect->child4);
	free(sect);
  }
}

void SECTOR_addpoint(SECTOR* sect, double x, double y) {
  if (sect) {	 
	if (sect->type == 0) {
	  sect->x = x;
	  sect->y = y;
	} else if (sect->type == 1) {
	  double middle_x = (sect->xmax + sect->xmin) / 2.;
	  double middle_y = (sect->ymax + sect->ymin) / 2.;

	  sect->child1 = SECTOR_createnew(sect->xmin, middle_x, sect->ymin, middle_y);
	  sect->child2 = SECTOR_createnew(middle_x, sect->xmax, sect->ymin, middle_y);
	  sect->child3 = SECTOR_createnew(sect->xmin, middle_x, middle_y, sect->ymax);
	  sect->child4 = SECTOR_createnew(middle_x, sect->xmax, middle_y, sect->ymax);

	  //reeds bestaande punt in een nieuwe subsector plaatsen
	  if (sect->x < middle_x) {
		if (sect->y < middle_y) {
		  SECTOR_addpoint(sect->child1, sect->x, sect->y);
		} else {
		  SECTOR_addpoint(sect->child3, sect->x, sect->y);
		}
	  } else {
		if (sect->y < middle_y) {
		  SECTOR_addpoint(sect->child2, sect->x, sect->y);
		} else {
		  SECTOR_addpoint(sect->child4, sect->x, sect->y);
		}
	  }

	  //nieuwe punt in de juiste subsector plaatsen.
	  if (x < middle_x) {
		if (y < middle_y) {
		  SECTOR_addpoint(sect->child1, x, y);
		} else {
		  SECTOR_addpoint(sect->child3, x, y);
		}
	  } else {
		if (y < middle_y) {
		  SECTOR_addpoint(sect->child2, x, y);
		} else {
		  SECTOR_addpoint(sect->child4, x, y);
		}
	  }
	} else if (sect->type == 2) {
	  double middle_x = (sect->xmax + sect->xmin) / 2.;
	  double middle_y = (sect->ymax + sect->ymin) / 2.;

	  //nieuwe punt in de juiste subsector plaatsen.
	  if (x < middle_x) {
		if (y < middle_y) {
		  SECTOR_addpoint(sect->child1, x, y);
		} else {
		  SECTOR_addpoint(sect->child3, x, y);
		}
	  } else {
		if (y < middle_y) {
		  SECTOR_addpoint(sect->child2, x, y);
		} else {
		  SECTOR_addpoint(sect->child4, x, y);
		}
	  }
	}
	
	if (sect->type < 2) sect->type = sect->type + 1;
  }
}

void SECTOR_showall(SECTOR* sect) {
  if (sect && (sect->type)) {
	if (sect->type == 1) {
	  //sector bevat slechts een punt, dus geen subsectoren
	  printf("<%1.10f, %1.10f>\r\n", sect->x, sect->y);
	} else {
      //sector bevat subsectoren, dus allemaal afgaan.
	  SECTOR_showall(sect->child1);
	  SECTOR_showall(sect->child2);
	  SECTOR_showall(sect->child3);
	  SECTOR_showall(sect->child4);
	}
  }
}

void SECTOR_findnearestpoint(SECTOR* sect, double* x, double* y, double* bestfound) {
  if (sect && sect->type) {
	if (sect->type == 1) {
	  //de sector bevat slechts een punt en geen subs, dus dit punt is het dichtstebij
	  double d = ((sect->x-*x)*(sect->x-*x)+(sect->y-*y)*(sect->y-*y));
	  if (d < *bestfound) *bestfound = d;
	} else {
	  //de sector bevat subsectoren
      if (*x < sect->xmin) {
		if (*y < sect->ymin) {
		  // X23
		  // 456
		  // 789
		  double d = (sect->xmin-*x)*(sect->xmin-*x)+(sect->ymin-*y)*(sect->ymin-*y);
		  if (d < *bestfound) {
			SECTOR_findnearestpoint(sect->child1, x, y, bestfound);
			SECTOR_findnearestpoint(sect->child2, x, y, bestfound);
			SECTOR_findnearestpoint(sect->child3, x, y, bestfound);
		  }
		} else if (*y >= sect->ymin && *y < sect->ymax) {
		  // 123
		  // X56
		  // 789
		  double d = (sect->xmin-*x)*(sect->xmin-*x);
		  if (d < *bestfound) {
		    SECTOR_findnearestpoint(sect->child1, x, y, bestfound);
		    SECTOR_findnearestpoint(sect->child3, x, y, bestfound);
		    SECTOR_findnearestpoint(sect->child2, x, y, bestfound);
		    SECTOR_findnearestpoint(sect->child4, x, y, bestfound);
		  }
	    } else if (*y >= sect->ymax) {
		  // 123
		  // 234
		  // X67
		  double d = (sect->xmin-*x)*(sect->xmin-*x)+(sect->ymax-*y)*(sect->ymax-*y);
		  if (d < *bestfound) {
			SECTOR_findnearestpoint(sect->child3, x, y, bestfound);
			SECTOR_findnearestpoint(sect->child1, x, y, bestfound);
			SECTOR_findnearestpoint(sect->child4, x, y, bestfound);
		  }
		}
	  } else if (*x >= sect->xmin && *x < sect->xmax) {
		if (*y < sect->ymin) {
		  // 1X3
		  // 456
		  // 789
		  double d = (sect->ymin-*y)*(sect->ymin-*y);
		  if (d < *bestfound) {
		    SECTOR_findnearestpoint(sect->child1, x, y, bestfound);
		    SECTOR_findnearestpoint(sect->child2, x, y, bestfound);
		    SECTOR_findnearestpoint(sect->child3, x, y, bestfound);
		    SECTOR_findnearestpoint(sect->child4, x, y, bestfound);
		  }
		} else if (*y >= sect->ymin && *y < sect->ymax) {
		  // 123
		  // 4X6
		  // 789

          //ingeval het punt binnen de sectorgrenzen ligt, check alle vier de subsectoren
          SECTOR_findnearestpoint(sect->child1, x, y, bestfound);
          SECTOR_findnearestpoint(sect->child2, x, y, bestfound);
          SECTOR_findnearestpoint(sect->child3, x, y, bestfound);
          SECTOR_findnearestpoint(sect->child4, x, y, bestfound);

		} else if (*y >= sect->ymax) {
		  // 123
		  // 456
		  // 7X9
		  double d = (sect->ymax-*y)*(sect->ymax-*y);
		  if (d < *bestfound) {
		    SECTOR_findnearestpoint(sect->child3, x, y, bestfound);
		    SECTOR_findnearestpoint(sect->child4, x, y, bestfound);
		    SECTOR_findnearestpoint(sect->child1, x, y, bestfound);
		    SECTOR_findnearestpoint(sect->child2, x, y, bestfound);
		  }
		}
	  } else if (*x >= sect->xmax) {
		if (*y < sect->ymin) {
		  // 12X
		  // 456
		  // 789
		  double d = (sect->xmax-*x)*(sect->xmax-*x)+(sect->ymin-*y)*(sect->ymin-*y);
		  if (d < *bestfound) {
			SECTOR_findnearestpoint(sect->child2, x, y, bestfound);
			SECTOR_findnearestpoint(sect->child1, x, y, bestfound);
			SECTOR_findnearestpoint(sect->child4, x, y, bestfound);
		  }
		} else if (*y >= sect->ymin && *y < sect->ymax) {
		  // 123
		  // 45X
		  // 789
		  double d = (sect->xmax-*x)*(sect->xmax-*x);
		  if (d < *bestfound) {
		    SECTOR_findnearestpoint(sect->child2, x, y, bestfound);
		    SECTOR_findnearestpoint(sect->child4, x, y, bestfound);
		    SECTOR_findnearestpoint(sect->child1, x, y, bestfound);
		    SECTOR_findnearestpoint(sect->child3, x, y, bestfound);
		  }
		} else if (*y >= sect->ymax) {
		  // 123
		  // 456
		  // 78X
		  double d = (sect->xmax-*x)*(sect->xmax-*x)+(sect->ymax-*y)*(sect->ymax-*y);
		  if (d < *bestfound) {
			SECTOR_findnearestpoint(sect->child4, x, y, bestfound);
			SECTOR_findnearestpoint(sect->child2, x, y, bestfound);
			SECTOR_findnearestpoint(sect->child3, x, y, bestfound);
		  }
		}
	  }
	}
  }
}