ultrafractal - translationg coloring algorithm ( ucl) to c

  • 4 Replies
  • 287 Views

0 Members and 1 Guest are viewing this topic.

Offline Adam Majewski

  • *
  • Strange Attractor
  • ******
  • Posts: 93
« on: January 04, 2018, 04:22:43 PM »
Hi,

I have mmf3.ucl http://www.fractalforums.com/programming/smooth-external-angle-of-mandelbrot-set/   for Ultra Fractal by David Makin, It works .


Code: [Select]
MMF3-FieldLines {
;
; This version of calculating field lines will work reasonably well with both
; Mandelbrot and Julia Sets for divergent fractals with a divergence that is
; close to being a positive integer >=2 but can also produce quite interesting
; results with fractals that do not fit that criteria - just not rendering the
; field lines correctly in other cases :-)


global:
  float twopi = 2.0*#pi
  float dtwopi = 0.5/#pi
  float dp = 1.0/@power
  float m = @power*dtwopi
  float dp2 = 1.0/(2.0+@power) ; Changed from originally trying 1/power
 
 
 
init:
  complex zv[#maxiter]
  int i = 0
  float a = 0.0
  float b = 0.0
  float c = 0.0
  float s = 0.0
  if @julia
    zv[i] = #z
    i = i + 1
  endif
 
 
 
loop:
  zv[i] = #z
  i = i + 1
 
 
 
final:

  if (s = atan2(#z))<0
    s = s + twopi
  endif
 
  while i > @skipiter
    i = i - 1
   
    if @version==0
      b = (@power*atan2(zv[i]))%twopi
    else
      b = atan2(zv[i]^@power)
    endif
   
   
    if b < 0.0
      b = b + twopi
    endif
   
    c = c + b - s
   
    if c >= #pi
      s = s + twopi
      if @fixit || @accident
        c = c - twopi
      endif
     
    elseif c < -#pi
      s = s - twopi
      if @fixit || @accident
        c = c + twopi
      endif
    endif
   
   
    if @fixit
      c = dp2*c
    elseif !@accident
      c = 0.0
    endif
   
    if (a = atan2(zv[i])) < 0.0
      a = a + twopi
    endif
   
    s = dp*(s + twopi*floor(a*m))
  endwhile
 
 
 
  s = s*dtwopi
 
  if @fixskip
    while s<0.0
      s = s + 1.0
    endwhile
   
    while s>=1.0
      s = s - 1.0
    endwhile
   
  endif
 
  #index = s
 
 
 
 
default:
  title = "MMF3 Field Lines (use high bailouts)"
  heading
    caption = "Information"
;    text = "This colouring is a little touchy around certain values, if you \
;            get obvious lines or dots that shouldn't be there but the 'Degree \
;            of Divergence' is set properly for the main formula and you are \
;            using a high bailout of say 1e20 or higher then try \
;            adjusting the location very slightly. For example if vertical \
;            lines appear on an unrotated fractal then try changing the x \
;            (real) location very slightly, if horizontal lines appear then try \
;            changing the y (imag) location very slightly."
  endheading
  heading
  endheading
  param version
    caption = "Version"
    default = true
    hint = ""
  endparam
  param power
    caption = "Degree of divergence"
    default = 2.0
    hint = "This is an estimate of the divergence of the main formula, e.g. 2 for z^2+c, 3 for z^3+c etc."
  endparam
  param julia
    caption = "Julia ?"
    default = false
    hint = "Enable for correct rendering of Julia Sets."
  endparam
  param skipiter
    caption = "'Start' iteration"
    default = 0
    min = 0
    hint = ""
  endparam
  param fixskip
    caption = "Fix 'Start' iteration"
    default = false
    hint = ""
  endparam
  param fixit
    caption = "Fix Field Lines"
    default = false
    hint = ""
  endparam
  param accident
    caption = "Happy Accident"
    default = false
    hint = ""
    visible = !@fixit
  endparam
}
I want to translate it to stand alone c program. Here is my version :
Code: [Select]

/*


   http://linas.org/art-gallery/escape/phase/phase.html
   image
   http://linas.org/art-gallery/escape/phase/phase.gif
 
   c console program:
   
   
   --------------------------------
   1. draws Mandelbrot set for Fc(z)=z*z +c
   using Mandelbrot algorithm ( boolean escape time )
   -------------------------------         
   2. technique of creating ppm file is  based on the code of Claudio Rocchini
   http://en.wikipedia.org/wiki/Image:Color_complex_plot.jpg
   create 24 bit color graphic file ,  portable pixmap file = PPM
   see http://en.wikipedia.org/wiki/Portable_pixmap
   to see the file use external application ( graphic viewer)
   
   
   complex point c -> virtual 2D array -> memory 1D array -> ppm file on the disc -> png file
   
   C -> pixel (iX,iY)  -> index k  -> 24bit color
   
   -----
   https://stackoverflow.com/questions/6418807/how-to-work-with-complex-numbers-in-c
   complex numbers are built in type
 
 --------------
 formated with emacs
   -------------
   to compile :

 
 
   gcc f.c -lm -Wall
 
 
   ./a.out
   
   
   to convert to png using ImageMagic

   convert phase.ppm phase.png 



   ----------------------
   git add  phase.png phase.c README.md
   git commit -m "phase"
   git push -u origin master



 
 
*/
#include <stdio.h>
#include <stdlib.h> // malloc
#include <math.h>
#include <complex.h> // https://stackoverflow.com/questions/6418807/how-to-work-with-complex-numbers-in-c
 


 
/* screen ( integer) coordinate */

const int iWidth  = 1000;
const int iHeight = 1000;


/* world ( double) coordinate = parameter plane*/
// double complex C =  Cx + Cy*I ;
const double CxMin=-2.4;
const double CxMax= 1.5;
const double CyMin=-1.95;
const double CyMax= 1.95;

/* */
double PixelWidth; //=(CxMax-CxMin)/iWidth;
double PixelHeight; // =(CyMax-CyMin)/iHeight;


/* color component ( R or G or B) is coded from 0 to 255 */
/* it is 24 bit color RGB file */
int ColorBytes = 3; // 3*8 = 24 bit color       

/* iterations  */
const int IterationMax=400;

/* bail-out value , radius of circle ;  */
const double EscapeRadius=3.0;

/*

param version     caption = "Version"     default = true     hint = ""   endparam

param power     caption = "Degree of divergence"     default = 2.0     hint = "This is an estimate of the divergence of the main formula, e.g. 2 for z^2+c, 3 for z^3+c etc."   endparam

param julia     caption = "Julia ?"     default = false     hint = "Enable for correct rendering of Julia Sets."   endparam

param skipiter     caption = "'Start' iteration"     default = 0     min = 0     hint = ""   endparam

param fixskip     caption = "Fix 'Start' iteration"      default = false    hint = ""   endparam


param fixit     caption = "Fix Field Lines"     default = false     hint = ""   endparam

param accident     caption = "Happy Accident"     default = false     hint = ""     visible = !@fixit   endparam

*/

double power = 2.0;
int version = 1;
int julia = 1;
int skipiter = 0;
int fixskip = 0;
int fixit = 0;
int accident = 0;

/*

 
global:
  float TwoPi = 2.0*#pi
  float dTwoPi = 0.5/#pi
  float dp = 1.0/@power
  float m = @power*dTwoPi
  float dp2 = 1.0/(2.0+@power) ; Changed from originally trying 1/power


 

*/

double TwoPi=2.0*M_PI; // float TwoPi = 2.0*#pi
double dTwoPi = 0.5/M_PI;  // float dTwoPi = 0.5/#pi
double dp ; // = 1.0/power;
double m; // float m = @power*dTwoPi
double dp2 ;  //float dp2 = 1.0/(2.0+@power) ; Changed from originally trying 1/power




       
// memmory virtual 1D array
unsigned char *data;       
size_t MemmorySize;       
   
   
void GiveLinasColor(double position , int k, unsigned char c[])
{
  /* based on the code by Linas Vepstas January 16 1994 : void make_cmap (void) */

   
  int i;
  int iMax = 240;
  i=(int)(iMax-1)*position; 
  c[0] = c[1] = c[2] = 0;
  /* set up a default look up table */
  /* ramp up to blue */
  if (i<60) {
    c[k] = 0;
    c[k+1] = 0;
    c[k+2] = (unsigned char) i*3;
  }
  /* ramp down from blue, up to green */
  if (i>=60 && i<120) {
    c[k] = 0;
    c[k+1] = (unsigned char) (i-60)*3;
    c[k+2] = (unsigned char) (120-i)*3;
  }
  /* ramp from green to yellow */
  if (i>=120 && i<180) {
    /* vlt[i].r = (char) (((i-120)*7) / 2); */
    c[k] = (unsigned char) (210 - (7*(180-i)*(180-i)) / 120);
    c[k+1] = (unsigned char) (210 -i/4);
    c[k+2] = 0;
  }
  /* ramp from yellow to red (pink) */
  if (i>=180 && i<iMax) {
    c[k] = (unsigned char) (210 + (3*(i-180))/4);
    c[k+1] = (unsigned char) (510 - 2*i);
    c[k+2] = (unsigned char) (i-180)/3;
  }
   
}

       
       
/*
   gives position ( index) in 1D virtual array  of 2D point (iX,iY) from ; uses also global variable iWidth
   without bounds check !!
*/
int f(int ix, int iy)
{ return ColorBytes*(ix + iy*iWidth); }
       
       
       
double complex give_c(int iX, int iY){
  double Cx,Cy;
  Cy=CyMax - iY*PixelHeight; // inverse y axis
 
  Cx=CxMin + iX*PixelWidth;
   
  return Cx+Cy*I;
 
 
}
 
 

double GiveTurn(double complex z)
{
  double argument;
 
  argument = carg(z); //   argument in radians from -pi to pi
  if (argument<0) argument=argument + TwoPi; //   argument in radians from 0 to 2*pi
  return argument/TwoPi ; // argument in turns from 0.0 to 1.0
}



double GiveAtan2(complex double Z){
return atan2(cimag(Z), creal(Z)); // double atan2( double y, double x );
}



double Give_t(complex double zv[], int final_i){

double a = 0.0;
double b = 0.0;
double c = 0.0; // case sensitive
double s = 0.0;
int i = final_i;

complex double Z;

Z = zv[i];
s = GiveAtan2(Z);
if (s<0) s = s + TwoPi;

while (i > skipiter){
         
    i = i - 1;
    Z = zv[i];
    a = GiveAtan2(Z);
    if (version==0)
    b = fmod(power*a, TwoPi);
    else b = GiveAtan2(cpow(Z,power));
   
    if (b < 0.0) b = b + TwoPi;
        c = c + b - s;
   
    if (c >= M_PI)
    {s = s + TwoPi;
   
    if (fixit || accident) c = c - TwoPi;}
   
          else
          {
          if (c < -M_PI)
          {s = s - TwoPi;
          if (fixit || accident) c = c + TwoPi;}}
     
        if (fixit) c = dp2*c;
    else {if (!accident) c = 0.0;}
   
        if (a < 0.0) a = a + TwoPi;
        s = dp*(s + TwoPi*floor(a*m));
  } // end while
 
 
  return s;
 



}



int ComputeAndSavePixelColor(int iX, int iY){


complex double zv[IterationMax];
complex double C; //
  int i; // iteration
  double complex Z= 0.0; // initial value for iteration Z0
  int k; // index of the 1D array
   
  C = give_c(iX, iY);
   
  // iteration = loop
  for(i=0;i<IterationMax;i++)
    {
     
      Z=Z*Z+C;
      zv[i] = Z;
      if(cabs(Z)>EscapeRadius) break; // bailout test
    }
   
     
  // index of 1D memory array
  k = f(iX, iY); 
   
 
  // final
  double t = Give_t(zv, i);
  GiveLinasColor(t , k,  data);
     
   
 
   
return 0;
}
 
 
 
int setup(){

  dp = 1.0/power;
  m  = power*dTwoPi;
  dp2 = 1.0/(2.0+power) ; //Changed from originally trying 1/power
 
  //
  PixelWidth=(CxMax-CxMin)/iWidth;
  PixelHeight=(CyMax-CyMin)/iHeight;
  //
   
  //
  MemmorySize = iWidth * iHeight * ColorBytes * sizeof (unsigned char); // https://stackoverflow.com/questions/492384/how-to-find-the-sizeof-a-pointer-pointing-to-an-array
       
  /* create dynamic 1D arrays for colors   */
  data = malloc (MemmorySize);
  if (data == NULL )
    { fprintf (stderr, " Error: Could not allocate memory !!! \n"); return 1; }

  printf (" No errors. End of setup \n");
  return 0;

}
 
 
 
 
 
// save dynamic "A" array to pgm file
int SaveArray_2_PPM_file (unsigned char A[])
{

  FILE *fp;
  const unsigned int MaxColorComponentValue = 255; /* color component is coded from 0 to 255 ;  it is 8 bit color file */
  char *filename = "fieldlines.ppm";
  char *comment = "# "; /* comment should start with # */

  /* save image to the pgm file  */
  fp = fopen (filename, "wb"); /*create new file,give it a name and open it in binary mode  */
  fprintf (fp, "P6\n %s\n %u %u\n %u\n", comment, iWidth, iHeight, MaxColorComponentValue); /*write header to the file */
  fwrite (A, MemmorySize, 1, fp); /*write A array to the file in one step */
  printf ("File %s saved. \n", filename);
  fclose (fp);

  return 0;
}


 
 
void CreateImage(){
  int iX,iY; // screen = integer coordinate in pixels       

  // fill the array = render image = scanline 2D  of virtual 2D array
  for(iY=0;iY<iHeight;iY++)
    for(iX=0;iX<iWidth;iX++)
      ComputeAndSavePixelColor(iX, iY);
     
     
  SaveArray_2_PPM_file (data);       
}
 
 
 
void info(){

  printf(" Parameter plane ( c plane) with Mandelbrot set for complex quadratic polynomial fc(z) = z^2 + c\n ");
  printf(" Rectangle part of 2D parameter plane: corners: \n CxMin = %f;   CxMax = %f;  CyMin = %f; CyMax = %f \n ", CxMin, CxMax, CyMin, CyMax);
  printf(" center and radius: \n CenterX = %f;   CenterY = %f;  radius = %f\n ", (CxMax+CxMin)/2.0, (CyMax+CyMin)/2.0, fabs(CyMax-CyMin)/2.0);
  printf(" Mag = zoom = %f\n ",  2.0/fabs(CyMax-CyMin));
  printf("PixelWidth = %f and PixelHeight =%f\n", PixelWidth, PixelHeight);
  printf(" Escape Radius = %f\n ", EscapeRadius);
  printf(" Iteration Max = %d\n ", IterationMax);



}
 
 
 
void close(){
 
  info();
  free(data);
}
 
 
 
int main()
{
 
  setup();     
  CreateImage();     
  close();
 
       
  return 0;
}


It copiles without error anf makes image , but image is bad .

Where is the error ?

TIA

Adam

Offline claude

  • *
  • Fractal Freak
  • **
  • Posts: 714
    • mathr.co.uk
« Reply #1 on: January 04, 2018, 04:52:33 PM »
I haven't tested this (I don't know what it should look like) but from comparing the two listings you are missing an update to the 'a' variable, hidden inside an if test:

Quote
Code: [Select]
        if (fixit) c = dp2*c;
    else {if (!accident) c = 0.0;}
   
        if (a < 0.0) a = a + TwoPi;
        s = dp*(s + TwoPi*floor(a*m));
  } // end while

