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

#define MY_MAX_BYTE 100000000

/******************************************************************/
int main( int argc , char * argv[] ){
  int n, i, level, numLevel;
  long double myPi, cos60, sin60;
  long double * xS0, * yS0;
  long double * xS1, * yS1;
  char fileName[100];
  FILE * gp, * sD;
  
  if ( argc != 2 ){
    printf("Use parameter: numLevel\n");
    exit(-9);
  }
  
  numLevel = (int) atoi(argv[1]);
  printf("# numLevel %d\n", numLevel);
  //gp = popen("gnuplot -persist", "w");
  gp = popen("gnuplot", "w");
  if(!gp){
    printf("ABORT failing popening gnuplot pipe\n");
    exit(-9);
  }
  
  myPi = 4. * atan(1);
  cos60 = cos(myPi/3.);
  sin60 = sin(myPi/3.);

  fprintf(gp,"set key bottom right\n");  fflush(gp);
  fprintf(gp,"set size ratio -1\n");  fflush(gp);
  level = -1;
  n = 1;
  // next not needed here but used to be coherent with the following malloc
  if((n+1)*sizeof(long double) > MY_MAX_BYTE){
    printf("ABORT cowardly refusing allocation size %ld max %ld level+1 %d\n",
	   (n+1)*sizeof(long double),(long int)MY_MAX_BYTE,level+1);
    exit(-9);
  }
  xS0 = malloc((n+1) * sizeof(long double));
  if(xS0==NULL){
    printf("ABORT xS0 malloc failed\n");
    exit(-9);
  }
  yS0 = malloc((n+1) * sizeof(long double));
  if(yS0==NULL){
    printf("ABORT yS0 malloc failed\n");
    exit(-9);
  }
  xS0[0] = 0.0;
  xS0[1] = 1.0;
  yS0[0] = 0.0;
  yS0[1] = 1.0;
  sprintf(fileName,"scratchserp.dat%2.2d",level+1);
  printf("%s\n",fileName);
  sD = fopen(fileName,"w");
  for(i=0; i<2*n; i++){
    fprintf(sD, "%Lg %Lg\n", xS0[i], yS0[i]);
  }
  fprintf(gp,"plot [-0.05:1.05][-0.05:1.05] \"%s\" with lines t \" level %d\"\n",fileName, level+1);
  fflush(gp);
  fprintf(gp,"pause 1 \"level %d waiting 1\"\n",level+1); fflush(gp);
  fflush(sD);
  fclose(sD);
  for(level = 0; level < numLevel; level++){
    int oldSegment;
    sprintf(fileName,"scratchserp.dat%2.2d",level+1);
    printf("# %s\n",fileName); fflush(stdout);
    sD = fopen(fileName,"w");
    n = 4 * n;
    printf("****************** level %d n %d B %ld\n",
	   level+1, n, (n+1) * sizeof(long double));
    fflush(stdout);
    if((n+1)*sizeof(long double) > MY_MAX_BYTE){
      printf("ABORT cowardly refusing allocation size %ld max %ld level+1 %d\n",
	     (n+1)*sizeof(long double), (long int)MY_MAX_BYTE,level+1);
      exit(-9);
    }
    xS1 = malloc((n + 1) * sizeof(long double));
    if(xS1==NULL){
      printf("ABORT xS1 malloc failed n = %d\n", n);
      exit(-9);
    }
    yS1 = malloc((n + 1) * sizeof(long double));
    if(yS1==NULL){
      printf("ABORT yS1 malloc failed n = %d\n", n);
      exit(-9);
    }
    xS1[0] = xS0[0];
    yS1[0] = yS0[0];
    for(oldSegment=0; oldSegment<n/4; oldSegment++){
      int a = oldSegment, b = a + 1, c = 1 + 4 * oldSegment;
      long double sXN, sYN;
      // next to minimize addressing and ops
      long double xS0A = xS0[a], xS0B = xS0[b];
      long double yS0A = yS0[a], yS0B = yS0[b];
      long double xD = xS0B - xS0A, yD = yS0B - yS0A;
      long double xD3 = xD/3., yD3 = yD/3.;
      xS1[c] = xS0A + xD3;
      yS1[c] = yS0A + yD3;
      sXN =   cos60 * xD3 - sin60 * yD3;
      sYN = + sin60 * xD3 + cos60 * yD3;
      sXN += xD3 + xS0A;
      sYN += yD3 + yS0A;
      xS1[c+1] = sXN;
      yS1[c+1] = sYN;
      xS1[c+2] = xS0A + 2. * xD3;
      yS1[c+2] = yS0A + 2. * yD3;
      xS1[c+3] = xS0B;
      yS1[c+3] = yS0B;
    }
    free(xS0);
    free(yS0);
    xS0 = xS1;
    yS0 = yS1;
    for(i=0; i<n+1; i++){
      fprintf(sD, "%Lg %Lg\n", xS1[i], yS1[i]);
    }
    fflush(sD);
    fclose(sD);
    fprintf(gp,
	    "plot [-0.05:1.05][-0.05:1.05] \"%s\" with lines t \" level %d\"\n",
	    fileName, level+1);
    fflush(gp);
    fprintf(gp,"pause 2 \"level %d waiting 2\"\n",level+1); fflush(gp);
    fflush(gp);
  }//level
  fprintf(gp,"exit\n");
  fflush(gp);
  return 0;
}