#include <u.h>
#include <libc.h>
#include <draw.h>
typedef struct Vec Vec; struct Vec{float x, y, z;};
Vec add(Vec a, Vec b){return (Vec){a.x+b.x, a.y+b.y, a.z+b.z};}
Vec sub(Vec a, Vec b){return (Vec){a.x-b.x, a.y-b.y, a.z-b.z};}
Vec mul(Vec p, float s){return (Vec){p.x*s, p.y*s, p.z*s};}
float dot(Vec a, Vec b){return a.x*b.x + a.y*b.y + a.z*b.z;}
float mag(Vec v){return sqrt(dot(v, v));}
Vec norm(Vec v){
	float len;
	len = mag(v);
	return len == 0 ? (Vec){0,0,0} : (Vec){v.x/len, v.y/len, v.z/len};
}
Vec cross(Vec a, Vec b){
	return (Vec){
		a.y*b.z - a.z*b.y,
		a.z*b.x - a.x*b.z,
		a.x*b.y - a.y*b.x,
	};
}
struct{Vec pos, right, up, fwd;}camera;
void
mvcam(float x0, float y0, float z0, float x1, float y1, float z1){
	camera.pos = (Vec){x0, y0, z0};
	camera.fwd = norm((Vec){x1 - x0, y1 - y0, z1 - z0});
	camera.right = cross((Vec){0.0, 1.0, 0.0}, camera.fwd);
	camera.up = cross(camera.fwd, camera.right);
}
float
ring(Vec p, float r1, float r2, float min){
	float d;
	d = sqrt(p.y * p.y + p.x * p.x) - r1;
	d = sqrt(d * d + p.z * p.z) - r2;
	return min < d ? min : d;
}
float
smoothmin(float a, float b, float k){
	float h;
	if(a > b){
		h = a;
		a = b;
		b =h;
	}
	h = (k - b + a) / k;
	if(h < 0.f)
		h = 0.f;
	return a - h * h * k * 0.25f;
}
float
tail(Vec p, float min){
	float d;
	p.y -= 0.9f;
	p.x += 0.35f;
	d = p.x;
	p.x = 0.56529953116 * p.x - 0.824885713338 * p.y;
	p.y = 0.824885713338 * d + 0.56529953116 * p.y;
	p.x = fabs(p.x);
	d = .5f * p.x > 0.866025f * p.y ? p.x * 0.866025f + p.y * 0.5f : sqrt(p.x * p.x + p.y * p.y);
	d = sqrt(dot(p, p) + 1.4f * 1.4f - 2.8f * d) - 0.5;
	return smoothmin(d, min, 0.5f);
}
float
dist(Vec p){
	float d;
	d = ring(p, 3.f, 0.5f, 99.f);
	p.y += 0.7;
	d = ring(p, 0.95f, 0.5f, d);
	return tail(p, d);
}
void
marchrow(uchar *out, int y){
	Vec p, n, lux;
	float d, tot;
	int x, i;
	for(x = -32; x < 32; x++){
		p = add(camera.fwd, mul(camera.right, x / 32.0));
		p = add(p, mul(camera.up, y / 32.0));
		n = norm(p);
		p = camera.pos;
		for(i = tot = 0.0; i < 50; i++){
			tot += d = dist(p);
			if(tot >= 80.0 || d < 0.05f)
				break;
			p = add(p, mul(n, d));
		}
		lux = norm(sub(camera.pos, p));
		n.x = (dist((Vec){p.x + 0.05f, p.y, p.z}) - d) / 0.05f;
		n.y = (dist((Vec){p.x, p.y + 0.05f, p.z}) - d) / 0.05f;
		n.z = (dist((Vec){p.x, p.y, p.z + 0.05f}) - d) / 0.05f;
		d = dot(lux, n) * 255;
		*out++ = d < 0.0 ? 0 : d > 255.0 ? 255 : d;
	}
}
Image *scanline;
Rectangle scanlinerect;
uchar hues[4][13];
hue;
void
quadruple(uchar *p, uchar *a, uchar *b, uchar *c, uchar *d){
	uchar out[512 * 2], *o;
	int i;
	o = out;
	for(i = 0; i < 64; i++){
		o[0] = o[1] = o[2] = o[3] = c[i] == a[i] && c[i] != d[i] && a[i] != b[i] ? a[i] : p[i];
		o[4] = o[5] = o[6] = o[7] = a[i] == b[i] && a[i] != c[i] && b[i] != d[i] ? b[i] : p[i];
		o[512] = o[513] = o[514] = o[515] = d[i] == c[i] && d[i] != b[i] && c[i] != a[i] ? c[i] : p[i];
		o[516] = o[517] = o[518] = o[519] = b[i] == d[i] && b[i] != a[i] && d[i] != c[i] ? d[i] : p[i];
		o[0] = o[1] = o[2] = o[3] = c[i] + a[i] + p[i] + p[i] >> 2;
		o[4] = o[5] = o[6] = o[7] = a[i] + b[i] + p[i] + p[i] >> 2;
		o[512] = o[513] = o[514] = o[515] = c[i] + d[i] + p[i] + p[i] >> 2;
		o[516] = o[517] = o[518] = o[519] = d[i] + b[i] + p[i] + p[i] >> 2;
		o += 8;
	}
	for(i = 0; i < 1024; i++)
		out[i] = hues[hue][out[i] * 13 >> 8];
	scanlinerect.max.y = scanlinerect.min.y + 4;
	loadimage(scanline, scanline->r, out, 512);
	draw(screen, scanlinerect, scanline, nil, ZP);
	scanlinerect.min.y = scanlinerect.max.y;
	scanlinerect.max.y = scanlinerect.min.y + 4;
	loadimage(scanline, scanline->r, o, 512);
	draw(screen, scanlinerect, scanline, nil, ZP);
	scanlinerect.min.y = scanlinerect.max.y;
}
void
nuscan(void){
	uchar *a, *d, *p, *t;
	uchar buf[(64 + 2) * 3];
	Rectangle r;
	int i, j;
	scanlinerect.min.x = screen->r.min.x + screen->r.max.x - 512 >> 1;
	scanlinerect.max.x = scanlinerect.min.x + 512;
	scanlinerect.min.y = screen->r.min.y + screen->r.max.y - 512 >> 1;
	memset(buf, 0, sizeof buf);
	a = buf + 1;
	p = a + 64 + 2;
	d = p + 64 + 2;
	marchrow(p, -32);
	for(i = 1; i < 64; i++){
		marchrow(d, i - 32);
		quadruple(p, a, p + 1, p - 1, d);
		t = a;
		a = p;
		p = d;
		d = t;
	}
	memset(d, 0, 64);
	quadruple(p, a, p + 1, p - 1, d);
	flushimage(display, 1);
}
void
main(int, char**){
	int t;
	if(initdraw(nil, nil, "⑨") < 0)
		sysfatal("initdraw: %r");
	scanline = allocimage(display, Rect(0, 0, 512, 1), CMAP8, 1, DNofill);
	if(scanline == nil)
		sysfatal("allocimage: %r");
	draw(screen, screen->r, display->black, nil, ZP);
	for(t = 1; t < 13; t++){
		hues[1][t] = 47 + t * 17;
		hues[2][t] = 50 + t * 17;
		hues[3][t] = -16 + t * 17;
		hues[0][t] = -13 + t * 17;
	}
	mvcam(0.0, 0.0, -5.5, 0, 0, 0);
	for(t = 0; t < 1200; t++){
		mvcam(-5.5 * sin(t * 2.0 * PI / 120.0), 0.0, -5.5 * cos(t * 2.0 * PI / 120.0), 0, 0, 0);
		if(!(t & 31))
			hue = hue + 1 & 3;
		nuscan();
		sleep(50);
	}
	exits(nil);
}