/* Visualisation of Mandelbrot Set using OpenMP
 * $Id: mb.c,v 1.3 2024/07/04 10:45:19 oj14ozun Exp $
 * https://wwwcip.cs.fau.de/~oj14ozun/src+etc/mb.c */

#include <assert.h>
#include <complex.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#ifndef RESOLUTION
#define RESOLUTION 1000
#endif
#define XSTART -2.0
#define XEND 1.0
#define YSTART -1.0
#define YEND 1.0
#define DEPTH 255
#define BITS 255

#define W  ((XEND-XSTART) * RESOLUTION)
#define Wu ((unsigned) W)
#define H  ((YEND-YSTART) * RESOLUTION)
#define Hu ((unsigned) H)

static unsigned mb(double complex c)
{
     unsigned n = 0;
     double complex d = 0, d1;
     do {
	  d1 = d;
	  d = d * d + c;
	  n++;
     } while (cabs(d-d1) > 1e-8 && n < DEPTH);
     return n;
}

static void compute(unsigned char matrix[])
{
#pragma omp parallel for shared(matrix) schedule(static)
     for (unsigned i = 0; i < Wu; ++i) {
	  unsigned n = i * Hu * 3;
	  for (unsigned j = 0; j < Hu; ++j) {
	       double complex c = 0;
	       c += (XSTART + i * 1./RESOLUTION);
	       c += (YSTART + j * 1./RESOLUTION) * I;

	       double v = (mb(c) / (double)DEPTH);
	       matrix[n++] = BITS * v * sinl(M_PI_2 * i/W);
	       matrix[n++] = BITS * v * erfl(M_PI_2 * j/H);
	       matrix[n++] = BITS * v * cosl(M_PI_2 * (j+i)/(W+H));
          }
     }
}

int main()
{
     unsigned char *matrix = calloc(Hu * Wu, 3);
     if (!matrix) {
	  perror("calloc");
	  return EXIT_FAILURE;
     }

     compute(matrix);
     printf("P6 %d %d %d\n", Hu, Wu, BITS);
     fwrite(matrix, H*W, 3, stdout);
     
     return EXIT_SUCCESS;
}

/* Local Variables: */
/* compile-command: "cc -o mb  mb.c -Ofast -lm -fopenmp" */
/* End: */