#include <unistd.h>
#include <stdint.h>
#include <stdlib.h>
#include <math.h>

#define P_SZ 4194304
#define M_SZ 262144
#define TOT_SPL (P_SZ / 16)
#define TPF 100

struct __attribute__((packed)) Splat {
    int16_t x, y, z;
    uint16_t t0;
    uint8_t dt;
    uint8_t op;
    uint8_t r, g, b, pad;
};

static struct Splat p_pool[TOT_SPL];
static int8_t m_pool[TOT_SPL];
static uint8_t tar_hdr[512];
static char path_buf[128];

static uint32_t seed = 9876;
static uint32_t djb_rand(void) {
    seed = seed * 1103515245 + 12345;
    return (seed / 65536) % 32768;
}

static void put_octal(uint8_t *buf, uint32_t val, int size) {
    int i;
    buf[size - 1] = ' ';
    for (i = size - 2; i >= 0; i--) {
        buf[i] = '0' + (val % 8);
        val /= 8;
    }
}

static void write_tar_header(const char *name, uint32_t size_dec) {
    int i; uint32_t chk = 0;
    for (i = 0; i < 512; i++) tar_hdr[i] = 0;
    for (i = 0; name[i] != '\0' && i < 100; i++) tar_hdr[i] = name[i];
    put_octal(tar_hdr + 100, 0000644, 8);
    put_octal(tar_hdr + 108, 0000000, 8);
    put_octal(tar_hdr + 116, 0000000, 8);
    if (size_dec == 4194304) put_octal(tar_hdr + 124, 020000000, 12);
    else if (size_dec == 262144) put_octal(tar_hdr + 124, 01000000, 12);
    else put_octal(tar_hdr + 124, size_dec, 12);
    put_octal(tar_hdr + 136, 00000000000, 12);
    tar_hdr[156] = '0';
    for (i = 0; i < 8; i++) tar_hdr[148 + i] = ' ';
    for (i = 0; i < 512; i++) chk += tar_hdr[i];
    put_octal(tar_hdr + 148, chk, 7);
    tar_hdr[155] = '\0';
    write(1, tar_hdr, 512);
}

static void build_filename(char *buf, int32_t sz, int32_t zeit, char sfx) {
    int32_t ux = 5000, uy = 5000, uz = sz + 5000;
    char *p = buf;
    *p++ = 'w'; *p++ = 'o'; *p++ = 'r'; *p++ = 'l'; *p++ = 'd'; *p++ = '/';
    *p++ = '0' + (ux / 1000); *p++ = '0' + ((ux / 100) % 10); *p++ = '0' + ((ux / 10) % 10); *p++ = '0' + (ux % 10); *p++ = '/';
    *p++ = '0' + (uy / 1000); *p++ = '0' + ((uy / 100) % 10); *p++ = '0' + ((uy / 10) % 10); *p++ = '0' + (uy % 10); *p++ = '/';
    *p++ = '0' + (uz / 1000); *p++ = '0' + ((uz / 100) % 10); *p++ = '0' + ((uz / 10) % 10); *p++ = '0' + (uz % 10); *p++ = '/';
    *p++ = '0' + (zeit / 10000); *p++ = '0' + ((zeit / 1000) % 10); *p++ = '0' + ((zeit / 100) % 10); *p++ = '0' + ((zeit / 10) % 10); *p++ = '0' + (zeit % 10);
    *p++ = '.'; *p++ = sfx; *p = '\0';
}

int main(void) {
    int32_t sz, tp; uint32_t i;
    for (sz = 0; sz <= 5; sz++) {
        for (tp = 0; tp <= 4; tp++) {
            int32_t zeit = tp * 100;
            double phase = (double)zeit * 0.02;

            for (i = 0; i < TOT_SPL; i++) {
                int32_t gx = djb_rand() % 500;
                int32_t gy = djb_rand() % 500;
                int32_t gz = djb_rand() % 500;

                p_pool[i].x = gx; p_pool[i].y = gy; p_pool[i].z = gz;
                p_pool[i].t0 = 0; p_pool[i].dt = 100; p_pool[i].op = 0; m_pool[i] = 0;

                double dx = (double)(gx - 250);
                double dy = (double)(gy - 250);
                double dz = (double)(gz - 250);
                double r2d = sqrt(dx*dx + dz*dz);
                double r3d = sqrt(dx*dx + dy*dy + dz*dz);

                if (r3d > 10.0 && r3d < 240.0) {
                    double angle = atan2(dz, dx);
                    double target_r = 90.0 + 40.0 * sin(dy * 0.04 + phase) + 25.0 * cos(angle * 5.0 + dy * 0.02);
                    double dist_to_surface = fabs(r2d - target_r);

                    if (dist_to_surface < 12.0) {
                        p_pool[i].op = 255;
                        p_pool[i].r = (uint8_t)(128.0 + 127.0 * sin(angle + phase));
                        p_pool[i].g = (uint8_t)(128.0 + 127.0 * cos(dy * 0.03));
                        p_pool[i].b = (uint8_t)(200.0 + 55.0 * sin(phase * 0.5));

                        double nx = dx / (r2d + 1.0);
                        double nz = dz / (r2d + 1.0);
                        double specular = fabs(nx * cos(angle) + nz * sin(angle));
                        m_pool[i] = (int8_t)(15.0 + 110.0 * specular);
                    }
                }

                if (p_pool[i].op == 0 && r3d < 40.0 && (djb_rand() % 1000) < 8) {
                    p_pool[i].op = 200;
                    p_pool[i].r = 0; p_pool[i].g = 255; p_pool[i].b = 255;
                    p_pool[i].t0 = djb_rand() % 50; p_pool[i].dt = 30 + (djb_rand() % 40);
                    m_pool[i] = 127;
                }
            }

            build_filename(path_buf, sz, zeit, 'p');
            write_tar_header(path_buf, P_SZ); write(1, p_pool, P_SZ);
            build_filename(path_buf, sz, zeit, 'm');
            write_tar_header(path_buf, M_SZ); write(1, m_pool, M_SZ);
        }
    }

    write_tar_header("finished", 512);
    for (i = 0; i < 512; i++) ((uint8_t*)p_pool)[i] = 0;
    write(1, p_pool, 512);
    for (i = 0; i < 512; i++) tar_hdr[i] = 0;
    write(1, tar_hdr, 512); write(1, tar_hdr, 512);
    return 0;
}