A PPM plotter with aspect ratios...

/*
    A Simple 2d Plane For Hugh, with proper aspect ratios!
    By: Chris M. Thomasson
*_____________________________________________________________*/


#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 5000000


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

bool
ct_canvas_create(
    struct ct_canvas* const self,
    unsigned long width,
    unsigned long height
){
    size_t size = width * height;

    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, "w");

    if (fout)
    {
        char const ppm_head[] =
            "P3\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)
        {
            unsigned int c = self->buf[i];
            fprintf(fout, "%u %u %u  ", c, 0U, 0U);
        }

        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_point(
    struct ct_plot* const self,
    double complex z,
    unsigned char 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;
}



// Compute the fractal
void
ct_ifs(
    struct ct_plot* const plot,
    unsigned long n
){
    // 2 sets
    double complex jp[] = {
        .5 + 0 * I,
        -5.5 + I*.0
    };

    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);

        z = root;

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

         ct_plot_point(plot, z, 255);
         ct_plot_point(plot, root, 255);

        if (! (i % 256))
        {
            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;

        ct_plot_point(plot, z, 255);
    }
}


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 = 2.0;

        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);

        // 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;
}