/* mandel.c - Mandelbrot image generator
* Written 6 Nov 2008 by John T. Wodder II
* Last edited 6 Nov 2008 by John Wodder
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <unistd.h>

_Bool escapes(double a, double b);

int main(int argc, char** argv) {
double minx = -2.0, maxx = 0.5, miny = -1.0, maxy = 1.0;
int precision = 100, ch;
while ((ch = getopt(argc, argv, "x:y:X:Y:p:")) != -1) {
 switch (ch) {
  case 'x': minx = strtod(optarg, NULL); break;
  case 'y': miny = strtod(optarg, NULL); break;
  case 'X': maxx = strtod(optarg, NULL); break;
  case 'Y': maxy = strtod(optarg, NULL); break;
  case 'p': precision = strtol(optarg, NULL, 0); break;
  default:
   fprintf(stderr, "Usage: mandel [-x minx] [-y miny] [-X maxx] [-Y maxy]"
    " [-p precision] [outfile]\n");
   return 2;
 }
}
char* filename = (optind < argc) ? argv[optind] : "mandelbrot.eps";
FILE* eps = fopen(filename, "w");
if (eps == NULL) {
 fprintf(stderr, "mandel: %s: ", filename);
 perror(NULL);
 return 5;
}
fprintf(eps, "%%!PS-Adobe-3.0 EPSF-3.0\n%%%%BoundingBox: %.1f %ld %.1f %ld\n",
 minx * precision - 0.5, lround(miny * precision), maxx * precision + 0.5,
 lround(maxy * precision + 1));
for (double x = minx; x <= maxx; x += 1.0 / precision) {
 for (double y = miny; y <= maxy; y += 1.0 / precision) {
  if (!escapes(x, y)) {
   fprintf(eps, "%ld %ld moveto 0 1 rlineto\n", lround(x * precision),
    lround(y * precision));
  }
 }
}
fprintf(eps, "stroke showpage\n");
fclose(eps);
return 0;
}

_Bool escapes(double a, double b) {
if (hypot(a, b) > 2.0) return 1;
double x = 0.0, y = 0.0;
for (int i=0; i<100; i++) {
 double xTmp = x;
 x = x*x - y*y + a;
 y = 2 * xTmp * y + b;
 if (hypot(x, y) > 2.0) return 1;
}
return 0;
}