/*

gcc d.c -Wall -lm
./a.out

program creates 2 files 
- svg file ( without labels which are in the text files)
- text file with text group ( labela) 
one can copy content of the text fil into svg file to have file with labels 

It is manual optimisation ( file size is smaller)
 
 

Draw a model of Mandelbrot set. 
Based on the paper M. Romera et al, Int. J. Bifurcation Chaos 13, 2279 (2003). 
Shrubs in the Mandelbrot Set Ordering 
www.tic.itefi.csic.es/gerardo/publica/Romera03.pdf

compare with : 
[[:File:Cactus_model_of_Mandelbrot_set.svg]]

angle is measured :
* in turns 
* counterclockwise 

https://en.wikipedia.org/wiki/Turn_(geometry)

angle can be :
* global ( 0 = 1 is on horizontal = x axis )
* local, where 0 is internal angle of parent 

CircleLevel :
0 = main circle (main cardioid) = period 1 cycle
1 = circles on the main circles
2 = 

stroke_width : A number without units is rendered the same as a percentage
https://css-tricks.com/almanac/properties/s/stroke-width/
Its default value is 1. 
If a <percentage> is used, the value represents a percentage of the current viewport. 
If a value of 0 is used the outline will never be drawn.

mutually and externally tangent circles
http://mathworld.wolfram.com/TangentCircles.html
Two circles are mutually and externally tangent if distance between their centers is equal to the sum of their radii

svg check
https://validator.w3.org/check

http://jsfiddle.net/j5ykdtgc/

cd existing_folder
git init
git remote add origin git@gitlab.com:adammajewski/mandelbrot-shrub-model.git
git add d.c
git commit -m "Initial commit"
git push -u origin master

*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h> //strncat
#include <math.h> // sin, M_PI

/* -  ------------------------- global variables ----------------------------- */

 // limits and main parameters of the program ~ detailes of the image
 // what is efficient relation between limits ?
 
 
 #define CircleLevelMax   10  // 
 const int dMax = 60;  // maximal denominator of internal angle 
 // svg 
 // viewBox
 double  Xmax = 1000.0;
 double  Ymax = 1000.0;
 double stroke_width; // = 0.5;
 double minimal_radius; /* radius of minimal circle to draw */
 double min_font_size;
 double min_length;
  
          
           
// RGB colors = 24 bit 
 // http://www.december.com/html/spec/color4.html  
char *black    = "#000000"; /* hexadecimal number as a string for svg color*/
char *white    = "#FFFFFF";
char *burgundy = "#9E0508";
char *SeaGreen = "#068481";
char *turquoise= "#40E0D0";
char *red      = "#FF0000";
char *navy     = "#000080";
char *blue     = "#0000FF";

	 	 
 // ----------- svg file  -----------------------------
 FILE * svg_fp;
 #define BUFSIZE 11 // 6+4+1 = name + .svg + eof
 char svg_name [BUFSIZE]; /* name of file */
 char *svg_file_name; //="shrub.svg";
 // text file
 FILE * txt_fp;
 char txt_name [BUFSIZE]; /* name of file */
 char *txt_file_name; //
 
 // ----------- svg parameters  -----------------
 char *comment = "<!-- sample comment in SVG file  \n can be multi-line -->";
 

 char *stroke_color ; 
 char *fill_color ;
 
 
 
 const double main_angle = 0.5; // angle in turns ; if 0.5 then main antenna is horizontal 
 int n1 = 1;
 int d1 = 1;
 
 int NoOfDrawnCircles = 0;
 int NoOfAllCircles = 0;
 int NoOfDrawnShrubs = 0;
 int NoOfAllShrubs = 0;
 
 
 // type for saving info about circle ( hyperbolic component of Mandelbrot set ) 
 
 typedef struct {
  int CircleLevel;
  
  // global (internal) angle in turns T = N/D
  double T;
  
  // local (internal) angle in turns t = n/d
  int n; // numerator
  int d; // denominator
  
  // ------- circle parameters -------
  // center of the circle c = cx + cy*I
  double cx;
  double cy;
  double radius;
  // ------  shrub parameters  -----------------
  // internal angle of the wake / limb / shrub  
  int N;
  int D;
  // apex = end point, it simulates M-Feigenbaum ponit, here shrub starts 
  // tangency point for n/d = 1/2 
  double ax;
  double ay;
  // 
  double length; // of the shrub
  
} DataType;

 
DataType Data = {0, 0.5,1,1 }; // declare and initialize the structure Data
DataType *DataP ;//= &Data = pointer to Data 
 
// --------------------------- functions -------------------------------===============

// GiveNewRadius(old,_p,_q):=old*abs(sin(%pi*_p/_q)/(_q*_q))

double GiveNewRadius(double OldRadius, int num, int den ){

  //double t = (double)num / den;
  double d2 = den*den;
  //printf( "t = %f ;  d2 = %f ",t, d2);
  //double temp =(num*num)/d2;
  //printf( "temp = %f ; ",temp);
  return 1.3*OldRadius/d2; //
}

//GiveGlobalAngle(old,new):=mod(old-(1/2-new),1);
double GiveNewAngle (double old, double n, double d) {
  double t =old-(main_angle-n/d);
  if (t > 1.0) t = t - 1.0;
  return t;
}

// pass the struct by reference. This means that your function can access the struct outside of the function and modify its values. You do this by passing a pointer to the structure to the function.
// https://stackoverflow.com/questions/10370047/passing-struct-to-function
void SetData(DataType *w, int CircleLevel, double T, int n, int d, double cx, double cy,  double radius, double ax, double ay,  double length ){

    w->CircleLevel = CircleLevel;
    w->T = T;
    w->n = n;
    w->d = d;
    w->cx = cx;
    w->cy = cy;
    w->radius = radius;
    w->N = 1;
    w->D = 1;
    
    w->ax = ax;
    w->ay = ay;
    w->length = length; }
       

     
  
void PrintData(DataType w){

  printf("CircleLevel = %d \n", w.CircleLevel);
  printf("global angle = T  = %f \n", w.T);
  printf("local angle  = n/d  = %d/%d \n", w.n, w.d);
  printf("center = (%f;%f) radius = %f \n", w.cx, w.cy, w.radius);
  printf("wake angle  = N/D  = %d/%d \n", w.N, w.D);
  
  printf("apex = (%f;%f) \n", w.ax,w.ay);
  printf("length = %f\n", w.length);
  
}

/*
w_new is changed and then it is used to draw = output
w_old is is an input 

it is the longest procedure !!!!
maybe some changes ?

*/
int UpdateData(DataType *w_new, DataType w_old, int n, int d){

 
 
  double distance ; // distance between centers of tangent circles
  double angle;         
  
  
  
  // global angle 
  w_new->T = GiveNewAngle(w_old.T, n, d);
  angle = 2.0*M_PI*w_new->T;
  w_new->radius = GiveNewRadius(w_old.radius, n, d );
  distance =  w_old.radius + w_new->radius; // distance between centers of tangent circles
  
  // new center
  w_new->cx = w_old.cx + distance* cos(angle);
  w_new->cy = w_old.cy + distance* sin(angle);
  
  
  if (w_old.CircleLevel < CircleLevelMax) 
    {w_new->CircleLevel = w_old.CircleLevel +1; // next circle 
     }
    
    else w_new->CircleLevel = w_old.CircleLevel; // go to the shrub not next circle
    
    
  // local angle in turns t = n/d
  w_new->n = n;
  w_new->d = d;
  
  
  
  
  
  
  
  
  NoOfAllCircles += 1;
  
  /* ------- shrub parameters -------------- */
  
  switch(w_old.CircleLevel) {
  
    case 0  : // w_new->CircleLevel = 1
      w_new->length = 2.0*w_new->radius;
      w_new->N = n;
      w_new->D = d;
      break; /*  */
	
  
  
   
   default : /* w_old.CircleLevel>0 && w_new->CircleLevel > 1 */
     // length
     if (w_new->n ==1 && w_new->d ==2)
        w_new->length = w_old.length + 2.0 * w_new->radius; // limb
        else w_new->length = 2.0 * w_new->radius; // sublimb
     // N/D  
     // if new angle = 1/2  
     if (n ==1 && d ==2) {
       // save old angle
       w_new->N = w_old.N;
       w_new->D = w_old.D;
       }
       else { // use new values 
         w_new->N = n;
         w_new->D = d;
       }
       
          
        
   } // switch

  
  // when shrub is drawn then no circle is drawn, so use previous circle parameters 
  // apex of the last circle aproximates MF point
  angle = 2.0*M_PI*w_old.T;
  w_new->ax = w_old.cx + w_old.radius * cos(angle);  //
  w_new->ay = w_old.cy + w_old.radius * sin(angle); //
          

 //if ( w_new->D > dMax) 
  // {printf("w_new->D = %d > dMax for CircleLevel = %d  d= %d \n",w_new->D,  w_new->CircleLevel, d);
  //  return 1; }
    
 return 0;   
 
  
}

double  min(double a, double b) {
  return a < b ? a : b;
}

void beginSVG(){ // 

 
  snprintf(svg_name, 6, "%02d_%02d", dMax, CircleLevelMax ); /*nn_nn + \0 = 5+1  */ 
  svg_file_name = strncat(svg_name,".svg",9 ); // 5+4 = length(name) + length(.ext)
  //
  
  

 svg_fp = fopen( svg_file_name,"w");
 if( svg_fp == NULL ) {printf (" svg file open error \n"); return ; }
 
 

 // SVG
 fprintf(svg_fp,
     "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n"
     "%s \n "
     "<svg width=\"40cm\" height=\"40cm\" viewBox=\"0 0 %.0f %.0f\"\n"
     " xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\">\n",
     comment,Xmax,Ymax);

  fprintf(svg_fp,"<g stroke=\"%s\" stroke-width=\"%f\" fill=\"%s\">\n",   stroke_color, stroke_width, fill_color); // group open

  printf(" beginSVG done \n");

}

int BeginTextFile(){

  snprintf(txt_name, 6, "%02d_%02d", dMax, CircleLevelMax ); /*nn_nn + \0 = 5+1  */
  txt_file_name = strncat(txt_name,".txt",9 ); // 5+4 = length(name) + length(.ext)
  //
  
  

 txt_fp = fopen( txt_file_name,"w");
 if( txt_fp == NULL ) {printf (" txt file open error \n"); return 1; } 
 fprintf(txt_fp,"<g text-anchor=\"middle\"  fill=\"black\" stroke =\"none\">\n"); // group open

 return 0;

};

void Begin(){

   // main circle = period 1 
  double  R1; // radius
 // center = x1 + y1*I;
 double X1;
 double Y1;

  // period 1 circle
  // center 
  X1 = Xmax - Xmax/3.5;
  Y1 = Ymax/2.0;
  R1 = min(X1, Y1)- min(X1, Y1)*0.5; // radius
  
  //
  DataP = &Data;
  SetData(DataP, 0, main_angle, n1, d1, X1, Y1, R1, X1-R1, Y1, 0.0 );
 
  stroke_color = black; // stroke
  fill_color = white; // fill
 
  // set limits : all depends on the stroke_width !!!!
  stroke_width = 0.5; //  a percentage of the current viewport
  minimal_radius = stroke_width/10.0; /* radius of minimal circle to draw */
  min_font_size = minimal_radius/2.0;
  min_length = minimal_radius;
 

  BeginTextFile();

  beginSVG();
}