Should be:

Quote
Code: [Select]
        if (fixit) c = dp2*c;
    else {if (!accident) c = 0.0;}
   
        if ((a = carg(zv[i]))< 0.0) a = a + TwoPi;
        s = dp*(s + TwoPi*floor(a*m));
  } // end while

Hopefully this is the only error.

Offline Adam Majewski

  • *
  • Strange Attractor
  • ******
  • Posts: 93
« Reply #2 on: January 04, 2018, 05:01:49 PM »
It shoud look like the externa angle, see attached ea_uf.png
but still looks : fl.png

to be true I do not see any updates of :
int version = 1;
int julia = 1;
int skipiter = 0;
int fixskip = 0;
int fixit = 0;
int accident = 0;

variables
« Last Edit: January 04, 2018, 05:18:19 PM by Adam Majewski »

Offline claude

  • *
  • Fractal Freak
  • **
  • Posts: 714
    • mathr.co.uk
« Reply #3 on: January 04, 2018, 10:41:29 PM »
Another issue is mismatch between s (based on 2PI) and givecolor (based on 0..1).  Replace

Code: [Select]
  } // end while
 
  return s;
 
}

with

Quote
Code: [Select]
  } // end while
 
  s /= TwoPi;
  s -= floor(s);
  return s;
 
}

