/*


It's a fun map. Here's one for the Fbm point defined by going to period 3 mini, then it's relative period 3 mini, etc.
 gerrit


https://fractalforums.org/fractal-image-gallery/18/feigenbaum-stretch/2112
https://fractalforums.org/programming/11/m-streching-feigenbaum-from-mandelbrot-graphics-library-by-claude/3077



m-streching-feigenbaum.c

Transform is c ←e^c+c0 
with c0=−1.401155189092051.





gcc s3.c -lm -Wall -fopenmp -lpari -lmpc -lmpfr -lgmp -lm $(pkg-config --cflags --libs cairo)




=======================
centers of satellites cardioids ( on the main antenna )
c = -1.754877666246693  period = 3
c = -1.785865646410674  period = 9
c = -1.786429858055761  period = 27
c = -1.786440067359078  period = 81
c = -1.786440252157044  period = 243
c = -1.7864402555041126 period = 729
c = -1.7864402555625307 period = 2187
c = -1.786440255563619017 period = 6561
c = -1.7864402555636388747
c = -1.7864402555;
-------------------------------------------------------------

--------------------------------------------------
on main antenna ( real slice of mandelbrot set ):

period 3 (3/7,4/7)  cusp c= -1.75 = -7/4 center c=-1.754877666246693
period 9 .(001001010) .(001010001)"



3
[3/7,4/7]
9
[228/511,283/511]
27
[59914011/134217727,74303716/134217727]
81
[1079314880428363211458788/2417851639229258349412351,       1338536758800895137953563/2417851639229258349412351]
243
[6309681859766642153488012706542361061031341827515755455785608209904056091/14134776518227074636666380005943348126619871175004951664972849610340958207, 7825094658460432483178367299400987065588529347489196209187241400436902116 /14134776518227074636666380005943348126619871175004951664972849610340958207]

(%o22) [1260623372715624251596526631747370048404057087743061433025671746894079267500591465669860569986636462919899527462685012516848740080247243839800204694110235562302586056652374438695747306836417796879240393278534689460439268
         /2824013958708217496949108842204627863351353911851577524683401930862693830361198499905873920995229996970897865498283996578123296865878390947626553088486946106430796091482716120572632072492703527723757359478834530365734911,
        1563390585992593245352582210457257814947296824108516091657730183968614562860607034236013351008593534050998338035598984061274556785631147107826348394376710544128210034830341681876884765656285730844516966200299840905295643
         /2824013958708217496949108842204627863351353911851577524683401930862693830361198499905873920995229996970897865498283996578123296865878390947626553088486946106430796091482716120572632072492703527723757359478834530365734911]
(%i23) t:h2(t)
(%o23) [10053540528705587606144343964571395857305739074823089969864347614424856359955693316255490127296397738095346834867196211074086432187862300532851025254229333089534629081303899791323025713592226561693490194031326260764758759553574780477087039336817214120161773323134578467312898656865849059615104468133377000807189815685871691674981838136633883396727308512724284079919417313824331180651588488371172994635146902928538559653853188574520399426562385056331207135050181559504990454855980461686322800390811892846857617711714719180724079051842494381856293521988999830348460009326896820106553672108357391760790443678125366138873892416952301039236023546099897030145095451
         /22521666186739810716017129863546474247925813985186306383480395759255528226899398031271474259176472311138130302909797370563984143657840869015723623128208458564711932390036470899999940399003451938771468124775412134017692556582989981820657376622292117781429636251681663364133246565240603802421309871873386179061056683850002592550658912005406539073633368208169912295100765107813056035892005304712416817301893705332191267615274619407764899319819287315641053762034009508858938727147166108668063720886082155921552955435591503022900749768195878181491707144967461286550140023171580051111358071521188926359524314904062599527129957121490685127011258936786658012354838527,
        12468125658034223109872785898975078390620074910363216413616048144830671866943704715015984131880074573042783468042601159489897711469978568482872597873979125475177303308732571108676914685411225377077977930744085873252933797029415201343570337285474903661267862928547084896820347908374754742806205403740009178253866868164130900875677073868772655676906059695445628215181347793988724855240416816341243822666746802403652707961421430833244499893256902259309846626983827949353948272291185646981740920495270263074695337723876783842176670716353383799635413622978461456201680013844683231004804399412831534598733871225937233388256064704538384087775235390686760982209743076
         /22521666186739810716017129863546474247925813985186306383480395759255528226899398031271474259176472311138130302909797370563984143657840869015723623128208458564711932390036470899999940399003451938771468124775412134017692556582989981820657376622292117781429636251681663364133246565240603802421309871873386179061056683850002592550658912005406539073633368208169912295100765107813056035892005304712416817301893705332191267615274619407764899319819287315641053762034009508858938727147166108668063720886082155921552955435591503022900749768195878181491707144967461286550140023171580051111358071521188926359524314904062599527129957121490685127011258936786658012354838527]

*/


#include "mandelbrot-numerics.h"
#include "mandelbrot-graphics.h"
//#include "mandelbrot-symbolics.h"
#include "m_d_transform.c"
#include "m_d_mat2.c"
#include "m_d_util.h"
#include "m_pixel.c"
#include "m_d_colour.c"
#include "m_d_render_scanline.c"
//#include "m_binangle.c"
#include "m_d_exray_in.c"
#include <stdio.h>
#include <string.h>		// strcat

// #include <pari/pari.h>



// colors
  m_pixel_t black; // = m_pixel_rgba(0, 0, 0, 1);
  m_pixel_t white; // = m_pixel_rgba(1, 1, 1, 1);
  m_pixel_t red  ; // = m_pixel_rgba(1, 0, 0, 1);
  
// m_d_transform.c
  m_d_transform *rect ; //= m_d_transform_rectangular(w, h, - m * r, r);
  m_d_transform *exponential; // = m_d_transform_exponential(c); // Here's one for the Fbm point defined by going to period 3 mini, then it's relative period 3 mini, etc.
  m_d_transform *transform ; //= m_d_transform_compose(rect, exponential);
  m_d_colour_t *colour ; //= m_d_colour_minimal(white, black, white);
  m_image *image ; //= m_image_new(w, h);  
  



/*

input
* s_b_angle = binary angle as a string ( for cairo_text_path)
* s_q_angle = the same angle as a decimal rational number ( for m_d_exray_in_new)
*/
void draw_external_ray(m_image *image, m_d_transform *transform, const char *s_b_angle, const char *s_q_angle, m_pixel_t colour, double dx, double dy, int iMax) {
  //int maxiters = 5024;
  double r = sqrt(2);

  // rational angle 
  mpq_t q_angle;
  mpq_init(q_angle);
  mpq_set_str(q_angle, s_q_angle, 10);
  m_d_exray_in *ray = m_d_exray_in_new(q_angle, 8);
  mpq_clear(q_angle);

  cairo_surface_t *surface = m_image_surface(image);
  cairo_t *cr = cairo_create(surface);
  cairo_set_source_rgba(cr, m_pixel_red(colour), m_pixel_green(colour), m_pixel_blue(colour), m_pixel_alpha(colour));
  bool first = true;
  
  
  for (int i = 0; i < iMax; ++i) {
    if (m_failed == m_d_exray_in_step(ray, 64)) {
      break;
    }
    double _Complex c = m_d_exray_in_get(ray);
    
    if (cabs(c + 0.75) > r) {
      continue;
    }
  
  
    //double t = carg(c +3.14 ); // +0.75
    double t = cimag(c) > 0 ? twopi / 4 : -twopi / 4;
 
    double _Complex dc = 1;
    m_d_transform_reverse(transform, &c, &dc);
    if (first) {
      cairo_save(cr);
      cairo_translate(cr, creal(c) + dx, cimag(c) + dy);
      cairo_rotate(cr, -t);
      cairo_select_font_face(cr, "LMMono10", CAIRO_FONT_SLANT_NORMAL, CAIRO_FONT_WEIGHT_NORMAL);
      cairo_set_font_size(cr, 15);
      //cairo_text_path(cr, s_b_angle); 
      cairo_fill(cr);
      cairo_restore(cr);
      cairo_move_to(cr, creal(c) + dx, cimag(c) + dy);
      first = false;
    } else {
      cairo_line_to(cr, creal(c), cimag(c));
    }
  }
  cairo_stroke(cr);
  cairo_destroy(cr);
}



void ray(const char *s_q_angle,int iMax){

	draw_external_ray(image, transform, "",s_q_angle, red, 0, 0, iMax);
}



extern int main(int argc, char **argv) {
  (void) argc; // https://stackoverflow.com/questions/21045615/what-does-voidvar-actually-do
  (void) argv;
  
  // image size in pixels
  int h = 1500;
  int w = 3*h;
  //
  double r = 4.0;
  double m = 2.5;
  
  // 
  double er = 100;
  int maxiters = 100100;
  double c =  -1.7864402555636388747; //-1.7864402555625307; //-1.7864402555041126;//-1.7864402556782493; // -1.78644029;//-1.786440254937780;//-1.786440254938960 ;//-1.786440254886202;//-1.786440691627496; // -1.786440690499556; //-1.786437960060477; // -1.786723048596194; 
  // file name = 3_c_r_m.png  
  char name[50];
  snprintf(name, 50, "s3_%d_%.16f_%.1f_%.1f", h, fabs(c), r, m );
  const char *filename = strncat(name,".png", 4); //"1781.png";
 
 
  // colors
  black = m_pixel_rgba(0, 0, 0, 1);
  white = m_pixel_rgba(1, 1, 1, 1);
  red   = m_pixel_rgba(1, 0, 0, 1);
  
  // m_d_transform.c
  rect = m_d_transform_rectangular(w, h, - m * r, r);
  exponential = m_d_transform_exponential(c); // Here's one for the Fbm point defined by going to period 3 mini, then it's relative period 3 mini, etc.
  transform = m_d_transform_compose(rect, exponential);
  colour = m_d_colour_minimal(white, black, white);
  image = m_image_new(w, h);
    

  m_d_render_scanline(image, transform, er, maxiters, colour);
  
   // period 1
  ray("0/1",  400);
  ray("1/1",  400);
  // period 3
  ray("3/7",  400);
  ray("4/7",  400);
  // period 9
  ray("228/511", 800);
  ray("283/511", 800);
  // period 27 
  ray("59914011/134217727", 2200);
  ray("74303716/134217727", 2200);
  // period 81 
  ray("1079314880428363211458788/2417851639229258349412351", 6000);
  ray("1338536758800895137953563/2417851639229258349412351", 6000);
  // period 243  
  ray("6309681859766642153488012706542361061031341827515755455785608209904056091/14134776518227074636666380005943348126619871175004951664972849610340958207", 18000);
  ray("7825094658460432483178367299400987065588529347489196209187241400436902116 /14134776518227074636666380005943348126619871175004951664972849610340958207", 18000);
  
  m_image_save_png(image, filename);
  //
  return 0;
}