void EndTxt(void){
 fprintf(txt_fp,"</g>\n"); // group close
 fclose(txt_fp); // file close 
 printf("EndTxt done,\n file %s saved \n",txt_file_name );
 
}

void EndSVG(void){
 fprintf(svg_fp,"</g>\n"); // group close
 fprintf(svg_fp,"</svg>\n"); // svg close
 fclose(svg_fp); // file close 
 printf("endSVG done\n file %s saved \n",svg_file_name );
 
}

/*

 draw text to the file in the svg format 

<text x="0" y="35" font-family="Verdana" font-size="35">
    Hello, out there
  </text>
  
  https://stackoverflow.com/questions/8865458/how-do-i-vertically-center-text-with-css?rq=1
  https://stackoverflow.com/questions/12250403/vertical-alignment-of-text-element-in-svg?rq=1
  https://www.w3.org/TR/SVG/text.html#TextElementDYAttribute

*/

void draw_txt(double x, double y,  double font_size, int n, int d  ){

       double dy; // vertical alligment inside a circle 
       
       if ( font_size > min_font_size){
         dy = font_size/3.0; // vertical alligment inside a circle 
         fprintf(txt_fp,"<text x=\"%.1f\" y=\"%.1f\" dy=\"%.2f\"  font-size=\"%.2f\">%d/%d</text>\n", x, y, dy, font_size, n,d);
        // font-family = \"serif\" font-weight=\"normal\" font-style=\"normal\"
        }
}

// draw circle to svg file
void draw_circle(double x, double y,  double radius ){

    if (radius > minimal_radius){
       fprintf(svg_fp,"<circle cx=\"%.1f\" cy=\"%.1f\" r=\"%.1f\" />\n", x, y, radius);
      // else printf ("radius < minimal_radius\n");
      NoOfDrawnCircles += 1; // count the number of drawn circles
      }
}

// draw circle to svg file
void draw_labeled_circle(DataType w){

    
      int n;
      double font_size;
      if (w.n==1 && w.d==1) 
               n=1 ;
               else n= w.d-w.n ; /* reverse y axis 
               because  in svg the initial coordinate system has the origin at the top/left 
               with the x-axis pointing to the right and the y-axis pointing down 
               https://www.w3.org/TR/SVGTiny12/coords.html#InitialViewport
               */
               
      draw_circle(w.cx, w.cy, w.radius);
      font_size = w.radius/2.0;
      if (font_size > min_font_size) 
        draw_txt   (w.cx, w.cy, font_size, n, w.d ); 
      
}

/*

<line x1="0" y1="0" x2="200" y2="200" style="stroke:rgb(255,0,0);stroke-width:2" />
*/

// draw line to svg file
void draw_line(double x1, double y1,   double T, double length ){

  double x2;
  double y2;
  double angle;
  
  angle  = 2.0*M_PI*T;
  
  x2 = x1 + length * cos(angle);
  y2 = y1 + length * sin(angle);
  fprintf(svg_fp,"<line x1=\"%.1f\" y1=\"%.1f\"  x2=\"%.1f\" y2=\"%.1f\"/>\n", x1, y1, x2, y2);
      
}

  
void draw_line_from_points(double x1, double y1,   double x2, double y2 ){

  fprintf(svg_fp,"<line x1=\"%.1f\" y1=\"%.1f\"  x2=\"%.1f\" y2=\"%.1f\"/>\n", x1, y1, x2, y2);
      
}

/* 

[[:File:N-th_arm_stars.svg]]
not working properly
---   bugs --------------
"Though there are some bugs (see lower side of 1/1->1/2, for
example)." Claude

*/
 
void DrawStar(DataType w){

 double cx,cy; // center of the star
 double length; // of the arms
 double angle; // between arms 
 double T; // global angle in turns 
 
 int n;
 
 length = w.length/3.0;
 T = 2.0*M_PI*w.T;
 cx = w.ax + length * cos(T);
 cy = w.ay + length * sin(T);
 angle = 1.0 / w.D; // angle in turns between w.D arms 
 
 // draw w.D arms = the star 
 // set up starting angle 
 T = w.T + 0.5;
 if (T > 1.0) T -= 1.0; // modulo 1 
  
 for(n=0; n<w.D; n++) {
    draw_line(cx, cy, T, length); 
   T = T + angle;
   if (T> 1.0) T -= 1.0; // modulo 1 
   }
}

/*

. The part of the Mandelbrot set connected to the main cardioid at this bifurcation point c_{p/q} is called the p/q-limb. 
*/

// shrub is a part of the limb after Myrberg-Feigenbaum point ( end of period doubling cascade )
void DrawShrub(DataType w){
// to do 

//PrintAddress(w);
NoOfAllShrubs += 1;
//(w.n == 1 && w.d == 2 &&
if( w.length > min_length) 
 {
   NoOfDrawnShrubs += 1;
   DrawStar(w);
  //else PrintData(w);
 }
  
}

/*
https://stackoverflow.com/questions/19738919/gcd-function-for-c
The GCD function uses Euclid's Algorithm. 
It computes A mod B, then swaps A and B with an XOR swap.
*/

int gcd(int a, int b)
{
    int temp;
    while (b != 0)
    {
        temp = a % b;

        a = b;
        b = temp;
    }
    return a;
}