You may want a larger escape radius for smoother appearance.

Offline Adam Majewski

  • *
  • Strange Attractor
  • ******
  • Posts: 93
« Reply #4 on: January 05, 2018, 05:39:20 PM »
Thx Claude. It works now

Is it a different method then :
https://gitlab.com/adammajewski/parameter_external_angle
?

Code: [Select]

/*


   http://linas.org/art-gallery/escape/phase/phase.html
   image
   http://linas.org/art-gallery/escape/phase/phase.gif
 
   c console program:
   
   
   --------------------------------
   1. draws Mandelbrot set for Fc(z)=z*z +c
   using Mandelbrot algorithm ( boolean escape time )
   -------------------------------         
   2. technique of creating ppm file is  based on the code of Claudio Rocchini
   http://en.wikipedia.org/wiki/Image:Color_complex_plot.jpg
   create 24 bit color graphic file ,  portable pixmap file = PPM
   see http://en.wikipedia.org/wiki/Portable_pixmap
   to see the file use external application ( graphic viewer)
   
   
   complex point c -> virtual 2D array -> memory 1D array -> ppm file on the disc -> png file
   
   C -> pixel (iX,iY)  -> index k  -> 24bit color
   
   -----
   https://stackoverflow.com/questions/6418807/how-to-work-with-complex-numbers-in-c
   complex numbers are built in type
 
 --------------
 formated with emacs
   -------------
   to compile :

 
 
   gcc f.c -lm -Wall
 
 
   ./a.out
   
   
   to convert to png using ImageMagic

   convert phase.ppm phase.png 



   ----------------------
   git add  phase.png phase.c README.md
   git commit -m "phase"
   git push -u origin master



 
 
*/
#include <stdio.h>
#include <stdlib.h> // malloc
#include <math.h>
#include <complex.h> // https://stackoverflow.com/questions/6418807/how-to-work-with-complex-numbers-in-c
 


 
/* screen ( integer) coordinate */

