#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <ncurses.h>

#define REALLOC_STRIDE 0x40

enum { X, Y, Z, Ndim };

typedef struct {
    double (*edge)[2][Ndim];
    int numedges, resedges;
    double xang, yang;
} model_t;

void loop(model_t *);
model_t *model_new(FILE *);
void model_delete(model_t *);
void model_addedge(model_t *, double, double, double, double, double, double);
void model_rotate(model_t *, double, double);
void model_print(model_t *, double, double, double);

int main(int argc, char **argv)
{
    FILE *file = fopen(argv[1], "r");
    model_t *model = model_new(file);
    fclose(file);

    initscr();
    cbreak();
    noecho();

    loop(model);
    model_delete(model);

    echo();
    nocbreak();
    endwin();
}

void loop(model_t *model)
{
    timeout(0);

    double dist = 100, zoom = 100;
    while (1) {
        napms(100);
        int key;
        while ((key = getch()) != ERR) {
            switch (key) {
                case 'q': dist += 5; break;
                case 'a': dist -= 5; break;
                case 'w': zoom *= 1.125; break;
                case 's': zoom /= 1.125; break;
                case 'e': model_rotate(model, 0.0,  0.1); break;
                case 'd': model_rotate(model, 0.0, -0.1); break;
                case 'x': goto nahui;
            }
        }
        model_rotate(model, 0.1, 0.0);
        model_print(model, dist, zoom, 0.5);
    }
    nahui:

    timeout(-1);
}

model_t *model_new(FILE *file)
{
    model_t *this = malloc(sizeof(model_t));

    this->edge = malloc(REALLOC_STRIDE * sizeof(*this->edge));
    this->numedges = 0;
    this->resedges = REALLOC_STRIDE;
    this->yang = this->xang = 0;

    double x0, y0, z0, x1, y1, z1;
    while (fscanf(file, "%lf%lf%lf%lf%lf%lf", &x0, &y0, &z0, &x1, &y1, &z1) == 6)
    {
        model_addedge(this, x0, y0, z0, x1, y1, z1);
    }

    return this;
}

void model_delete(model_t *this)
{
    free(this->edge);
    free(this);
}

void model_addedge(model_t *this,
                   double x0, double y0, double z0,
                   double x1, double y1, double z1)
{
    if (this->numedges == this->resedges) {
        this->edge = realloc(this->edge, (this->resedges += REALLOC_STRIDE)
                                         * sizeof(*this->edge));
    }
    this->edge[this->numedges][0][X] = x0;
    this->edge[this->numedges][0][Y] = y0;
    this->edge[this->numedges][0][Z] = z0;
    this->edge[this->numedges][1][X] = x1;
    this->edge[this->numedges][1][Y] = y1;
    this->edge[this->numedges][1][Z] = z1;
    ++this->numedges;
}

void model_rotate(model_t *this, double yang, double xang)
{
    this->yang += yang;
    this->xang += xang;
}

void model_print(model_t *this, double dist, double zoom, double aspect)
{
    erase();
    move(0, 0);
    for (int i = 0; i < this->numedges; ++i) {
        double xcos = cos(this->xang),
               xsin = sin(this->xang),
               ycos = cos(this->yang),
               ysin = sin(this->yang),
               x0t = ycos * this->edge[i][0][X] + ysin * this->edge[i][0][Z],
               y0t = this->edge[i][0][Y],
               z0t = ycos * this->edge[i][0][Z] - ysin * this->edge[i][0][X],
               x0 = x0t,
               y0 = xcos * y0t + xsin * z0t,
               z0 = xcos * z0t - xsin * y0t + dist,
               x1t = ycos * this->edge[i][1][X] + ysin * this->edge[i][1][Z],
               y1t = this->edge[i][1][Y],
               z1t = ycos * this->edge[i][1][Z] - ysin * this->edge[i][1][X],
               x1 = x1t,
               y1 = xcos * y1t + xsin * z1t,
               z1 = xcos * z1t - xsin * y1t + dist;

        if (z0 < 5 || z1 < 5) continue;

        int cx = getmaxx(stdscr) / 2,
            cy = getmaxy(stdscr) / 2,
            x0_ = cx + zoom * x0 / z0,
            y0_ = cy + zoom * y0 / z0 * aspect,
            x1_ = cx + zoom * x1 / z1,
            y1_ = cy + zoom * y1 / z1 * aspect,
            len = (abs(x0_ - x1_) > abs(y0_ - y1_)) ? abs(x0_ - x1_) : abs(y0_ - y1_);
        for (int p = 0; p <= len; ++p)
            mvaddch(y0_ + (y1_ - y0_) * p / len, x0_ + (x1_ - x0_) * p / len, '$');
    }
    refresh();
}