Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
/* $ cc -std=c99 -Ofast -fopenmp -o julia julia.c -lm
* $ ./julia > output.ppm
* $ ppmtoy4m -F 60:1 < output.ppm > output.y4m
* $ x264 --qp 0 -o output.mp4 output.y4m
*
* Due to parallel compute, stdout must be seekable. Therefore it
* cannot be piped directly into ppmtoy4m.
*/
#define IMAGE_SIZE 600
#define IMAX 512
#define FRAMES 540
#define RADIUS 1.0
#define HEADER_LEN 21
#ifdef _WIN32
#define WIN32_LEAN_AND_MEAN
#include <windows.h>
static void
write_frame(const void *buf, int n)
{
size_t len = HEADER_LEN + 3UL * IMAGE_SIZE * IMAGE_SIZE;
DWORD out;
OVERLAPPED o = {.Offset = len * n};
WriteFile(GetStdHandle(STD_OUTPUT_HANDLE), buf, len, &out, &o);
}
#else /* POSIX */
#define _POSIX_C_SOURCE 200809L
#include <unistd.h>
static void
write_frame(const void *buf, int n)
{
size_t len = HEADER_LEN + 3UL * IMAGE_SIZE * IMAGE_SIZE;
pwrite(1, buf, len, len * n);
}
#endif
#include <math.h>
#include <stdio.h>
#include <complex.h>
int
main(void)
{
#pragma omp parallel for schedule(dynamic, 1)
for (int f = 0; f < FRAMES; f++) {
double theta = f * 3.141592653589793 * 2 / FRAMES;
unsigned char buf[HEADER_LEN + 3UL * IMAGE_SIZE * IMAGE_SIZE];
sprintf((char *)buf, "P6 % 6d % 6d 255\n", IMAGE_SIZE, IMAGE_SIZE);
unsigned char *p = buf + HEADER_LEN;
for (int y = 0; y < IMAGE_SIZE; y++) {
for (int x = 0; x < IMAGE_SIZE; x++) {
complex double z =
((x / (IMAGE_SIZE - 1.0) - 0.5) * RADIUS * 2) +
((y / (IMAGE_SIZE - 1.0) - 0.5) * RADIUS * 2) * I;
complex double c =
(-0.80 + 0.1 * sin(theta)) +
(+0.00 + 0.1 * cos(theta)) * I;
int i = 0;
for (; i < IMAX && cabs(z) <= 2.0; i++)
z = z * z + c;
int v = powf(i / (float)IMAX, 1 / 2.2f) * 255;
*p++ = v;
*p++ = v;
*p++ = v;
}
}
write_frame(buf, f);
}
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment