/*

gcc d.c -Wall -lm


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 

compare with : 
[[:File:Circle 7 4.svg|Circle 7 4.svg]]


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

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

*/

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




/* -  ------ global variables ---------- */
 
          
           
// RGB colors
 // http://www.december.com/html/spec/color4.html  
char *black    = "#000"; /* hexadecimal number as a string for svg color*/
char *white    = "#FFF";
char *burgundy = "#9E0508";
char *SeaGreen = "#068481";
char *turquoise= "#40E0D0";
char *red      = "#F00";
char *navy     = "#000080";
char *blue     = "#00F";

	 	 
 // svg file 
 FILE * fp;
 char *filename; //="shrub.svg";
 // svg parameters 
 char *comment = "";
 double stroke_width = 0.5;
 double minimal_radius; /* radius of minimal circle to draw */

 char *stroke_color ; 
 char *fill_color ;
 
 double  Xmax = 1000.0;
 double  Ymax = 1000.0;
 
 
 
 
 // main circle = period 1 
  double  r1; // radius
 // center = x1 + y1*I;
 double X1;
 double Y1;
 
 
 const int iLevelMax = 4; // 
 const int dMax = 14;  // maximal denominator of internal angle 
 const double main_angle = 0.5; // in turns
 
// ----------------- functions -------------------------------

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



void beginSVG(void){


  char name [30]; /* name of file */
  snprintf(name,20,"%d_%d",dMax, iLevelMax ); /*  */
  filename =strncat(name,".svg", 24);



 fp = fopen( filename,"w");
 if( fp == NULL ) {printf (" file open error \n"); return ; }
 stroke_color = black; // stroke
 fill_color = white; // fill
 minimal_radius = stroke_width; /* radius of minimal circle to draw */
 
 // period 1 circle
 // center 
 X1 = Xmax - Xmax/3.0;
 Y1 = Ymax/2.0;
 r1 = min(X1, Y1)- 2.0*min(X1, Y1)/3.0; // radius
 


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

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

printf(" beginSVG done \n");

}




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



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

    
       fprintf(fp,"<circle cx=\"%.1f\" cy=\"%.1f\" r=\"%.1f\" />\n", x, y, radius);
      // else printf ("radius < minimal_radius\n");
}





  
/*
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;
}



// 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-=1.0;
  return t;
}





/*
    test_v =Test(x1,y1, r1, x,y, Radius);
    printf(" n/d = %d / %d \n", num, den );
    printf( "main circle center = (%f ; %f) radius = %f \n", x1,y1, r1);
    printf( "new  circle center = (%f ; %f) radius = %f \n", x,y, Radius);
    printf( "test = %f \n", test_v);
    */
long double Test(double x1, double y1, double r1, double x2, double y2, double r2){

return (sqrt((x1-x2)*(x1-x2)+ (y1-y2)*(y1-y2)) - (r1+r2));

}




// shrub is a part of the limb after Myrberg-Feigenbaum point ( end of period doubling cascade )
void DrawShrub(double n, double d){

}

/*

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




// t = global angle in turns
void drawSVG(double t,  int iLevel, double x, double y, double r)
{
 // internal angle = n/d 
 int n;  // numerator
 int d;  // denominator
 
 
 double r_new; // Radius;
 double x_new,y_new; // center
 double distance; // new_radius + old_radius
 double t_local;
 int iLevel_new;

    
 
 
 
  draw_circle(x,y,r); 
 
 
 
   for (d = 2; d <= dMax; ++ d ){
     
     for (n = 1; n < d; ++ n ){
       if (gcd(n,d)==1 )// irreducible fraction
            {
            
           if ( iLevel > iLevelMax) { 
               DrawShrub(n,d);
               break;}
           iLevel_new = iLevel +1;
           
           // num/den * 1/2  where 1/2 is a local angle 
           // period doubling cascade 
           // new radius
           r_new = GiveNewRadius(r, n, d );
           if (r_new > minimal_radius) {
           //
           distance =  r + r_new  ;
           //printf("n = %d; d = %d iRadius = %d \n", num, den, iRadius);
           // new center
           t_local = GiveNewAngle(t, n, d);
           x_new = x + distance* cos(2.0*M_PI*t_local);
           y_new = y + distance* sin(2.0*M_PI*t_local);
           //
           drawSVG(t_local, iLevel_new, x_new, y_new, r_new);  
           }
           
           
           
           } //  if(gcd
     } // for (n
   } // for(d 
   
   
  
  
  
}


// ------------------------------------- main ---------------------------------------------------
int main(){
    
 

 
 
 
 
   
 beginSVG();  
 
 drawSVG(main_angle, 0, X1, Y1, r1);  
 
       
 EndSVG();
 
 
 return 0;
}