fractal/julia.c

/*
 * fractal
 * Written in 2020 by Lucas
 * CC0 1.0 Universal/Public domain - No rights reserved
 *
 * To the extent possible under law, the author(s) have dedicated all
 * copyright and related and neighboring rights to this software to the
 * public domain worldwide. This software is distributed without any
 * warranty. You should have received a copy of the CC0 Public Domain
 * Dedication along with this software. If not, see
 * <http://creativecommons.org/publicdomain/zero/1.0/>.
 */
/* gepopt visibility for GLIBC */
#if !defined(__OpenBSD__)
#define _DEFAULT_SOURCE
#endif

#include <complex.h>
#include <limits.h>
#include <math.h>
#include <stdlib.h>
#include <unistd.h>

#include "criticals.h"
#include "palette.h"
#include "pnmout.h"
#include "util.h"

#include "fun.h"

#define RITNUM(c, i) (c[i].it - c[i].rit)

#define OPTSTRING "a:B:b:c:D:d:e:h:Ll:N:n:w:x:y:z:"
#define RADIAN 3.14159265358979323846L

struct color contour_color = { 0 };
static long double bailout = 1e100,
    epsilon = 1e-8,
    sec_w = 4.0,
    contour = 0.0005L;
static long niters = 250, newton_iters = 100;

static void
print_critical_points(void)
{
	long double complex z;
	size_t i, n;

	n = criticals_size();
	fprintf(stderr, "Critical points:\n");
	for (i = 0; i < n; i++) {
		z = criticals_get(i);
		fprintf(stderr, "%4zu. (%.20Lf, %.20Lf)\n", i, creall(z),
		    cimagl(z));
	}
}

static int
tends_to_critical(long double complex z, long double *t)
{
	long double d;
	size_t i, n;

	n = criticals_size();
	for (i = 0; i < n; i++) {
		d = cabsl(z - criticals_get(i));
		if (d < epsilon) {
			*t = d;
			return 1;
		}
	}

	return 0;
}

static long
iterate(long double complex z, long double complex c, long double *rit)
{
	long double complex dz = 1.0L;
	long double t;
	long i;

	for (i = 0; i < niters && cabsl(z) < bailout; i++) {
		dz = df(z) * dz;
		z = f(z) + c;

		if (tends_to_critical(z, &t)) {
			*rit = logl(logl(t) / logl(epsilon)) / logl(2.0L);
			return fabsl(logl(t)) * t < sec_w * contour * cabsl(dz)
			    ? -1 : i;
		}
	}

	t = cabsl(z);
	if (t >= bailout) {
		*rit = logl(logl(t) / logl(bailout)) / logl(2.0L);
		return fabsl(logl(t)) * t < sec_w * contour * cabsl(dz)
		    ? -1 : i;
	}

	*rit = 0;
	return -1;
}

static void
usage(void)
{
	fprintf(stderr, "Usage: %s julia [-L] [-a center_x] [-B bailout]\n"
	    "    [-b center_y] [-c contour] [-D displacement] [-d density]\n"
	    "    [-e epsilon] [-h height] [-l light_intensity]\n"
	    "    [-N newton_iters] [-n iters] [-w width] [-x julia_x]\n"
	    "    [-y julia_y] [-z zoom]\n",
	    xgetprogname());
	exit(1);
}

