Multi-Julia Buddha Plotter by Chris M. Thomasson

/*
    Multi-Julia 2-ary Buddha Plotter
    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 <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <complex.h>
#include <tgmath.h>
#include <stdbool.h>


#define CT_WIDTH 1920
#define CT_HEIGHT 1080
#define CT_N 24000000



struct ct_rgb
{
    unsigned char r;
    unsigned char g;
    unsigned char b;
};


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
){
    size_t size = width * height * sizeof(*self->buf);

    self->buf = 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
){
    free(self->buf);
}

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

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

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

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

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

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

    return false;
}




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

struct ct_axes
ct_axes_from_point(
    double complex z,
    double radius
){
    struct ct_axes axes = {
        creal(z) - radius, creal(z) + radius,
        cimag(z) - radius, cimag(z) + 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_add(
    struct ct_plot* const self,
    double complex z,
    struct ct_rgb const* color
){
    long x = (creal(z) - self->plane.axes.xmin) / self->plane.xstep;
    long y = (self->plane.axes.ymax - cimag(z)) / self->plane.ystep;

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

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

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

        unsigned long addend = 3;

        // Sorry for the if ladder! Need to refine this... Yikes!
        if (exist.r < 255 - addend)
        {
             exist.r += addend;
        }

        else
        {
            exist.r = 255;

            if (exist.g < 255 - addend)
            {
                 exist.g += addend;
            }

            else
            {
                exist.g = 255;

                if (exist.b < 255 - addend)
                {
                     exist.b += addend;
                }

                else
                {
                    exist.b = 255;
                }
            }
        }

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

    return true;
}


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

    if (x > -1 && x < (long)self->canvas->width &&
        y > -1 && y < (long)self->canvas->height)
    {
        // Now, we can convert to index.
        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;
}


// The All Plotting Buddha
double complex
ct_all_plotting_buddha(
    struct ct_plot* const plot,
    double complex z,
    double complex c,
    unsigned long n
){
    //double complex origin_z = z;

    for (unsigned long i = 0; i < n; ++i)
    {
        z = z * z + c; // Mbrot formula

        // plot point
        struct ct_rgb rgb = { 255, 0, 0 };

        ct_plot_add(plot, z, &rgb);
    }

    return z;
}



// Compute the fractal
void
ct_ifs(
    struct ct_plot* const plot,
    unsigned long n
){
    // 2 sets

    /*
    // Magic Carpet
    double complex jp[] = {
        0 + 1 * I,
        0 + -1 * I
    };
    */

    double complex jp[] = {
        0.726514 + 0.240242 * I,
        0.171039 + 0.235043 * I
    };

    double complex z = 0+0*I;
    double complex c = 0+0*I;

    for (unsigned long i = 0; i < n; ++i)
    {
        double rn0 = rand() / (RAND_MAX - .0);
        double rn1 = rand() / (RAND_MAX - .0);

        if (rn0 > .5)
        {
            c = jp[0];
        }

        else
        {
            c = jp[1];
        }


        double complex d = z - c;
        double complex root = csqrt(d);

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

        ct_plot_add(plot, z, &rgb);
        ct_plot_add(plot, -z, &rgb);
        ct_plot_add(plot, root, &rgb);
        ct_plot_add(plot, -root, &rgb);

        ct_all_plotting_buddha(plot, z, z, 3);

        z = root;

        if (rn1 > .5)
        {
            z = -root;
        }

        if (! (i % (n / 3)))
        {
            printf("rendering: %lu of %lu\r", i + 1, n);
        }
    }

    printf("rendering: %lu of %lu\n", n, n);
}





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

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

        double complex z =
            (creal(c) + cos(angle) * radius) +
            (cimag(c) + sin(angle) * radius) * I;

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

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


int main(void)
{
    struct ct_canvas canvas;

    if (ct_canvas_create(&canvas, CT_WIDTH, CT_HEIGHT))
    {
        double complex plane_origin = 0+0*I;
        double plane_radius = 1.45;

        struct ct_axes axes = ct_axes_from_point(plane_origin, plane_radius);

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


        ct_ifs(&plot, CT_N);

        //ct_para(&plot, plane_origin, 1.0);

        // Our unit circle
        ct_circle(&plot, 0+0*I, 1.0, 2048);
        ct_circle(&plot, 2+0*I, 2.0, 2048);
        ct_circle(&plot, -2+0*I, 2.0, 2048);
        ct_circle(&plot, 0+2*I, 2.0, 2048);
        ct_circle(&plot, 0-2*I, 2.0, 2048);

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

        ct_canvas_destroy(&canvas);

        return EXIT_SUCCESS;
    }

    return EXIT_FAILURE;
}