mand.d

import std.stdio : writeln;
import std.complex : complex, Complex, abs, arg;
import std.algorithm.comparison : clamp;
import conv = std.conv;
import math = std.math;
import IF = imageformats;

alias Image = IF.IFImage;

immutable int size = 64; // 画像サイズ
immutable double cx = 0.001643721971153; // 描画域の中心X座標
immutable double cy = 0.822467633298876; // 描画域の中心Y座標
immutable double scale = 4.0 / size / 1e11 / 4; // 1 ピクセルは座標平面においていくつか？

immutable int max_iter = 2048;

Complex!double exp (Complex!double z) {
	return math.exp(z.re) * complex!double(math.cos(z.im), math.sin(z.im));
}

Complex!double cos (Complex!double z) {
	return (exp(complex!double(-z.im, z.re)) + exp(complex!double(z.im, -z.re))) / 2.0;
}

Complex!double sin (Complex!double z) {
	return (exp(complex!double(-z.im, z.re)) - exp(complex!double(z.im, -z.re))) / complex!double(2.0, 1.0);
}

int mand (Complex!double z, Complex!double c, int n) {
	if (n >= max_iter) {
		return 0;
	}
	// z = complex!double(math.abs(z.re), math.abs(z.im));
	auto z1 = z * z + c;
	if (z1.re * z1.re + z1.im * z1.im > 4.0) {
		return n;
	}
	return mand(z1, c, n + 1);
}

void main () {
	ubyte[] pixels;
	pixels.length = size * size;
	for (int y = 0; y < size; y++) {
		double ry = (y - size / 2) * scale + cy;
		for (int x = 0; x < size; x++) {
			double rx = (x - size / 2) * scale + cx;
			// auto z = complex!double(0, 0);
			auto c = complex!double(rx, ry);
			int n = mand(c, c, 0);
			pixels[y * size + x] = conv.to!ubyte(
				clamp(n * 0xFF / max_iter, 0, 0xFF)
			);
		}
	}
	IF.write_image("out.png", size, size, pixels, IF.ColFmt.Y);
}