/*
    Spiral Fractal Experiment for Alf
    By: Chris M. Thomasson
    Version: Pre-Alpha (0.0.0)


    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <https://www.gnu.org/licenses/>.

*_____________________________________________________________*/


#include <cstdlib>
#include <cstdio>
#include <cassert>
#include <complex>
#include <cmath>


#define CT_WIDTH (1920 * 1)
#define CT_HEIGHT (1080 * 1)

#define CT_PI 3.1415926535897932384626433832795
#define CT_PI2 (CT_PI * 2)


typedef std::complex<double> ct_complex;


struct ct_rgb_meta
{
    double red;
    double blue;
    double green;
    double alpha;
};


#define CT_RGB_META(mp_red, ct_blue, ct_green, ct_alpha) \
    { (mp_red), (ct_blue), (ct_green), (ct_alpha) }



struct ct_rgb
{
    unsigned char r;
    unsigned char g;
    unsigned char b;
    struct ct_rgb_meta meta; // meta data for this color
};


#define CT_RGB(mp_red, mp_green, mp_blue) \
    { (mp_red), (mp_green), (mp_blue), CT_RGB_META(0., 0., 0., 0.) }


struct ct_canvas
{
    unsigned long width;
    unsigned long height;
    struct ct_rgb* buf;
};

bool
ct_canvas_create(
    struct ct_canvas* const self,
    unsigned long width,
    unsigned long height
) {
    std::size_t size = width * height * sizeof(*self->buf);

    self->buf = (ct_rgb*)std::calloc(1, size);

    if (self->buf)
    {
        self->width = width;
        self->height = height;

        return true;
    }

    return false;
}

void
ct_canvas_destroy(
    struct ct_canvas const* const self
) {
    std::free(self->buf);
}



double
ct_canvas_log_density_post_processing_get_largest(
    struct ct_canvas const* const self
) {
    double largest = 0;

    std::size_t size = self->width * self->height;

    for (std::size_t i = 0; i < size; ++i)
    {
        struct ct_rgb const* c = self->buf + i;

        if (largest < c->meta.alpha)
        {
            largest = c->meta.alpha;
        }
    }

    return largest;
}


void
ct_canvas_log_density_post_processing(
    struct ct_canvas* const self
) {
    double largest = ct_canvas_log_density_post_processing_get_largest(self);
    std::size_t size = self->width * self->height;

    printf("largest = %f\n", largest);

    for (std::size_t i = 0; i < size; ++i)
    {
        struct ct_rgb* color = self->buf + i;

        if (color->meta.alpha > 0)
        {
            struct ct_rgb_meta c = color->meta;

            double dense = log10(c.alpha) / c.alpha;

            c.red *= dense;
            if (c.red > 1.0) c.red = 1.0;

            c.green *= dense;
            if (c.green > 1.0) c.green = 1.0;

            c.blue *= dense;
            if (c.blue > 1.0) c.blue = 1.0;

            color->r = (unsigned char)(c.red * 255);
            color->g = (unsigned char)(c.green * 255);
            color->b = (unsigned char)(c.blue * 255);
        }
    }
}




bool
ct_canvas_save_ppm(
    struct ct_canvas const* const self,
    char const* fname
) {
    std::FILE* fout = std::fopen(fname, "wb");

    if (fout)
    {
        char const ppm_head[] =
            "P6\n"
            "# Chris M. Thomasson Simple 2d Plane ver:0.0.0.0 (pre-alpha)";

        std::fprintf(fout, "%s\n%lu %lu\n%u\n",
            ppm_head,
            self->width, self->height,
            255U);

        std::size_t size = self->width * self->height;

        for (std::size_t i = 0; i < size; ++i)
        {
            struct ct_rgb* c = self->buf + i;
            std::fprintf(fout, "%c%c%c", c->r, c->g, c->b);
        }

        if (!std::fclose(fout))
        {
            return true;
        }
    }

    return false;
}




struct ct_axes
{
    double xmin;
    double xmax;
    double ymin;
    double ymax;
};

struct ct_axes
    ct_axes_from_point(
        ct_complex z,
        double radius
    ) {
    struct ct_axes axes = {
        z.real() - radius, z.real() + radius,
        z.imag() - radius, z.imag() + radius
    };

    return axes;
}


struct ct_plane
{
    struct ct_axes axes;
    double xstep;
    double ystep;
};


void
ct_plane_init(
    struct ct_plane* const self,
    struct ct_axes const* axes,
    unsigned long width,
    unsigned long height
) {
    self->axes = *axes;

    double awidth = self->axes.xmax - self->axes.xmin;
    double aheight = self->axes.ymax - self->axes.ymin;

    assert(width > 0 && height > 0 && awidth > 0.0);

    double daspect = fabs((double)height / width);
    double waspect = fabs(aheight / awidth);

    if (daspect > waspect)
    {
        double excess = aheight * (daspect / waspect - 1.0);
        self->axes.ymax += excess / 2.0;
        self->axes.ymin -= excess / 2.0;
    }

    else if (daspect < waspect)
    {
        double excess = awidth * (waspect / daspect - 1.0);
        self->axes.xmax += excess / 2.0;
        self->axes.xmin -= excess / 2.0;
    }

    self->xstep = (self->axes.xmax - self->axes.xmin) / width;
    self->ystep = (self->axes.ymax - self->axes.ymin) / height;
}



struct ct_plot
{
    struct ct_plane plane;
    struct ct_canvas* canvas;
};


void
ct_plot_init(
    struct ct_plot* const self,
    struct ct_axes const* axes,
    struct ct_canvas* canvas
) {
    ct_plane_init(&self->plane, axes, canvas->width - 1, canvas->height - 1);
    self->canvas = canvas;
}


bool
ct_plot_point(
    struct ct_plot* const self,
    ct_complex z,
    struct ct_rgb const* color
) {
    long x = (long)((z.real() - self->plane.axes.xmin) / self->plane.xstep);
    long y = (long)((self->plane.axes.ymax - z.imag()) / self->plane.ystep);

    if (x > -1 && x < (long)self->canvas->width &&
        y > -1 && y < (long)self->canvas->height)
    {
        // Now, we can convert to index.
        std::size_t i = x + y * self->canvas->width;

        assert(i < self->canvas->height * self->canvas->width);

        self->canvas->buf[i] = *color;
        return true;
    }

    return false;
}



bool
ct_plot_add(
    struct ct_plot* const self,
    ct_complex z,
    struct ct_rgb const* color
) {
    long x = (long)((z.real() - self->plane.axes.xmin) / self->plane.xstep);
    long y = (long)((self->plane.axes.ymax - z.imag()) / self->plane.ystep);

    if (x > -1 && x < (long)self->canvas->width &&
        y > -1 && y < (long)self->canvas->height)
    {
        // Now, we can convert to index.
        std::size_t i = x + y * self->canvas->width;

        assert(i < self->canvas->height * self->canvas->width);

        struct ct_rgb* exist = &self->canvas->buf[i];

        exist->meta.red += color->meta.red;
        exist->meta.green += color->meta.green;
        exist->meta.blue += color->meta.blue;
        exist->meta.alpha += 1.0;

        return true;
    }

    return true;
}


// slow, so what for now... ;^)
void ct_circle(
    struct ct_plot* const plot,
    ct_complex c,
    double radius,
    unsigned int n
) {
    double abase = CT_PI2 / n;

    for (unsigned int i = 0; i < n; ++i)
    {
        double angle = abase * i;

        ct_complex z = {
            c.real() + cos(angle) * radius,
            c.imag() + sin(angle) * radius
        };

        struct ct_rgb rgb = { 255, 255, 255 };

        ct_plot_point(plot, z, &rgb);
    }
}


