Mandelbrot set rendered from a first-person viewpoint

hipparchus.c

// I don't know what I'm looking at here.

#define _USE_MATH_DEFINES
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>

#define OUTPUT			"hipparchus.pgm"
#define WIDTH			1024
#define HEIGHT			256
#define NPARTICLES		10000000
#define NITER			100
#define EXPOSURE		0.0001

uint64_t seed = 1;

double random( double a, double b ) // this is xorshift64, but any old random number source will work
{
	seed ^= seed << 13;
	seed ^= seed >> 7;
	seed ^= seed << 17;
	return a + seed * (b-a) / (double)(UINT64_MAX);
}

int main( int argc, char *argv[] )
{
	int* histogram = calloc( WIDTH*HEIGHT, sizeof(int) );

	for ( int i = 0; i<NPARTICLES; i++ )
	{
	again:;
		// pick a random point on the 2D complex plane
		double cx = random( -2, +2 );
		double cy = random( -2, +2 );

		// iterate it until we get bored
		double zx = cx, zy = cy;
		for ( int iter = 0; iter<NITER; iter++ )
		{
			// z' = z*z + c
			double zx2 = zx*zx - zy*zy + cx;
			double zy2 = zx*zy + zy*zx + cy;
			zx = zx2;
			zy = zy2;

			// cartesian -> polar unwrap
			double x = atan2( zy, zx );
			double y = sqrt( zx*zx + zy*zy );
			if ( y > 2 )
				goto again; // reject non-stable orbits

			// find pixel coordinates
			int ix = (int)floor( x / ( 2.0*M_PI )  * WIDTH );
			int iy = (int)floor( ( 1.0 - y / 2.0 ) * HEIGHT );
			if ( ix < 0 )
				ix += WIDTH;

			// plot every point we pass through on our stable orbit
			histogram[ iy*WIDTH+ix ]++;
		}
	}

	// save image to PGM format
	FILE* fp = fopen( OUTPUT, "wb" );
	if ( !fp ) {
		perror( OUTPUT );
		return 1;
	}
	fprintf( fp, "P5 %i %i 255\n", WIDTH, HEIGHT );
	for ( int i = 0; i<WIDTH*HEIGHT; i++ )
	{
		// apply simple camera exposure curve
		int lum = (int)( 255 * ( 1 - exp( histogram[ i ] * -EXPOSURE ) ) );
		fputc( lum, fp );
	}
	fclose( fp );

	return 0;
}