278 lines
6.3 KiB
C
278 lines
6.3 KiB
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:c:D:d:h:Ll:N:n:w:x:y:z:"
|
|
|
|
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 = ncriticals();
|
|
fprintf(stderr, "Critical points:\n");
|
|
for (i = 0; i < n; i++) {
|
|
z = get_critical(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 = ncriticals();
|
|
for (i = 0; i < n; i++) {
|
|
d = cabsl(z - get_critical(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] [-x julia_x] [-y julia_y]\n"
|
|
" [-n niters] [-h height] [-w width] [-z zoom]\n"
|
|
" [-a center_x] [-b center_y] [-l light_intensity]\n"
|
|
" [-d density] [-c contour] [-D displacement]\n"
|
|
" [-N newton_iters]\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 = 720,
|
|
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':
|
|
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 'h':
|
|
height = parse_integer(optarg, 1, 65536);
|
|
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, 65536);
|
|
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 * 3.14159265358979323846L;
|
|
angle_z = (0.2L + 0.6L * drand48()) * 3.14159265358979323846L;
|
|
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 (behaves_critical(df, ddf, z, &crit, newton_iters,
|
|
bailout, epsilon))
|
|
add_critical_point(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);
|
|
}
|