/* 
  main drawing procedure 
  surprisingly short, but recurrent so maybe inefficient
  it differs from paper algorithm
  where first main cardioid is drawn then all n/d limbs ( families )
  here all 1/2-"limb"  with sublimbs ( approx)
  

*/

int drawSVG(DataType w)
{
 // internal angle = n/d 
 int n;  // numerator
 int d;  // denominator
 //
 DataType w_new; // 
 DataType *w_newP = &w_new; // to change w in a function (UpdateData) one must pass it by reference not value 
 
  
 if (w.CircleLevel > CircleLevelMax ||  w.radius < minimal_radius)  /* when to stop the recursion, 
 if  ||  w.radius < minimal_radius is removed then error */
   {  
       DrawShrub(w); // chaotic part of the set
       return 0;} // end of the recursion
 
 
 // periodic part of the set
 draw_labeled_circle(w); // draw one circle 
 // and smaller mutually and externally tangent circles 
 // n/d = local angle in turns
 for (d = 2; d <= dMax; ++d )
   for (n = 1; n < d; ++n )
     if (gcd(n,d)==1 )// irreducible fraction
       {
         UpdateData(w_newP, w, n, d);
         drawSVG(w_new); // recurence
       } //  if(gcd
     
  
  
  return 0;
}

void PrintInfo(DataType w){

 PrintData(w);  
 printf("limits \n");
 printf("Maximal Circle level = %d \n", CircleLevelMax);
 printf("Maximal denominator = %d \n", dMax);
 printf("width of the circle stroke = %f\n", stroke_width);
 printf("minimal radius of the circle = %f\n", minimal_radius);
 printf("min_font_size = %f\n", min_font_size);
 printf("min_length = %f\n", min_length);
  
 printf("\n");
 printf("there are %d circles drawn from all %d possible \n", NoOfDrawnCircles, NoOfAllCircles);
 printf("there are %d shrubs drawn from all %d possible \n", NoOfDrawnShrubs, NoOfAllShrubs);
}

void End(DataType w){

  EndTxt();
  EndSVG();
  PrintInfo(Data);

}

// ------------------------------------- main ---------------------------------------------------

int main(){
    

 
 
  
  Begin();
   
 
  drawSVG(Data);  
   
  End(Data);     
  
 
 
 return 0;
}

