#include <mandelbrot-numerics.h>
#include <mandelbrot-graphics.h>
#include <stdio.h>
#include <string.h> // strncat


/*

cd ~
export LD_LIBRARY_PATH=${HOME}/opt/lib
export PATH=${HOME}/opt/bin:${PATH}
make -C mandelbrot-graphics/c/bin prefix=${HOME}/opt install
m-stretching-cusps 0 0 1 1/2 1/3 0

*/
static const double twopi = 6.283185307179586;




extern int main(int argc, char **argv) {
  // image width and height in pixels
  int side = 2000; // to big ( 10000) gives an error : Segmentation fault (core dumped)
 
  int w = side * 6; // 
  int h = side;
  // 
  double r = 0.5; // proportional to the number of components on the strip, 
  /*
   r = 0.5 gives 4 prominent components counted from period 1 to one side only 
   r = 1.0 gives 10 components
   r = 1.5 gives 15
   r = 2.0 gives 20 ( one can see 2 sides of cardioid ?? because it is near cusp)
   r = 2.5 gives 26
   r = 5.0 gives 50 
  */
  double m = 0.99; // it moves image up and down
  // colors
  m_pixel_t black = m_pixel_rgba(0, 0, 0, 1);
  m_pixel_t white = m_pixel_rgba(1, 1, 1, 1);
  
  double er = 100; //escape radius
  int maxiters = 8192;
  
  
  // check input  
  if (! (argc == 7)) { // usage
  	printf("no input\n");
  	printf("example usage :  \n");
    	printf("%s re(nucleus) im(nucleus) period t_zero t_one t_infinity  \n", argv[0] );
     	printf("%s 0 0 1 1/2 1/3 0\n", argv[0] ); // / (forward slash) is forbidden in Linux filename
    return 1;
  }
  
  // read the input
  complex double nucleus = atof(argv[1]) + I * atof(argv[2]); // center of period componnet 
  int period = atoi(argv[3]);
  // 3 internal angles t of points c defining Moebius transformations 
  mpq_t tzero, tone, tinfinity;
  //
  mpq_init(tzero);
  mpq_init(tone);
  mpq_init(tinfinity);
  //
  mpq_set_str(tzero, argv[4], 10);
  mpq_set_str(tone, argv[5], 10);
  mpq_set_str(tinfinity, argv[6], 10);
  //
  mpq_canonicalize(tzero);
  mpq_canonicalize(tone);
  mpq_canonicalize(tinfinity);

  // name of the file 
  char name [100]; /* name of file */
  int n = snprintf(name, sizeof name, "%d_%f", period, r); /*  */
  if ( n<0 || n>100) 
  	{printf("return value of snprinf : n = %d\n", n);  return 1;}
  char *filename =strncat(name,".png", 4);
  //const char *filename = "06out094.png";	




  m_d_nucleus(&nucleus, nucleus, period, 64);
  // 3 complex points defining Moebius transformation
  complex double zero, one, infinity, z;
  m_d_interior(&z, &zero, nucleus, nucleus, cexp(I * twopi * mpq_get_d(tzero)), period, 64);
  m_d_interior(&z, &one, nucleus, nucleus, cexp(I * twopi * mpq_get_d(tone)), period, 64);
  if (mpq_sgn(tinfinity) == 0) {
    m_d_parent(tinfinity, &infinity, &z, nucleus, period, 64);
  } else {
    m_d_interior(&z, &infinity, nucleus, nucleus, cexp(I * twopi * mpq_get_d(tinfinity)), period, 64);
  }

  printf("parent component with period = %d and nucleus = %.16e %.16e\n", period, creal(nucleus), cimag(nucleus));	
  printf("child component with with internal angle tzero = %s and nucleus c = zero = %.16e %.16e  \n", argv[4], creal(zero), cimag(zero));
  printf("child component with with internal angle tone = %s and nucleus c = one = %.16e %.16e\n", argv[5], creal(one), cimag(one));
  printf("child component with with internal angle tinfinity = %s and nucleus c = infinity = %.16e %.16e\n", argv[6], creal(infinity), cimag(infinity));


  m_d_transform *rect = m_d_transform_rectangular(w, h, I * m * r, r);
  m_d_transform *transform = 0;
  switch (m_d_shape(nucleus, period)) {
    case m_cardioid: {
      complex double size = m_d_size(nucleus, period);
      complex double half, cusp;
      m_d_interior(&z, &half, nucleus, nucleus, -1, period, 64);
      m_d_interior(&z, &cusp, nucleus, nucleus, 1, period, 64);
//      complex double size = cusp - half;
      m_d_transform *cardioid = m_d_transform_cardioid();
      m_d_transform *linear = m_d_transform_linear(-nucleus / size, 1 / size);
      m_d_transform *t0 = m_d_transform_compose(cardioid, linear);
      m_d_transform_reverse(t0, &zero, &z);
      m_d_transform_reverse(t0, &one, &z);
      m_d_transform_reverse(t0, &infinity, &z);
      m_d_transform *moebius = m_d_transform_moebius3(zero, one, infinity);
      printf("Moebius coeeficients a = %f\n", moebius->m.a);
      m_d_transform *t1 = m_d_transform_compose(rect, moebius);
      
      transform = m_d_transform_compose(t1, t0);
      break;
    }
    case m_circle: {
      m_d_transform *moebius = m_d_transform_moebius3(zero, one, infinity);
      transform = m_d_transform_compose(rect, moebius);
      break;
    }
    default: {
      return 1;
    }
  }

  m_d_colour_t *colour = m_d_colour_minimal(white, black, white);
  m_image *image = m_image_new(w, h);
  m_d_render_scanline(image, transform, er, maxiters, colour);
  if (m_image_save_png(image, filename))
  	{printf("image %s saved\n", filename);
  	 printf("filename = period_r\n");}
  	else {printf("image %s not saved\n", filename);}

  return 0;
}