attractor_visualizer.c

#include <math.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>

#define W 1920
#define H 1080

const double start_x = 0.1;
const double start_y = 0;
const double start_z = 0;

const double time_step = 0.0005;
const int steps_per_frame = 10000;

const double brightness_smoothing = 0.99;
const double decay = 0.99;

const int zoom = H;

const int header_size = 128;
const int buff_size = (W * H + header_size);

typedef struct {
    double dx, dy, dz;
    double x, y, z;

    double t, dt;
} point_t;

void action(point_t *p) {
    double x = p->x;
    double y = p->y;
    double z = p->z;

    double t = p->t;
    double dt = p->dt;

    double b = 0.185;

    p->dx = sin(y) - b * x;
    p->dy = sin(z) - b * y;
    p->dz = sin(x) - b * z;

    p->x += p->dx * dt;
    p->y += p->dy * dt;
    p->z += p->dz * dt;

    p->t += dt;
}

double project(point_t p, double point[2], double direction[3]) {
    double x = p.x;
    double y = p.y;

    double dx = p.dx;
    double dy = p.dy;

    point[0] = x / 5;
    point[1] = y / 5;

    double norm = hypot(dx, dy);
    direction[0] = dx / norm;
    direction[1] = dy / norm;

    return norm;
}

double saturate(double x) { return x / (1.0 + fabs(x)); }

double squared_magnitude(double x, double y) { return x * x + y * y; }

double velocity_blending(double old, double new, double velocity) {
    velocity += 0.5;
    return fmax(old / velocity, 1.0 / (1.0 + 0.5 * new)) * velocity;
}

int main(void) {
    char *out_buff = malloc(buff_size);
    int img_start = snprintf(out_buff, header_size, "P5\n%d %d\n255\n", W, H);
    char(*img_data)[W] = (char(*)[W])(out_buff + img_start);
    int out_length = img_start + W * H;

    static double plane[W][H] = {0};

    point_t state = {0, 0, 0, start_x, start_y, start_z, 0, time_step};

    double avg_val_smooth = 0.1;

    while (true) {
        for (int i = 0; i < steps_per_frame; i++) {
            double draw_point[2];
            double direction[2];
            double velocity;

            action(&state);
            velocity = project(state, draw_point, direction);

            double sx = (draw_point[0] * zoom + W) * 0.5;
            double sy = (draw_point[1] * zoom + H) * 0.5;

            for (int j = -8; j <= 8; j++) {
                for (int k = -8; k <= 8; k++) {
                    int ix = sx + (direction[1]) * j;
                    int iy = sy + (direction[0]) * k;

                    ix += W;
                    iy += H;
                    ix %= W;
                    iy %= H;

                    plane[ix][iy] = velocity_blending(plane[ix][iy], squared_magnitude(ix - sx, iy - sy), velocity);
                }
            }
        }

        double avg_val = 0;
        int count = 1;

        for (int ix = 0; ix < W; ix++) {
            for (int iy = 0; iy < H; iy++) {
                if (plane[ix][iy] > avg_val_smooth / 5)
                    plane[ix][iy] *= decay;

                if (plane[ix][iy] < 0.00001)
                    continue;

                avg_val += plane[ix][iy];
                count += 1;
            }
        }

        avg_val = fmax(avg_val / count, 0.001);
        avg_val_smooth = avg_val * brightness_smoothing + avg_val * (1 - brightness_smoothing);

        for (int i = 0; i < W; i++) {
            for (int j = 0; j < H; j++) {
                img_data[j][i] = (int)(0.5 + saturate(plane[i][j] / avg_val_smooth) * 254.0);
            }
        }

        fwrite(out_buff, sizeof(char), out_length, stdout);
    }
}