int
julia_main(int argc, char *argv[])
{
	struct color *palette;
	struct {
		long long it;
		long double rit;
	} *rescache;
	struct color color;
	long double complex crit, julia, z;
	long double julia_x = 0.0,
	    julia_y = 0.0,
	    cen_x = 0.0,
	    cen_y = 0.0,
	    density = 28.3L,
	    light_intensity = 84.3L,
	    displacement = 0,
	    zoom = 1.0,
	    light, light_x, light_y, light_z, rit, rn, sec_h, zx, zy;
	long long it, idx, ritnum;
	int palette_size = 1000,
	    width = 640,
	    height = 400,
	    rescache_width, rescache_height, i, j;
	int ch, do_light = 0;

	while ((ch = getopt(argc, argv, OPTSTRING)) != -1)
		switch (ch) {
		case 'a':
			cen_x = parse_double(optarg);
			break;
		case 'B':
			bailout = parse_double(optarg);
			break;
		case 'b':
			cen_y = parse_double(optarg);
			break;
		case 'c':
			contour = parse_double(optarg);
			break;
		case 'D':
			displacement = parse_double(optarg);
			break;
		case 'd':
			density = parse_double(optarg);
			break;
		case 'e':
			epsilon = parse_double(optarg);
			break;
		case 'h':
			height = parse_integer(optarg, 1, INT_MAX);
			break;
		case 'L':
			do_light = 1;
			break;
		case 'l':
			light_intensity = parse_double(optarg);
			break;
		case 'N':
			newton_iters = parse_integer(optarg, 1, LONG_MAX);
			break;
		case 'n':
			niters = parse_integer(optarg, 1, LONG_MAX);
			break;
		case 'w':
			width = parse_integer(optarg, 1, INT_MAX);
			break;
		case 'x':
			julia_x = parse_double(optarg);
			break;
		case 'y':
			julia_y = parse_double(optarg);
			break;
		case 'z':
			zoom = parse_double(optarg);
			break;
		default:
			usage();
		}
	argc -= optind;
	argv += optind;
	if (argc > 1)
		usage();

	palette = palette_gen(palette_size);

	rescache_width = width + (do_light ? 1 : 0);
	rescache_height = height + (do_light ? 1 : 0);

	if ((rescache = malloc(sizeof(*rescache) * rescache_width
	    * rescache_height)) == NULL)
		errx(1, "out of memory");

	if (do_light) {
		long double angle_xy, angle_z;
		angle_xy = drand48() * 2.0L * RADIAN;
		do
			angle_z = drand48();
		while (angle_z == 0.0L || angle_z == 1.0L);
		angle_z = angle_z * 0.5L * RADIAN;
		light_x = cosl(angle_xy) * sinl(angle_z);
		light_y = sinl(angle_xy) * sinl(angle_z);
		light_z = cosl(angle_z);
	}

	sec_h = sec_w * height / (long double)width;
	julia = julia_x + I * julia_y;

	for (i = 0; i < height; i++) {
		zy = cen_y + (0.5 - i / (long double)height) * sec_h;
		for (j = 0; j < width; j++) {
			zx = cen_x + (j / (long double)width - 0.5L) * sec_w;
			z = zx + I * zy;
			if (criticals_test(df, ddf, z, &crit, newton_iters,
			    bailout, epsilon))
				criticals_add(crit);
		}
	}

	print_critical_points();

	sec_w /= zoom; sec_h /= zoom;
	for (i = 0; i < rescache_height; i++) {
		zy = cen_y + (0.5L - i / (long double)height) * sec_h;
		for (j = 0; j < rescache_width; j++) {
			zx = cen_x + (j / (long double)width - 0.5L) * sec_w;
			z = zx + I * zy;
			it = iterate(z, julia, &rit);
			rescache[i * rescache_width + j].it = it;
			rescache[i * rescache_width + j].rit = rit;
		}
	}

	pnmout_header(stdout, width, height);
	for (i = 0; i < height; i++) {
		for (j = 0; j < width; j++) {
			rn = RITNUM(rescache, i * rescache_width + j);

			if (rescache[i * rescache_width + j].it < 0)
				color = contour_color;
			else {
				ritnum = rn * density;
				if (do_light) {
					long double lxrn, lyrn, l, lx, ly, lz;
					lxrn = RITNUM(rescache,
					    i * rescache_width + j + 1);
					lyrn = RITNUM(rescache,
					    (i + 1) * rescache_width + j);
					lx = -sec_h * (lxrn - rn);
					ly = sec_w * (lyrn - rn);
					lz = sec_w * sec_h;
					l = sqrtl(lx * lx + ly * ly + lz * lz);
					light = (light_x * lx + light_y * ly
					    + light_z * lz) / l;
					ritnum += light * light_intensity;
				}
				idx = (long long)ritnum
				    + 1000000 * palette_size + displacement;
				color = palette[idx % palette_size
				    + (idx < 0 ? palette_size : 0)];
			}

			pnmout_pixel(stdout, color);
		}
	}

	free(rescache);
	free(palette);

	return !(fflush(stdout) == 0);
}