const int iWidth  = 1000;
const int iHeight = 1000;


/* world ( double) coordinate = parameter plane*/
// double complex C =  Cx + Cy*I ;
const double CxMin=-2.4;
const double CxMax= 1.5;
const double CyMin=-1.95;
const double CyMax= 1.95;

/* */
double PixelWidth; //=(CxMax-CxMin)/iWidth;
double PixelHeight; // =(CyMax-CyMin)/iHeight;


/* color component ( R or G or B) is coded from 0 to 255 */
/* it is 24 bit color RGB file */
int ColorBytes = 3; // 3*8 = 24 bit color


#define iBlack  0
#define iWhite  255
#define iRed 1
#define iGreen  2   
#define iBlue 
#define iCyan  4   

       

/* iterations  */
const int IterationMax=4000;

/* bail-out value , radius of circle ;  */
const double EscapeRadius=100.0;

/*

param version     caption = "Version"     default = true     hint = ""   endparam

param power     caption = "Degree of divergence"     default = 2.0     hint = "This is an estimate of the divergence of the main formula, e.g. 2 for z^2+c, 3 for z^3+c etc."   endparam

param julia     caption = "Julia ?"     default = false     hint = "Enable for correct rendering of Julia Sets."   endparam

param skipiter     caption = "'Start' iteration"     default = 0     min = 0     hint = ""   endparam

param fixskip     caption = "Fix 'Start' iteration"      default = false    hint = ""   endparam


param fixit     caption = "Fix Field Lines"     default = false     hint = ""   endparam

param accident     caption = "Happy Accident"     default = false     hint = ""     visible = !@fixit   endparam

*/

