Complex Iteration Example in C

/* Crude, simple code for Sparky
   Showing how to code complex numbers. This implments:
   (Z^2 + ((RE(Z) + IM(Z)) / 2.1)) - 0.83
_____________________________________________________________*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <bmpfile.h>



static bmpfile_t* g_bmp = NULL;




/* A 2d point
_____________________________________________________________*/
struct ffe_point
{
    double x;
    double y;
};

struct ffe_point
ffe_complex_add(
    struct ffe_point n0,
    struct ffe_point n1
){
    struct ffe_point ret = { n0.x + n1.x, n0.y + n1.y };
    return ret;
}

struct ffe_point
ffe_complex_add_real(
    struct ffe_point n0,
    double real
){
    struct ffe_point ret = { n0.x + real, n0.y };
    return ret;
}

struct ffe_point
ffe_complex_sub_real(
    struct ffe_point n0,
    double real
){
    struct ffe_point ret = { n0.x - real, n0.y };
    return ret;
}

struct ffe_point
ffe_complex_mul(
    struct ffe_point n0,
    struct ffe_point n1
){
    struct ffe_point ret = {
        n0.x * n1.x - n0.y * n1.y,
        n0.x * n1.y + n0.y * n1.x
    };

    return ret;
}




/* Axes describing a 2d plane
_____________________________________________________________*/
struct ffe_axes
{
    double xmin;
    double xmax;
    double ymin;
    double ymax;
};



/* Complex Iteration
_____________________________________________________________*/
void iterate_point(struct ffe_point const* c, unsigned int x, unsigned int y)
{
    struct ffe_point z = *c;

    /*(Z^2 + ((RE(Z) + IM(Z)) / 2.1)) - 0.83*/

    for (unsigned int i = 0; i < 256; ++i)
    {
        double temp = ((z.x + z.y) / 2.1);

        z = ffe_complex_mul(z, z);
        z = ffe_complex_add_real(z, temp);
        z = ffe_complex_sub_real(z, 0.83);

        double dis = sqrt(z.x * z.x + z.y * z.y);

        if (dis > 10.0)
        {
            rgb_pixel_t color = { 255, ((i + 1) * 5) % 255, 0, 0 };
            bmp_set_pixel(g_bmp, x, y, color);
            return;
        }
    }

    rgb_pixel_t color = { 255, 0, 0, 0 };
    bmp_set_pixel(g_bmp, x, y, color);
}


void iterate_axes(struct ffe_axes const* axes, unsigned int width, unsigned int height)
{
    double xstep = (axes->xmax - axes->xmin) / (width - 1.0);
    double ystep = (axes->ymax - axes->ymin) / (height - 1.0);

    for (unsigned int y = 0; y < height; ++y)
    {
        for (unsigned int x = 0; x < width; ++x)
        {
            struct ffe_point c = {
                axes->xmin + (x * xstep),
                axes->ymin + (y * ystep)
            };

            iterate_point(&c, x, y);
        }
    }
}


int
main(void)
{
    #define OUTBMP "C:\\Users\\Chris\\Desktop\\hello_fractal.bmp"

    {
        unsigned int width = 300;
        unsigned int height = 300;

        g_bmp = bmp_create(width, height, 24);

        struct ffe_axes axes = { -2.0, 2.0, -2.0, 2.0 };

        iterate_axes(&axes, width, height);

        bmp_save(g_bmp, OUTBMP);

        bmp_destroy(g_bmp);
    }

    puts("Done, hit <ENTER> to exit...");
    getchar();

    return 0;
}