Julia set rendering code

/*
 * Asks for a value of c, then prints the filled-in Julia set of f(z) = z^2 + c to the Windows terminal.
 */

#define _CRT_SECURE_NO_WARNINGS

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

// The resolution of the grid. Make sure the terminal is wide enough to print this many characters.
#define X_PIX 620
#define Y_PIX 310

// The maximum number of iterations to run.
// I find that increasing this value too greatly makes the output skinnier.
#define N_MAX 100

// Struct for complex numbers, in rectangular form. z = re + im * I.
struct Complex {
    double re;
    double im;
};

// Squares a complex number.
struct Complex square_complex(struct Complex z) {
    struct Complex output = { z.re * z.re - z.im * z.im, 2.00 * z.re * z.im };
    return output;
}

// Calculates the modulus of a complex number.
double modulus(struct Complex z) {
    return sqrt(z.re * z.re + z.im * z.im);
}

// Iterates a complex number through the function f(z) = z^2 + c.
struct Complex iterate(struct Complex z, struct Complex c) {
    struct Complex output = square_complex(z);
    output.re += c.re;
    output.im += c.im;
    return output;
}

int main(void) {
    // Initialize c and prompt the user for its real and imaginary parts.
    struct Complex c;

    printf("Enter the real part of c: ");
    scanf("%lf", &c.re);

    printf("\nEnter the imaginary part of c: ");
    scanf("%lf", &c.im);

    // The value rad is R, the threshold for determining whether the sequence |A_n| is bounded.
    double rad = 0.5 + 0.5 * sqrt(1 + 4 * modulus(c));

    printf("\n");

    double re, im;

    // Set up a nested for loop.
    for (int y = 1; y < Y_PIX; y++) {
        for (int x = 1; x < X_PIX; x++) {

            // Calculate the real and imaginary parts of the complex plane at the point (x, y) according to the grid resolution.
            // The grid displayed is from -rad to rad on the real and imaginary axes. Points outside of this area are guaranteed to not be in the filled-in Julia set.
            re = rad * -1.0 + ((rad *  2.0) / (double)X_PIX * x);
            im = rad *  1.0 + ((rad * -2.0) / (double)Y_PIX * y);

            struct Complex z = { re, im };

            // Iterate z through the function f(z) = z^2 + c up to n_MAX times, stopping as soon as we find a term whose distance from the origin is further than rad.
            int iter = 0;

            while (iter < N_MAX && modulus(z) <= rad) {
                z = iterate(z, c);

                iter++;
            }

            if (iter == N_MAX) {
                // The point is likely in the filled-in Julia set. Print a $.
                printf("$");
            }
            else {
                // The point is not in the filled-in Julia set. Print a period.
                printf(".");
            }

        }
        printf("\n");
    }

    system("PAUSE");
    return 0;
}