void ct_line(
    struct ct_plot* const plot,
    ct_complex p0,
    ct_complex p1,
    unsigned int n
) {
    ct_complex dif = p1 - p0;

    double normal_base = 1. / n;

    for (unsigned int i = 0; i < n; ++i)
    {
        double normal = normal_base * i;

        ct_complex z = p0 + dif * normal;

        struct ct_rgb rgb = { 255, 255, 0 };

        ct_plot_point(plot, z, &rgb);
    }
}





/* Vector Field
   By: Chris M. Thomasson
______________________________________________________*/
#include <vector>

struct ct_vfield_point
{
    ct_complex p;
    double m;
    ct_rgb color;
};

typedef std::vector<ct_vfield_point> ct_vfield_vector;


ct_complex
ct_normalize(
    ct_complex p
) {
    double dis = std::abs(p);

    if (!std::isnan(dis) && dis > 0.0000001)
    {
        p /= dis;
    }

    return p;
}

ct_complex
ct_vfield_normal(
    ct_vfield_vector const& self,
    ct_complex p,
    double npow
) {
    ct_complex g = { 0.0, 0.0 };

    const int imax = self.size();

    for (int i = 0; i < imax; ++i)
    {
        ct_complex dif = self[i].p - p;

        double sum = dif.real() * dif.real() + dif.imag() * dif.imag();
        double mass = std::pow(sum, npow);

        if (!std::isnan(mass) && mass > 0.000001)
        {
            g.real(g.real() + self[i].m * dif.real() / mass);
            g.imag(g.imag() + self[i].m * dif.imag() / mass);
        }
    }

    return ct_normalize(g);
}


double ct_pi_normal(
    ct_complex p,
    double start_angle = 0.0
) {
    double angle = std::arg(p) + start_angle;
    if (angle < 0.0) angle = angle + CT_PI2;
    angle = angle / CT_PI2;
    return angle;
}



void ct_vfield_plot_line(
    ct_vfield_vector const& self,
    struct ct_plot* const plot,
    ct_complex p0,
    double power,
    double astart,
    unsigned int n,
    ct_rgb const& color
) {
    double radius = .028; // 1 / (n - 0.0) * CT_PI;// .01;
    ct_complex z = p0;

    for (unsigned int i = 0; i < n; ++i)
    {
        //double ir = i / (n - 1.0);

        ct_complex vn = ct_vfield_normal(self, z, power);


        double angle = astart + arg(vn);

        //double xxxx = ct_pi_normal(angle, 0);

        ct_complex zn = {
            z.real() + std::cos(angle) * radius,
            z.imag() + std::sin(angle) * radius
        };

        //ct_rgb color_angle = color;// CT_RGBF(1, 1, 1);

       // plot.line(z, zn, 1.9, color_angle);

        ct_line(plot, z, zn, 128);

        z = zn;
    }
}


void ct_vfield_plot_line_para(
    ct_vfield_vector const& self,
    struct ct_plot* const plot,
    ct_complex p0,
    double power,
    double astart,
    double astart_radius,
    unsigned int n,
    ct_rgb lcolor
) {
    double abase = CT_PI2 / n;
    double radius = astart_radius;

    //ct_rgb bcolor = { 1, 0, 0, { 0, 0, 0 } };

    ct_complex xxxSTART = ct_vfield_normal(self, p0, power);
    double angleXXX = std::arg((p0 + p0 * xxxSTART) - p0);

    for (unsigned int i = 0; i < n; ++i)
    {
        //double ir = i / (n - 1.);

        double angle = abase * i + .0001 + angleXXX;

        ct_complex p = {
            std::cos(angle) * radius + p0.real(),
            std::sin(angle) * radius + p0.imag()
        };

        ct_rgb xcolor = lcolor;

        ct_line(plot, p0, p, 256);

        ct_vfield_plot_line(self, plot, p, power, astart, 1024, xcolor);
    }
}