double power = 2.0;
int version = 1;
int julia = 1;
int skipiter = 0;
int fixskip = 0;
int fixit = 0;
int accident = 0;

/*

 
global:
  float TwoPi = 2.0*#pi
  float dTwoPi = 0.5/#pi
  float dp = 1.0/@power
  float m = @power*dTwoPi
  float dp2 = 1.0/(2.0+@power) ; Changed from originally trying 1/power


 

*/

double TwoPi=2.0*M_PI; // float TwoPi = 2.0*#pi
double dTwoPi = 0.5/M_PI;  // float dTwoPi = 0.5/#pi
double dp ; // = 1.0/power;
double m; // float m = @power*dTwoPi
double dp2 ;  //float dp2 = 1.0/(2.0+@power) ; Changed from originally trying 1/power




       
// memmory virtual 1D array
unsigned char *data;       
size_t MemmorySize;   



void ColorPixel(int iColor, int k, unsigned char c[])
{
  switch (iColor)
    {
    case iBlack : c[k]   = 0; c[k+1] = 0; c[k+2] = 0; break;
    case iRed:    c[k]   = 255; c[k+1] = 0; c[k+2] = 0;   break;
    case iGreen : c[k]   = 0; c[k+1] = 255; c[k+2] = 0; break;
    case iBlue : c[k]   = 0; c[k+1] = 0; c[k+2] = 255; break;
    case iCyan : c[k]   = 0; c[k+1] = 255; c[k+2] = 255; break;
    case iWhite : c[k]   = 255; c[k+1] = 255; c[k+2] = 255; break;
    }
}
   
   
   
void GiveLinasColor(double position , int k, unsigned char c[])
{
  /* based on the code by Linas Vepstas January 16 1994 : void make_cmap (void) */

   
  int i;
  int iMax = 240;
  i=(int)(iMax-1)*position; 
  c[0] = c[1] = c[2] = 0;
  /* set up a default look up table */
  /* ramp up to blue */
  if (i<60) {
    c[k] = 0;
    c[k+1] = 0;
    c[k+2] = (unsigned char) i*3;
  }
  /* ramp down from blue, up to green */
  if (i>=60 && i<120) {
    c[k] = 0;
    c[k+1] = (unsigned char) (i-60)*3;
    c[k+2] = (unsigned char) (120-i)*3;
  }
  /* ramp from green to yellow */
  if (i>=120 && i<180) {
    /* vlt[i].r = (char) (((i-120)*7) / 2); */
    c[k] = (unsigned char) (210 - (7*(180-i)*(180-i)) / 120);
    c[k+1] = (unsigned char) (210 -i/4);
    c[k+2] = 0;
  }
  /* ramp from yellow to red (pink) */
  if (i>=180 && i<iMax) {
    c[k] = (unsigned char) (210 + (3*(i-180))/4);
    c[k+1] = (unsigned char) (510 - 2*i);
    c[k+2] = (unsigned char) (i-180)/3;
  }
   
}

       
       
