#include <math.h>
#include "svg.h"

#define CENTER_X  148.5
#define CENTER_Y  105.0
#define SCALE  4

#define COLOR_DELTA  17
#define COLOR_MAX  0xFF
#define COLOR_MIN  0x00

/*
  #define P  10.0
  #define R  28.0
  #define B  2.6666667

  #define dfXdT(t, x, y, z)  ( \
  (P)*((y) - (x)) \
  )

  #define dfYdT(t, x, y, z)  ( \
  (R)*(x) - (x)*(z) - (y) \
  )

  #define dfZdT(t, x, y, z)  ( \
  (x)*(y) - (B)*(z) \
  )
//*/

#define P  10.0
#define R  28.0
#define B  2.6666667

#define dfXdT(t, x, y, z)  (                            \
                            (P)*((y) - (x)) + (t)       \
                            )

#define dfYdT(t, x, y, z)  (                            \
                            (R)*(x) - (x)*(z) - (y)     \
                            )

#define dfZdT(t, x, y, z)  (                            \
                            (x)*(y) - (B)*(z) - (t)     \
                            )

double mySin;
double myCos;

// 座標変換と出力
void output(
            double xold, double yold, double zold,
            double xnew, double ynew, double znew
            ) {
    static int cg = COLOR_MAX;
    static int cb = COLOR_MIN;
    static int mode = 0;

    double xo = xold*myCos - yold*mySin;
    double yo = xold*mySin + yold*myCos;
    double xn = xnew*myCos - ynew*mySin;
    double yn = xnew*mySin + ynew*myCos;

    xo *= SCALE;
    yo *= SCALE;
    xn *= SCALE;
    yn *= SCALE;

    xo += CENTER_X;
    yo += CENTER_Y;
    xn += CENTER_X;
    yn += CENTER_Y;

    switch (mode) {
    case 0: cb += COLOR_DELTA; if(cb == COLOR_MAX) mode++; break;
    case 1: cg -= COLOR_DELTA; if(cg == COLOR_MIN) mode++; break;
    case 2: cb -= COLOR_DELTA; if(cb == COLOR_MIN) mode++; break;
    case 3: cg += COLOR_DELTA; if(cg == COLOR_MAX) mode=0; break;
    }

    stroke(rgb255(0, cg,  cb));
    line(xo, yo, xn, yn);
}


// ルンゲクッタ法
void rungekutta(
                double t0, double tx, int cnt,
                double x0, double y0, double z0
                ) {
    int i;
    double dt = (tx - t0) / cnt;

    double tnew, told = t0;
    double xnew, xold = x0;
    double ynew, yold = y0;
    double znew, zold = z0;

    for (i = 1; i <= cnt; i++) {
        double kx1 = dt*dfXdT(told, xold, yold, zold);
        double ky1 = dt*dfYdT(told, xold, yold, zold);
        double kz1 = dt*dfZdT(told, xold, yold, zold);

        double kx2 = dt*dfXdT(told + dt/2, xold + kx1/2, yold + ky1/2, zold + kz1/2);
        double ky2 = dt*dfYdT(told + dt/2, xold + kx1/2, yold + ky1/2, zold + kz1/2);
        double kz2 = dt*dfZdT(told + dt/2, xold + kx1/2, yold + ky1/2, zold + kz1/2);

        double kx3 = dt*dfXdT(told + dt/2, xold + kx2/2, yold + ky2/2, zold + kz2/2);
        double ky3 = dt*dfYdT(told + dt/2, xold + kx2/2, yold + ky2/2, zold + kz2/2);
        double kz3 = dt*dfZdT(told + dt/2, xold + kx2/2, yold + ky2/2, zold + kz2/2);

        double kx4 = dt*dfXdT(told + dt/2, xold + kx3, yold + ky3, zold + kz3);
        double ky4 = dt*dfYdT(told + dt/2, xold + kx3, yold + ky3, zold + kz3);
        double kz4 = dt*dfZdT(told + dt/2, xold + kx3, yold + ky3, zold + kz3);

        xnew = xold + (kx1 + kx2*2 + kx3*2 + kx4)/6;
        ynew = yold + (ky1 + ky2*2 + ky3*2 + ky4)/6;
        znew = zold + (kz1 + kz2*2 + kz3*2 + kz4)/6;

        output(xold, yold, zold, xnew, ynew, znew);

        xold = xnew;
        yold = ynew;
        zold = znew;

        told = t0 + dt*i;
    }
}

int main(void) {
    mySin = sin(-M_PI/4.0);
    myCos = cos(-M_PI/4.0);

    start();
    rulers();

    strokeWeight(0.5);

    translate(25, 50);
    scale(0.75, 0.75);

    rungekutta(0, 200, 40000, 1, 1, 1);

    finish();
    return 0;
}