/*

gcc e.c -Wall -lm
./a.out

program creates 2 files 
- svg file ( without labels which are in the text files)
- text file with text group ( labela) 
one can copy content of the text fil into svg file to have file with labels 

It is manual optimisation ( file size is smaller)
 
 



Draw a model of Mandelbrot set. 
Based on the paper M. Romera et al, Int. J. Bifurcation Chaos 13, 2279 (2003). 
Shrubs in the Mandelbrot Set Ordering 
www.tic.itefi.csic.es/gerardo/publica/Romera03.pdf

image:
[[:File:Shrub_model_of_Mandelbrot_set_60_10_labelled.png]]


compare with : 
[[:File:Cactus_model_of_Mandelbrot_set.svg]]


angle is measured :
* in turns 
* counterclockwise 

https://en.wikipedia.org/wiki/Turn_(geometry)

angle can be :
* global ( 0 = 1 is on horizontal = x axis )
* local, where 0 is internal angle of parent 

CircleLevel :
0 = main circle (main cardioid) = period 1 cycle
1 = circles on the main circles
2 = 


stroke_width : A number without units is rendered the same as a percentage
https://css-tricks.com/almanac/properties/s/stroke-width/
Its default value is 1. 
If a <percentage> is used, the value represents a percentage of the current viewport. 
If a value of 0 is used the outline will never be drawn.




mutually and externally tangent circles
http://mathworld.wolfram.com/TangentCircles.html
Two circles are mutually and externally tangent if distance between their centers is equal to the sum of their radii

svg check
https://validator.w3.org/check

svg test
http://jsfiddle.net/j5ykdtgc/


cd existing_folder
git init
git remote add origin git@gitlab.com:adammajewski/mandelbrot-shrub-model.git
git add d.c
git commit -m "Initial commit"
git push -u origin master


----------------

file 06_02.txt saved (dMax_CircleLevelMax.txt)
file 06_02.svg saved (dMax_CircleLevelMax.svg)
w:
circle parameters : 
	CircleLevel = 0 
	global angle = T  = 0.500000 
	local angle  = n/d  = 1/1 
	center = (714.285714;500.000000) radius = 250.000000 
	wake angle  = N/D  = 1/1 
shrub/star parameters : 
	apex of the last circle = (464.285714;500.000000) = end point , it simulates M-Feigenbaum ponit, here shrub starts (base of the shrub/star))
	length of the shrub = 0.000000

limits 
Maximal Circle level = 2 
Maximal denominator = 6 
width of the circle stroke = 0.500000
minimal radius of the circle = 0.050000
min_font_size = 0.025000
min_length of the shrub= 0.050000

there are 2754 circles drawn from all 30294 possible 
there are 4235 shrubs drawn from all 27541 possible 

------------

limb = 1/2  
limb = 1/3  
limb = 2/3  
limb = 1/4  
limb = 3/4  
limb = 1/5  
limb = 2/5  
limb = 3/5  
limb = 4/5  
limb = 1/6  
limb = 5/6  
limb = 1/7  
limb = 2/7  
limb = 3/7  
limb = 4/7  
limb = 5/7  
limb = 6/7  
limb = 1/8  
limb = 3/8  
limb = 5/8  
limb = 7/8  
limb = 1/9  
limb = 2/9  
limb = 4/9  
limb = 5/9  
limb = 7/9  
limb = 8/9  
limb = 1/10  
limb = 3/10  
limb = 7/10  
limb = 9/10  
limb = 1/11  
limb = 2/11  
limb = 3/11  
limb = 4/11  
limb = 5/11  
limb = 6/11  
limb = 7/11  
limb = 8/11  
limb = 9/11  
limb = 10/11  
limb = 1/12  
limb = 5/12  
limb = 7/12  
limb = 11/12  
limb = 1/13  
limb = 2/13  
limb = 3/13  
limb = 4/13  
limb = 5/13  
limb = 6/13  
limb = 7/13  
limb = 8/13  
limb = 9/13  
limb = 10/13  
limb = 11/13  
limb = 12/13  
limb = 1/14  
limb = 3/14  
limb = 5/14  
limb = 9/14  
limb = 11/14  
limb = 13/14  
limb = 1/15  
limb = 2/15  
limb = 4/15  
limb = 7/15  
limb = 8/15  
limb = 11/15  
limb = 13/15  
limb = 14/15  
limb = 1/16  
limb = 3/16  
limb = 5/16  
limb = 7/16  
limb = 9/16  
limb = 11/16  
limb = 13/16  
limb = 15/16  
limb = 1/17  
limb = 2/17  
limb = 3/17  
limb = 4/17  
limb = 5/17  
limb = 6/17  
limb = 7/17  
limb = 8/17  
limb = 9/17  
limb = 10/17  
limb = 11/17  
limb = 12/17  
limb = 13/17  
limb = 14/17  
limb = 15/17  
limb = 16/17  
limb = 1/18  
limb = 5/18  
limb = 7/18  
limb = 11/18  
limb = 13/18  
limb = 17/18  
limb = 1/19  
limb = 2/19  
limb = 3/19  
limb = 4/19  
limb = 5/19  
limb = 6/19  
limb = 7/19  
limb = 8/19  
limb = 9/19  
limb = 10/19  
limb = 11/19  
limb = 12/19  
limb = 13/19  
limb = 14/19  
limb = 15/19  
limb = 16/19  
limb = 17/19  
limb = 18/19  
limb = 1/20  
limb = 3/20  
limb = 7/20  
limb = 9/20  
limb = 11/20  
limb = 13/20  
limb = 17/20  
limb = 19/20  
file 20_04.txt saved (dMax_CircleLevelMax.txt)
file 20_04.svg saved (dMax_CircleLevelMax.svg)

level 0 parameters from w data structure:

component/circle parameters : 
	CircleLevel = 0 
	global angle = T  = 0.500000 
	local angle  = n/d  = 1/1 
	center = (714.285714;500.000000) radius = 250.000000 
	wake angle  = N/D  = 1/1 

shrub/star parameters : 
	apex of the last circle = (464.285714;500.000000) = end point , it simulates M-Feigenbaum ponit, here shrub starts (base of the shrub/star))
	length of the shrub = 0.000000

limits 
	Maximal Circle level = 4 
	Maximal denominator = 20 
	width of the circle stroke = 0.500000
	minimal radius of the circle = 0.050000
	min_font_size = 0.025000
	min_length of the shrub= 0.050000
summary 
	there are 9768 circles drawn from all 262518652 possible 
	there are 23537 shrubs drawn from all 260451577 possible 

real	0m41,573s
user	0m41,561s
sys	0m0,009s




*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h> //strncat
#include <math.h> // sin, M_PI




/* -  ------------------------- global variables ----------------------------- */


 // limits and main parameters of the program ~ detailes of the image
 // what is efficient relation between limits ?
 
 const int dMax = 20;  //  60 maximal denominator of internal angle 
 #define CircleLevelMax  4// 10  // 
 
 // svg 
 // viewBox
 double  Xmax = 1000.0;
 double  Ymax = 1000.0;
 double stroke_width; // = 0.5;
 double minimal_radius; /* radius of minimal circle to draw */
 double min_font_size;
 double min_length;
  
          
           
// RGB colors = 24 bit 
 // http://www.december.com/html/spec/color4.html  
char *black    = "#000000"; /* hexadecimal number as a string for svg color*/
char *white    = "#FFFFFF";
char *burgundy = "#9E0508";
char *SeaGreen = "#068481";
char *turquoise= "#40E0D0";
char *red      = "#FF0000";
char *navy     = "#000080";
char *blue     = "#0000FF";

	 	 
 // ----------- svg file  -----------------------------
 FILE * svg_fp;
 #define BUFSIZE 11 // 6+4+1 = name + .svg + eof
 char svg_name [BUFSIZE]; /* name of file */
 char *svg_file_name; //="shrub.svg";
 // text file
 FILE * txt_fp;
 char txt_name [BUFSIZE]; /* name of file */
 char *txt_file_name; //
 
 // ----------- svg parameters  -----------------
 char *comment = "<!-- sample comment in SVG file  \n can be multi-line -->";
 

 char *stroke_color ; 
 char *fill_color ;
 
 
 
 const double main_angle = 0.5; // angle in turns ; if 0.5 then main antenna is horizontal 
 int n1 = 1;
 int d1 = 1;
 
 int NoOfDrawnCircles = 0;
 int NoOfAllCircles = 0;
 int NoOfDrawnShrubs = 0;
 int NoOfAllShrubs = 0;
 
 
 // type for saving info about circle ( hyperbolic component of Mandelbrot set ) 
 
 typedef struct {
  int CircleLevel;
  
  // global (internal) angle in turns T = N/D
  double T;
  
  // local (internal) angle in turns t = n/d
  int n; // numerator
  int d; // denominator
  
  // ------- circle parameters -------
  // center of the circle c = cx + cy*I
  double cx;
  double cy;
  double radius;
  // ------  shrub parameters  -----------------
  // internal angle of the wake / limb / shrub  
  int N;
  int D;
  // apex = end point, it simulates M-Feigenbaum ponit, here shrub starts 
  // tangency point for n/d = 1/2 
  double ax;
  double ay;
  // 
  double length; // of the shrub
  
} DataType;

 
DataType Data = {0, 0.5,1,1 }; // declare and initialize the structure Data
DataType *DataP ;//= &Data = pointer to Data 
 
// --------------------------- functions -------------------------------===============



// GiveNewRadius(old,_p,_q):=old*abs(sin(%pi*_p/_q)/(_q*_q))

double GiveNewRadius(double OldRadius, int num, int den ){

  //double t = (double)num / den;
  double d2 = den*den;
  //printf( "t = %f ;  d2 = %f ",t, d2);
  //double temp =(num*num)/d2;
  //printf( "temp = %f ; ",temp);
  return 1.3*OldRadius/d2; //
}


//GiveGlobalAngle(old,new):=mod(old-(1/2-new),1);
double GiveNewAngle (double old, double n, double d) {
  double t =old-(main_angle-n/d);
  if (t > 1.0) t = t - 1.0;
  return t;
}



// pass the struct by reference. This means that your function can access the struct outside of the function and modify its values. You do this by passing a pointer to the structure to the function.
// https://stackoverflow.com/questions/10370047/passing-struct-to-function
void SetData(DataType *w, int CircleLevel, double T, int n, int d, double cx, double cy,  double radius, double ax, double ay,  double length ){

    w->CircleLevel = CircleLevel;
    w->T = T;
    w->n = n;
    w->d = d;
    w->cx = cx;
    w->cy = cy;
    w->radius = radius;
    w->N = 1;
    w->D = 1;
    
    w->ax = ax;
    w->ay = ay;
    w->length = length; }
       



     
  
void PrintData(DataType w){

  printf("\nlevel 0 parameters from w data structure:\n\n");
  printf("component/circle parameters : \n");
  printf("\tCircleLevel = %d \n", w.CircleLevel);
  printf("\tglobal angle = T  = %f \n", w.T);
  printf("\tlocal angle  = n/d  = %d/%d \n", w.n, w.d);
  printf("\tcenter = (%f;%f) radius = %f \n", w.cx, w.cy, w.radius);
  printf("\twake angle  = N/D  = %d/%d \n", w.N, w.D);
  printf("\nshrub/star parameters : \n");
  printf("\tapex of the last circle = (%f;%f) = end point , it simulates M-Feigenbaum ponit, here shrub starts (base of the shrub/star))\n", w.ax,w.ay);
  printf("\tlength of the shrub = %f\n", w.length);
  printf("\n");
  
}



/*
w_new is changed and then it is used to draw = output
w_old is is an input 

it is the longest procedure !!!!
maybe some changes ?

*/
int UpdateData(DataType *w_new, DataType w_old, int n, int d){

 
 
  double distance ; // distance between centers of tangent circles
  double angle;         
  
  
  
  // global angle 
  w_new->T = GiveNewAngle(w_old.T, n, d);
  angle = 2.0*M_PI*w_new->T;
  w_new->radius = GiveNewRadius(w_old.radius, n, d );
  distance =  w_old.radius + w_new->radius; // distance between centers of tangent circles
  
  // new center
  w_new->cx = w_old.cx + distance* cos(angle);
  w_new->cy = w_old.cy + distance* sin(angle);
  
  
  if (w_old.CircleLevel < CircleLevelMax) 
    {w_new->CircleLevel = w_old.CircleLevel +1; // next circle 
     }
    
    else w_new->CircleLevel = w_old.CircleLevel; // go to the shrub not next circle
    
    
  // local angle in turns t = n/d
  w_new->n = n;
  w_new->d = d;
  
  
  
  
  
  
  
  
  NoOfAllCircles += 1;
  
  /* ------- shrub parameters -------------- */
  
  switch(w_old.CircleLevel) {
  
    case 0  : // w_new->CircleLevel = 1
      w_new->length = 2.0*w_new->radius;
      w_new->N = n;
      w_new->D = d;
      break; /*  */
	
  
  
   
   default : /* w_old.CircleLevel>0 && w_new->CircleLevel > 1 */
     // length
     if (w_new->n ==1 && w_new->d ==2)
        w_new->length = w_old.length + 2.0 * w_new->radius; // limb
        else w_new->length = 2.0 * w_new->radius; // sublimb
     // N/D  
     // if new angle = 1/2  
     if (n ==1 && d ==2) {
       // save old angle
       w_new->N = w_old.N;
       w_new->D = w_old.D;
       }
       else { // use new values 
         w_new->N = n;
         w_new->D = d;
       }
       
          
        
   } // switch

  
  // when shrub is drawn then no circle is drawn, so use previous circle parameters 
  // apex of the last circle aproximates MF point
  angle = 2.0*M_PI*w_old.T;
  w_new->ax = w_old.cx + w_old.radius * cos(angle);  //
  w_new->ay = w_old.cy + w_old.radius * sin(angle); //
          

 //if ( w_new->D > dMax) 
  // {printf("w_new->D = %d > dMax for CircleLevel = %d  d= %d \n",w_new->D,  w_new->CircleLevel, d);
  //  return 1; }
    
 return 0;   
 
  
}



double  min(double a, double b) {
  return a < b ? a : b;
}





void beginSVG(){ // 

 
  snprintf(svg_name, 6, "%02d_%02d", dMax, CircleLevelMax ); /*nn_nn + \0 = 5+1  */ 
  svg_file_name = strncat(svg_name,".svg",9 ); // 5+4 = length(name) + length(.ext)
  //
  
  


 svg_fp = fopen( svg_file_name,"w");
 if( svg_fp == NULL ) {printf (" svg file open error \n"); return ; }
 
 

 // SVG
 fprintf(svg_fp,
     "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n"
     "%s \n "
     "<svg width=\"40cm\" height=\"40cm\" viewBox=\"0 0 %.0f %.0f\"\n"
     " xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\">\n",
     comment,Xmax,Ymax);

  fprintf(svg_fp,"<g stroke=\"%s\" stroke-width=\"%f\" fill=\"%s\">\n",   stroke_color, stroke_width, fill_color); // group open

  //printf(" beginSVG done \n");

}


int BeginTextFile(){

  snprintf(txt_name, 6, "%02d_%02d", dMax, CircleLevelMax ); /*nn_nn + \0 = 5+1  */
  txt_file_name = strncat(txt_name,".txt",9 ); // 5+4 = length(name) + length(.ext)
  //
  
  


 txt_fp = fopen( txt_file_name,"w");
 if( txt_fp == NULL ) {printf (" txt file open error \n"); return 1; } 
 fprintf(txt_fp,"<g text-anchor=\"middle\"  fill=\"black\" stroke =\"none\">\n"); // group open



 return 0;

};

void Begin(){


   // main circle = period 1 
  double  R1; // radius
 // center = x1 + y1*I;
 double X1;
 double Y1;

  // period 1 circle
  // center 
  X1 = Xmax - Xmax/3.5;
  Y1 = Ymax/2.0;
  R1 = min(X1, Y1)- min(X1, Y1)*0.5; // radius
  
  //
  DataP = &Data;
  SetData(DataP, 0, main_angle, n1, d1, X1, Y1, R1, X1-R1, Y1, 0.0 );
 
  stroke_color = black; // stroke
  fill_color = white; // fill
 
  // set limits : all depends on the stroke_width !!!!
  stroke_width = 0.5; //  a percentage of the current viewport
  minimal_radius = stroke_width/10.0; /* radius of minimal circle to draw */
  min_font_size = minimal_radius/2.0;
  min_length = minimal_radius;
 

  BeginTextFile();

  beginSVG();
}



void EndTxt(void){
 fprintf(txt_fp,"</g>\n"); // group close
 fclose(txt_fp); // file close 
 printf("file %s saved (dMax_CircleLevelMax.txt)\n",txt_file_name );
 
}


void EndSVG(void){
 fprintf(svg_fp,"</g>\n"); // group close
 fprintf(svg_fp,"</svg>\n"); // svg close
 fclose(svg_fp); // file close 
 printf("file %s saved (dMax_CircleLevelMax.svg)\n",svg_file_name );
 
}



/*

 draw text to the file in the svg format 

<text x="0" y="35" font-family="Verdana" font-size="35">
    Hello, out there
  </text>
  
  https://stackoverflow.com/questions/8865458/how-do-i-vertically-center-text-with-css?rq=1
  https://stackoverflow.com/questions/12250403/vertical-alignment-of-text-element-in-svg?rq=1
  https://www.w3.org/TR/SVG/text.html#TextElementDYAttribute

*/

void draw_txt(DataType w ){





       
      double x = w.cx;
      double y = w.cy;
      double font_size = w.radius/2.0;
      double dy; // vertical alligment inside a circle
      int n ;
      int d = w.d;
      
            
      if (w.n==1 && w.d==1) 
               n=1 ;
               else n= w.d-w.n ; /* reverse y axis 
               because  in svg the initial coordinate system has the origin at the top/left 
               with the x-axis pointing to the right and the y-axis pointing down 
               https://www.w3.org/TR/SVGTiny12/coords.html#InitialViewport
               */
      
           
       if ( font_size > min_font_size){
         dy = font_size/3.0; // vertical alligment inside a circle 
         fprintf(txt_fp,"<text x=\"%.1f\" y=\"%.1f\" dy=\"%.2f\"  font-size=\"%.2f\">%d/%d</text>\n", x, y, dy, font_size, n,d);
        // font-family = \"serif\" font-weight=\"normal\" font-style=\"normal\"
        }
}



// draw circle to svg file
void draw_circle(DataType w ){


     double x = w.cx;
     double y = w.cy;
     double radius = w.radius;
    

    if (radius > minimal_radius){
       // inf about progres 
       if (w.CircleLevel == 1) 
       		printf("limb = %d/%d  \n", w.n, w.d);		    
       		
       	
       fprintf(svg_fp,"<circle cx=\"%.1f\" cy=\"%.1f\" r=\"%.1f\" />\n", x, y, radius);
      // else printf ("radius < minimal_radius\n");
      NoOfDrawnCircles += 1; // count the number of drawn circles
      }
}


// https://www.w3.org/TR/SVG/paths.html#PathElement

/*
smooth curves in svg:
* Bezier curves, cubic ( C)  and quadratic (Q)
* arc = sections of circles or ellipses ( A)
* path 

*/



/* 
	https://en.wikipedia.org/wiki/Cardioid
	cusp at 0.0, 
	x(t) = 2a (1 - cos(t))* cos(t)
	y(t) = 2a (1 - cos(t))* sin(t)
*/
// compute points of the cardioid

// save points as a svg path : 

// <path id="cardioid" d="M 100 100 L 300 100 L 200 300 z" fill="orange" stroke="black" stroke-width="3" />
 // fprintf(svg_fp,"<circle cx=\"%.1f\" cy=\"%.1f\" r=\"%.1f\" />\n", x, y, radius);

/*
void draw_cardioid(DataType w){

 draw_circle(w); // temporary solition

}

*/


// draw circle to svg file
void draw_labeled_component(DataType w){

    
      
      
      
       draw_circle(w);// component
       draw_txt   ( w ); // label
      
}






/*

<line x1="0" y1="0" x2="200" y2="200" style="stroke:rgb(255,0,0);stroke-width:2" />
*/

// draw line to svg file
void draw_line(double x1, double y1,   double T, double length ){

  double x2;
  double y2;
  double angle;
  
  angle  = 2.0*M_PI*T;
  
  x2 = x1 + length * cos(angle);
  y2 = y1 + length * sin(angle);
  fprintf(svg_fp,"<line x1=\"%.1f\" y1=\"%.1f\"  x2=\"%.1f\" y2=\"%.1f\"/>\n", x1, y1, x2, y2);
      
}




  
void draw_line_from_points(double x1, double y1,   double x2, double y2 ){

  fprintf(svg_fp,"<line x1=\"%.1f\" y1=\"%.1f\"  x2=\"%.1f\" y2=\"%.1f\"/>\n", x1, y1, x2, y2);
      
}




/* 

[[:File:N-th_arm_stars.svg]]
not working properly
---   bugs --------------
"Though there are some bugs (see lower side of 1/1->1/2, for
example)." Claude


*/
 
void DrawStar(DataType w){

 double cx,cy; // center of the star
 double length; // of the arms
 double angle; // between arms 
 double T; // global angle in turns 
 
 int n;
 
 length = w.length/3.0;
 T = 2.0*M_PI*w.T;
 cx = w.ax + length * cos(T);
 cy = w.ay + length * sin(T);
 angle = 1.0 / w.D; // angle in turns between w.D arms 
 
 // draw w.D arms = the star 
 // set up starting angle 
 T = w.T + 0.5;
 if (T > 1.0) T -= 1.0; // modulo 1 
  
 for(n=0; n<w.D; n++) {
    draw_line(cx, cy, T, length); 
   T = T + angle;
   if (T> 1.0) T -= 1.0; // modulo 1 
   }
}










/*

. The part of the Mandelbrot set connected to the main cardioid at this bifurcation point c_{p/q} is called the p/q-limb. 
*/

// shrub is a part of the limb after Myrberg-Feigenbaum point ( end of period doubling cascade )
void DrawShrub(DataType w){
// to do 

//PrintAddress(w);
NoOfAllShrubs += 1;
//(w.n == 1 && w.d == 2 &&
if( w.length > min_length) 
 {
   NoOfDrawnShrubs += 1;
   DrawStar(w);
  //else PrintData(w);
 }
  
}




/*
https://stackoverflow.com/questions/19738919/gcd-function-for-c
The GCD function uses Euclid's Algorithm. 
It computes A mod B, then swaps A and B with an XOR swap.
*/

int gcd(int a, int b)
{
    int temp;
    while (b != 0)
    {
        temp = a % b;

        a = b;
        b = temp;
    }
    return a;
}



/* 
  main drawing procedure 
  surprisingly short, but recurrent so maybe inefficient
  it differs from paper algorithm
  where first main cardioid is drawn then all n/d limbs ( families )
  here all 1/2-"limb"  with sublimbs ( approx)
  


*/

int drawSVG(DataType w)
{
 // internal angle = n/d 
 int n;  // numerator
 int d;  // denominator
 //
 DataType w_new; // 
 DataType *w_newP = &w_new; // to change w in a function (UpdateData) one must pass it by reference not value 
 
  
     
  if (w.CircleLevel < CircleLevelMax ||  w.radius > minimal_radius)  /* when to stop the recursion,  if  it is removed then error */
   {
 
 	// periodic part of the set
 	draw_labeled_component(w); // draw one base component
 	// and smaller mutually and externally tangent circles 
 	// n/d = local angle in turns
 	for (d = 2; d <= dMax; ++d )
   		for (n = 1; n < d; ++n )
     			if (gcd(n,d)==1 )// irreducible fraction
       			{
         			UpdateData(w_newP, w, n, d);
         			drawSVG(w_new); // recurence
       			} //  if(gcd
       
   } 
   // stop    
   else {  
   	DrawShrub(w); // chaotic part of the set
   	return 0;} // end of the recursion    
 
 
  
  return 0;
}


void PrintInfo(DataType w){

 PrintData(w);  
 printf("limits \n");
 printf("\tMaximal Circle level = %d \n", CircleLevelMax);
 printf("\tMaximal denominator = %d \n", dMax);
 printf("\twidth of the circle stroke = %f\n", stroke_width);
 printf("\tminimal radius of the circle = %f\n", minimal_radius);
 printf("\tmin_font_size = %f\n", min_font_size);
 printf("\tmin_length of the shrub= %f\n", min_length);
  
 printf("summary \n");
 printf("\tthere are %d circles drawn from all %d possible \n", NoOfDrawnCircles, NoOfAllCircles);
 printf("\tthere are %d shrubs drawn from all %d possible \n", NoOfDrawnShrubs, NoOfAllShrubs);
}


void End(DataType w){

  EndTxt();
  EndSVG();
  PrintInfo(Data);

}



// ------------------------------------- main ---------------------------------------------------


int main(){
    

 
 
  
  Begin();
   
 
  drawSVG(Data);  // 
   
  End(Data);     
  
 
 
 return 0;
}