/*
   gives position ( index) in 1D virtual array  of 2D point (iX,iY) from ; uses also global variable iWidth
   without bounds check !!
*/
int f(int ix, int iy)
{ return ColorBytes*(ix + iy*iWidth); }
       
       
       
double complex give_c(int iX, int iY){
  double Cx,Cy;
  Cy=CyMax - iY*PixelHeight; // inverse y axis
 
  Cx=CxMin + iX*PixelWidth;
   
  return Cx+Cy*I;
 
 
}
 
 

double GiveTurn(double complex z)
{
  double argument;
 
  argument = carg(z); //   argument in radians from -pi to pi
  if (argument<0) argument=argument + TwoPi; //   argument in radians from 0 to 2*pi
  return argument/TwoPi ; // argument in turns from 0.0 to 1.0
}



double GiveAtan2(complex double Z){
return atan2(cimag(Z), creal(Z)); // double atan2( double y, double x );
}



double Give_t(complex double zv[], int final_i){

double a = 0.0;
double b = 0.0;
double c = 0.0; // case sensitive
double s = 0.0;
int i = final_i;

complex double Z;

Z = zv[i];
s = GiveAtan2(Z);
if (s<0) s = s + TwoPi;

while (i > skipiter){
         
    i = i - 1;
    Z = zv[i];
    a = GiveAtan2(Z);
    if (version==0)
    b = fmod(power*a, TwoPi);
    else b = GiveAtan2(cpow(Z,power));
   
    if (b < 0.0) b = b + TwoPi;
        c = c + b - s;
   
    if (c >= M_PI)
    {s = s + TwoPi;
   
    if (fixit || accident) c = c - TwoPi;}
   
          else
          {
          if (c < -M_PI)
          {s = s - TwoPi;
          if (fixit || accident) c = c + TwoPi;}}
     
        if (fixit) c = dp2*c;
    else {if (!accident) c = 0.0;}
   
        //if (a < 0.0) a = a + TwoPi;
        //s = dp*(s + TwoPi*floor(a*m));
                if ((a = carg(zv[i]))< 0.0) a = a + TwoPi;
        s = dp*(s + TwoPi*floor(a*m));

       
  } // end while
 
  s /= TwoPi;
  s -= floor(s);
 
  return s;
 



}



int ComputeAndSavePixelColor(int iX, int iY){


complex double zv[IterationMax];
complex double C; //
  int i; // iteration
  double complex Z= 0.0; // initial value for iteration Z0
  int k; // index of the 1D array
   
  C = give_c(iX, iY);
   
  // iteration = loop
  for(i=0;i<IterationMax;i++)
    {
     
      Z=Z*Z+C;
      zv[i] = Z;
      if(cabs(Z)>EscapeRadius) break; // bailout test
    }
   
     
  // index of 1D memory array
  k = f(iX, iY); 
   
 
  // final
  if (i==IterationMax)
  ColorPixel( iBlack, k, data);
  else {
 
  double t = Give_t(zv, i);
  GiveLinasColor(t , k,  data);
     
 
  }
 
   
 
   
return 0;
}
 
 
 
int setup(){

  dp = 1.0/power;
  m  = power*dTwoPi;
  dp2 = 1.0/(2.0+power) ; //Changed from originally trying 1/power
 
  //
  PixelWidth=(CxMax-CxMin)/iWidth;
  PixelHeight=(CyMax-CyMin)/iHeight;
  //
   
  //
  MemmorySize = iWidth * iHeight * ColorBytes * sizeof (unsigned char); // https://stackoverflow.com/questions/492384/how-to-find-the-sizeof-a-pointer-pointing-to-an-array
       
  /* create dynamic 1D arrays for colors   */
  data = malloc (MemmorySize);
  if (data == NULL )
    { fprintf (stderr, " Error: Could not allocate memory !!! \n"); return 1; }

  printf (" No errors. End of setup \n");
  return 0;

}
 
 
 
 
 
// save dynamic "A" array to pgm file
int SaveArray_2_PPM_file (unsigned char A[])
{

  FILE *fp;
  const unsigned int MaxColorComponentValue = 255; /* color component is coded from 0 to 255 ;  it is 8 bit color file */
  char *filename = "fieldlines_4000.ppm";
  char *comment = "# "; /* comment should start with # */

  /* save image to the pgm file  */
  fp = fopen (filename, "wb"); /*create new file,give it a name and open it in binary mode  */
  fprintf (fp, "P6\n %s\n %u %u\n %u\n", comment, iWidth, iHeight, MaxColorComponentValue); /*write header to the file */
  fwrite (A, MemmorySize, 1, fp); /*write A array to the file in one step */
  printf ("File %s saved. \n", filename);
  fclose (fp);

  return 0;
}


 
 
void CreateImage(){
  int iX,iY; // screen = integer coordinate in pixels       

  // fill the array = render image = scanline 2D  of virtual 2D array
  for(iY=0;iY<iHeight;iY++)
    for(iX=0;iX<iWidth;iX++)
      ComputeAndSavePixelColor(iX, iY);
     
     
  SaveArray_2_PPM_file (data);       
}
 
 
 
void info(){

  printf(" Parameter plane ( c plane) with Mandelbrot set for complex quadratic polynomial fc(z) = z^2 + c\n ");
  printf(" Rectangle part of 2D parameter plane: corners: \n CxMin = %f;   CxMax = %f;  CyMin = %f; CyMax = %f \n ", CxMin, CxMax, CyMin, CyMax);
  printf(" center and radius: \n CenterX = %f;   CenterY = %f;  radius = %f\n ", (CxMax+CxMin)/2.0, (CyMax+CyMin)/2.0, fabs(CyMax-CyMin)/2.0);
  printf(" Mag = zoom = %f\n ",  2.0/fabs(CyMax-CyMin));
  printf("PixelWidth = %f and PixelHeight =%f\n", PixelWidth, PixelHeight);
  printf(" Escape Radius = %f\n ", EscapeRadius);
  printf(" Iteration Max = %d\n ", IterationMax);



}
 
 
 
void close(){
 
  info();
  free(data);
}
 
 
 
int main()
{
 
  setup();     
  CreateImage();     
  close();
 
       
  return 0;
}


clip
Can Ultrafractal orbit trap animated images?

Started by abbaddon1001 on UltraFractal

0 Replies
189 Views
Last post May 22, 2018, 10:32:35 PM
by abbaddon1001
xx
Triangle Inequality Average Algorithm

Started by mikoval on Fractal Mathematics And New Theories

10 Replies
339 Views
Last post August 05, 2018, 02:36:30 PM
by FractalDave
xx
Chaos to speed up numerical algorithm

Started by gerrit on Physics

1 Replies
254 Views
Last post February 15, 2018, 12:40:28 PM
by knighty
xx
Understanding the Algorithm for Sierpinski Tetrahedron Distance

Started by Ebanflo on Fractal Mathematics And New Theories

1 Replies
126 Views
Last post March 12, 2018, 11:07:47 AM
by knighty
xx
modifying Mandelbox algorithm with shape inversions in Fragmentarium

Started by CozyG on Fragmentarium

1 Replies
70 Views
Last post September 09, 2018, 02:21:17 AM
by mclarekin