Skip to content

Instantly share code, notes, and snippets.

@jabbalaci
Created May 9, 2020 16:34
Show Gist options
  • Save jabbalaci/5a6a5f75f6ef207fda9975f2bbd8c50d to your computer and use it in GitHub Desktop.
Save jabbalaci/5a6a5f75f6ef207fda9975f2bbd8c50d to your computer and use it in GitHub Desktop.
Lorenz system animation
/* Lorenz system animation
* $ cc -Ofast lorenz.c -lm
* $ ./a.out | mpv --no-correct-pts --fps=60 --fs -
* $ ./a.out | x264 --fps 60 --frames 3600 -o lorenz.mp4 /dev/stdin
* This is free and unencumbered software released into the public domain.
* Ref: https://nullprogram.com/blog/2017/11/03/
* Ref: https://nullprogram.com/blog/2018/02/15/
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#define N 256
#define W 1920
#define H 1080
#define SCALE 20.0f
#define MAX 128
#define PI 3.1415927f
static unsigned long long u32s;
static unsigned long u32(void)
{
return (u32s = u32s*0xfc5434fdb4a9e74d + 1) >> 32 & 0xffffffff;
}
static float rnd(float min, float max)
{
return ldexpf(u32(), -32)*(max - min) + min;
}
struct v2 { int x, y; };
struct v3 { float x, y, z; };
static struct v2
ortho(struct v3 p, float scale, float az, float el)
{
float x = p.x*cosf(az) - p.y*sinf(az);
float y = p.x*sinf(az) + p.y*cosf(az);
return (struct v2){
roundf(scale*x + W/2),
roundf(scale*(y*sinf(el) + p.z*cosf(el)) + H/2)
};
}
static unsigned char buf[H][W][3];
static void
clear(void)
{
memset(buf, 0, sizeof(buf));
}
static void
point(struct v2 p, long color)
{
if (p.x >= 0 && p.x < W && p.y >= 0 && p.y < H) {
buf[p.y][p.x][0] = color >> 16;
buf[p.y][p.x][1] = color >> 8;
buf[p.y][p.x][2] = color >> 0;
}
}
static void
line(struct v2 a, struct v2 b, long color)
{
int dx = abs(b.x - a.x);
int dy = abs(b.y - a.y);
int sx = b.x < a.x ? -1 : 1;
int sy = b.y < a.y ? -1 : 1;
if (dx > dy) {
int d = 2 * dy - dx;
struct v2 p = {0, a.y};
for (p.x = a.x; p.x != b.x; p.x += sx) {
point(p, color);
if (d > 0) {
p.y += sy;
d -= 2 * dx;
}
d += 2 * dy;
}
} else {
int d = 2 * dx - dy;
struct v2 p = {a.x, 0};
for (p.y = a.y; p.y != b.y; p.y += sy) {
point(p, color);
if (d > 0) {
p.x += sx;
d -= 2 * dy;
}
d += 2 * dx;
}
}
point(b, color);
}
static void
frame(void)
{
#define xstr(s) str(s)
#define str(s) #s
const char header[] = "P6\n" xstr(W) " " xstr(H) "\n255\n";
if (!fwrite(header, sizeof(header) - 1, 1, stdout) ||
!fwrite(buf, sizeof(buf), 1, stdout)) {
exit(1);
}
}
static struct v3
lorenz(struct v3 s)
{
#define SIGMA 10.0f
#define BETA 8.0f/3.0f
#define RHO 28.0f
#define DT 0.004f
#define DX(x, y, z) (SIGMA*(y - x))
#define DY(x, y, z) (x*(RHO - z) - y)
#define DZ(x, y, z) (x*y - BETA*z)
float x = s.x;
float y = s.y;
float z = s.z;
float k1dx = DX(x, y, z);
float k1dy = DY(x, y, z);
float k1dz = DZ(x, y, z);
float k2x = x + k1dx * DT / 2;
float k2y = y + k1dy * DT / 2;
float k2z = z + k1dz * DT / 2;
float k2dx = DX(k2x, k2y, k2z);
float k2dy = DY(k2x, k2y, k2z);
float k2dz = DZ(k2x, k2y, k2z);
float k3x = x + k2dx * DT / 2;
float k3y = y + k2dy * DT / 2;
float k3z = z + k2dz * DT / 2;
float k3dx = DX(k3x, k3y, k3z);
float k3dy = DY(k3x, k3y, k3z);
float k3dz = DZ(k3x, k3y, k3z);
float k4x = x + k3dx * DT;
float k4y = y + k3dy * DT;
float k4z = z + k3dz * DT;
float k4dx = DX(k4x, k4y, k4z);
float k4dy = DY(k4x, k4y, k4z);
float k4dz = DZ(k4x, k4y, k4z);
s.x = x + (k1dx + 2*k2dx + 2*k3dx + k4dx) * DT / 6;
s.y = y + (k1dy + 2*k2dy + 2*k3dy + k4dy) * DT / 6;
s.z = z + (k1dz + 2*k2dz + 2*k3dz + k4dz) * DT / 6;
return s;
};
static struct v3
translate(struct v3 v)
{
v.z -= RHO;
return v;
}
int
main(void)
{
#ifdef _WIN32
int _setmode(int, int);
_setmode(1, 0x8000); /* Set stdin/stdout to binary mode. */
#endif
struct {
long color;
int head, tail;
struct v3 pos[MAX];
} points[N];
u32s = time(0);
for (int i = 0; i < N; i++) {
points[i].color = u32()>>8 | 0x404040;
points[i].head = points[i].tail = 0;
points[i].pos[0].x = rnd(-0.01f, +0.01f) + 25;
points[i].pos[0].y = rnd(-0.01f, +0.01f) + 25;
points[i].pos[0].z = rnd(-0.01f, +0.01f);
}
float az = rnd(-PI, PI);
float elv = rnd(-PI, PI);
for (;;) {
clear();
float el = sinf(elv)*0.5f;
for (int i = 0; i < N; i++) {
int head = (points[i].head + 1)%MAX;
int tail = points[i].tail;
if (head == tail) {
tail = (tail + 1)%MAX;
}
points[i].pos[head] = lorenz(points[i].pos[points[i].head]);
points[i].head = head;
points[i].tail = tail;
for (int j = head; j != tail; j = (j + MAX - 1)%MAX) {
struct v3 h = translate(points[i].pos[j]);
struct v3 t = translate(points[i].pos[(j+MAX-1)%MAX]);
struct v2 a = ortho(h, SCALE, az, el);
struct v2 b = ortho(t, SCALE, az, el);
line(a, b, points[i].color);
}
struct v3 h = translate(points[i].pos[head]);
struct v2 r = ortho(h, SCALE, az, el);
point((struct v2){r.x+1, r.y+0}, points[i].color);
point((struct v2){r.x-1, r.y+0}, points[i].color);
point((struct v2){r.x+0, r.y+1}, points[i].color);
point((struct v2){r.x+0, r.y-1}, points[i].color);
}
frame();
az += 0.0027f;
elv += 0.01f;
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment