implemented julia sets

4 files changed, 55 insertions, 8 deletions

diff --git a/src/fractals.c b/src/fractals.c
index a2a43db..2c5b42c 100644
--- a/src/fractals.c
+++ b/src/fractals.c
@@ -1,6 +1,3 @@
-#include <complex.h>
-#include <stddef.h>
-#include <stdint.h>
 #include <stdio.h>
 #include <unistd.h>
 #include <stdlib.h>
@@ -71,8 +68,23 @@ int main(const int argc, char *argv[]) {
     }
 
 
-    mandelbrot_grid(grid, iterations);
+    //const long double complex c = 0.285L + 0.01L*I;
+    const long double complex c = -0.835L -0.321L* I;
+    const double R = 100;
+    julia_grid(grid, iterations, c, R);
 
     if(verbose)print_grid_info(grid);
     print_grid(grid, iterations);
+
+    // //write grid to file
+    // FILE* write_file = fopen("test.grid", "wb");
+    // write_grid(write_file , grid);
+    // fclose(write_file);
+    //
+    // //attempt to read grid from file
+    // FILE* read_file = fopen("test2.grid", "rb");
+    // grid_t* grid2 = read_grid(read_file);
+    // fclose(read_file);
+    //
+    // printf("Grids are %s equal\n", grid_equal(grid, grid2) ? "exactly" :"not exactly");
 }
diff --git a/src/fractals.h b/src/fractals.h
index 5506745..d27eb00 100644
--- a/src/fractals.h
+++ b/src/fractals.h
@@ -11,4 +11,5 @@ void mandelbrot_grid(grid_t* grid,  const size_t max_iterations);
 size_t multibrot(const long double complex z0, const size_t max_iterations, const double d);
 void multibrot_grid(grid_t* grid, const size_t max_iterations, const double d);
 
-size_t julia(const double R, const long double complex z0, const long double complex c, const size_t max_iterations);
+size_t julia(const long double complex z0, const long double complex c, const size_t max_iterations, const double R);
+void julia_grid(grid_t* grid, const size_t max_iterations, const long double complex c, const double R);
diff --git a/src/serial-fractals.c b/src/serial-fractals.c
index 718fccb..fe420d6 100644
--- a/src/serial-fractals.c
+++ b/src/serial-fractals.c
@@ -59,10 +59,10 @@ void multibrot_grid(grid_t* grid, const size_t max_iterations, const double d){
 /*
  * Computes ????? for a julia set
  * implementation of https://en.wikipedia.org/wiki/Julia_set#Pseudocode
+ *
+ * This behaves weirdly, needs a very small number of iterations to be visibile
  */
-size_t julia(const double R, const long double complex z0, const long double complex c, const size_t max_iterations){
-    //FIXME: I'm notsure if this is currently implemented correctly
-    if(R*R - R >= cabsl(z0)) return 0;
+size_t julia(const long double complex z0, const long double complex c, const size_t max_iterations, const double R){
     long double complex z = z0;
 
     size_t iteration = 0;
@@ -73,3 +73,10 @@ size_t julia(const double R, const long double complex z0, const long double com
     return iteration;
 }
 
+void julia_grid(grid_t* grid, const size_t max_iterations, const long double complex c, const double R){
+    const size_t size = grid->size;
+    size_t* data = grid->data;
+    for(size_t i = 0; i <size; i++){
+        data[i] = julia(grid_to_complex(grid, i), c, max_iterations, R);
+    }
+}
diff --git a/src/shared-fractals.c b/src/shared-fractals.c
index 9f65e80..49a034e 100644
--- a/src/shared-fractals.c
+++ b/src/shared-fractals.c
@@ -60,3 +60,30 @@ void multibrot_grid(grid_t* grid, const size_t max_iterations, const double d){
         data[i] = multibrot(grid_to_complex(grid, i), max_iterations, d);
     }
 }
+
+/*
+ * Computes ????? for a julia set
+ * implementation of https://en.wikipedia.org/wiki/Julia_set#Pseudocode
+ *
+ * This behaves weirdly, needs a very small number of iterations to be visibile
+ */
+
+size_t julia(const long double complex z0, const long double complex c, const size_t max_iterations, const double R){
+    long double complex z = z0;
+
+    size_t iteration = 0;
+    while(cabsl(z) < R && iteration < max_iterations){
+        z = z * z + c;
+        iteration++;
+    }
+    return iteration;
+}
+
+void julia_grid(grid_t* grid, const size_t max_iterations, const long double complex c, const double R){
+    const size_t size = grid->size;
+    size_t* data = grid->data;
+    #pragma omp parallel for default(none) shared(data, size, grid, max_iterations, c, R) schedule(dynamic)
+    for(size_t i = 0; i <size; i++){
+        data[i] = julia(grid_to_complex(grid, i), c, max_iterations, R);
+    }
+}