// Plotter Lines
void ct_vfield_plot_lines(
    ct_vfield_vector const& self,
    struct ct_plot* const plot,
    double power,
    double astart,
    double astart_radius,
    unsigned int n
) {
    std::printf("\n\nct_vfield_plot_lines plotting...\n\n");

    for (unsigned int i = 0; i < self.size(); ++i)
    {
        //double ir = i / (self.size() - 1.);
        if (self[i].m > 0.0) continue;
        ct_rgb lcolor = self[i].color;
        ct_vfield_plot_line_para(self, plot, self[i].p, power, astart, astart_radius, n, lcolor);
        std::printf("line:(%u of %llu)\r", i, (unsigned long long)self.size() - 1);
    }
}








/* Spiral Fractal
   By: Chris M. Thomasson
______________________________________________________*/

void ct_spiral_fractal(
    unsigned int ri,
    unsigned int rn,
    struct ct_plot* const plot,
    ct_vfield_vector& vfield,
    ct_complex p0,
    ct_complex p1,
    unsigned int n
) {
    if (ri >= rn) return;

    ct_complex dif = p1 - p0;

    double sangle = std::arg(dif);
    double sradius = std::abs(dif);

    double rbase = sradius / n;
    double abase = CT_PI / 2 / n;

    for (unsigned long i = 0; i < n; ++i)
    {
        double radius = rbase * i;

        double angle = 0 + std::sin(abase * i) * CT_PI2 + sangle;
        double anglenx = 0 + std::sin(abase * (i + 1)) * CT_PI2 + sangle;

        ct_complex s0_raw = { std::cos(angle), std::sin(angle) };
        ct_complex s1_raw = { std::cos(anglenx), std::sin(anglenx) };

        ct_complex s0 = p0 + s0_raw * radius;
        ct_complex s1 = p0 + s1_raw * (radius + rbase);

        ct_rgb color = { 0, 0, 0, { .7, .6, .5 } };

        ct_plot_add(plot, s0, &color);
        ct_plot_add(plot, s1, &color);

        //ct_plot_add(plot, -s0, &color);
       // ct_plot_add(plot, -s1, &color);


        if (ri == rn - 1)
        {
            //ct_line(plot, s0, s1, 128);

            vfield.push_back({ s0, -1, { 1, 0, 0, { 0, 0, 0 } } });
            //vfield.push_back({ s1, -1, { 1, 0, 0, { 0, 0, 0 } } });
        }


        if (!((i + 1) % 2))
        {
            ct_spiral_fractal(ri + 1, rn, plot, vfield, s0, s1, n * 2);
        }

        else
        {
            ct_spiral_fractal(ri + 1, rn, plot, vfield, s1, s0, n * 2);
        }
    }
}




/* Entry
______________________________________________________*/
int main(void)
{
    std::printf("\n\nChris M. Thomasson's Vector Field based on Spiral Plot\n\n");


    struct ct_canvas canvas;

    if (ct_canvas_create(&canvas, CT_WIDTH, CT_HEIGHT))
    {
        ct_complex plane_origin = { 0, 0 };
        double plane_radius = 2.0;

        struct ct_axes axes = ct_axes_from_point(plane_origin, plane_radius);

        struct ct_plot plot;
        ct_plot_init(&plot, &axes, &canvas);

        // Our unit circle
        ct_circle(&plot, { 0, 0 }, 1.0, 2048);
        ct_circle(&plot, { 2, 0 }, 2.0, 2048);
        ct_circle(&plot, { -2, 0 }, 2.0, 2048);
        ct_circle(&plot, { 0, 2 }, 2.0, 2048);
        ct_circle(&plot, { 0, -2 }, 2.0, 2048);

        {
            ct_vfield_vector vfield;

            std::printf("\n\nct_spiral_fractal plotting...\n\n");
            ct_spiral_fractal(0, 3, &plot, vfield, { -1, 0 }, { 1, 0 }, 3);

            ct_vfield_plot_lines(vfield, &plot, 2., 1.2, 0.04, 6);
        }

        ct_canvas_save_ppm(&canvas, "ct_plane.ppm");

        ct_canvas_destroy(&canvas);

        std::printf("\n\nPlot Complete! :^)\n\n");

        return EXIT_SUCCESS;
    }

    return EXIT_FAILURE;
}