/*
* apollonian_gasket_exterior.c
*/
#include
// "Compile with -std=c99; link with -lm."
#include
#define maximum_iterations 255
typedef struct {
double x;
double y;
} point;
typedef struct {
point center;
double radius_squared;
} circle;
/*
Below are some sample "trestleworks". To make your own, you need to define
number_of_mirrors, mirror, offset, scale, and fold().
The "mirrors" are circles that generate the fractal. (That's right, the
trestlework of mirrors supports the gasket. Just be careful not to mix your
metaphors.)
"number_of_mirrors" has to be a macro, because we initialize a static
variable with it later.
circle mirror[] is an array of circles. They should not overlap, and each one
should be tangent to at least three others. (The gasket lies entirely
inside the mirrors, except it intersects them at the tangent points. It's
perpendicular to them at those points.)
offset and scale should be the coordinates and zoom level you'd need to get
the whole gasket (or one image of the quotient space, if your gasket fills
the plane) in frame. All user-supplied values will be added to these, for
convenience. (This isn't strictly necessary -- scaling and translating the
generators would have the same effect -- but it makes it more convenient to
write some trestleworks.)
The fold function defines an equivalence relation on (x, y). By tiling or
rotating, for example, you can get away with fewer than four circles.
*/
// This trestlework is a lattice of unit circles tiling the plane.
/*
// Thanks to fold(), we only need to define one mirror.
//const unsigned int number_of_mirrors = 1;
#define number_of_mirrors 1
circle mirror[] = {{{0, 0}, 1}};
// The quotient space looks like [-1, 1)^2, with the topology of a torus.
const point offset = {0, 0};
const double scale = 0;
void fold(double *x, double *y) {
*x = remainder(*x, 2);
*y = remainder(*y, 2);
}
*/
/*
// This is the same thing with smaller circles crammed into the gaps.
//const unsigned int number_of_mirrors = 2;
#define number_of_mirrors 3
circle mirror[] = {
{{0, 0}, 1},
// {{1, 1}, 3 - sqrt(8)},
{{1, 1}, .17157287525380990240}
};
const point offset = {0, 0};
const double scale = 0;
// The quotient space looks like [0, 1]^2.
void fold(double *x, double *y) {
*x = fabs(remainder(*x, 2));
*y = fabs(remainder(*y, 2));
}
*/
// This is a classic four-circle gasket based on the smallest
// Pythagorean triple.
//const unsigned int number_of_mirrors = 2;
#define number_of_mirrors 4
#define square(x) (x*x)
// Scale by 32 and translate, so "0 0 0" gives us (-9, 55)^2, which
// includes the entire gasket.
const point offset = {23, 23};
const double scale = -5;
circle mirror[] = {
{{92, 0}, square(69)},
{{ 0, 69}, square(46)},
{{ 0, 0}, square(23)},
{{20, 21}, square( 6)}
};
// The quotient space is the whole plane.
void fold(double *x, double *y) {}
// End of generator definitions
// If (x, y) is inside mirror[m] and not at its center, map it to its
// reflection and return true. Return false otherwise.
double norm_squared(double x, double y) {
return x*x + y*y;
}
int reflect_in(double *x, double *y, unsigned int m) {
circle M = mirror[m];
point center = M.center;
double Tx = *x - center.x;
double Ty = *y - center.y;
double radius_squared = M.radius_squared;
// It's possible, if Tx and Ty are close enough to zero, for n
// to round off to zero even when Tx and Ty aren't both zero.
double ns = norm_squared(Tx, Ty);
if (ns >= radius_squared || ns == 0)
return 0;
// (x, y) is strictly inside this circle, and not at the center.
double scale = radius_squared/ns;
*x = scale*Tx + center.x;
*y = scale*Ty + center.y;
return 1;
}
// If (x, y) is inside a mirror and not at its center, map it to its
// reflection and return true. Return false otherwise.
int reflect(double *x, double *y) {
// Remember which mirror we last used.
static unsigned int m = number_of_mirrors;
unsigned int previous_mirror = m;
// Map (x, y) to a representative point of its equivalence
// class in the quotient space.
fold(x, y);
// Find the mirror, if any, containing (x, y), and calculate
// the reflection.
// If fold() didn't change (x, y), then it's a waste of time
// to try previous_mirror; every reflection is its own inverse.
// We could test this on every iteration, but I think it's
// cheaper to just save previous_mirror for last.
++m;
while (m < number_of_mirrors) {
if (reflect_in(x, y, m)) {
return 1;
}
++m;
}
m = 0;
while (m <= previous_mirror) {
if (reflect_in(x, y, m)) {
return 1;
}
++m;
}
// There's no exterior reflection of (x, y).
// We'll pass a new (x, y) next time; try all the mirrors on it.
m = number_of_mirrors;
return 0;
}
// Reflect (x, y) until it's on or outside all the mirrors.
// Return the number of reflections required.
unsigned int reflections(double x, double y) {
unsigned int count = 0;
while (reflect(&x, &y) && count < maximum_iterations) {
++count;
}
return 2*count;
}
void write_pgm(
double first_x, double last_x, double first_y, double last_y,
unsigned int rows, unsigned int columns
) {
// The PGM header:
printf(
"P5\n# (%g, %g)x(%g, %g)\n%d %d\n%d\n",
first_x, last_x, first_y, last_y,
columns, rows, maximum_iterations
);
// These can be negative, which is fine.
double width = last_x - first_x;
double height = last_y - first_y;
unsigned int half_rows = 2*rows, half_columns = 2*columns;
unsigned int half_row = 1;
while (half_row < half_rows) {
double y = first_y + height*half_row/half_rows;
unsigned int half_column = 1;
while (half_column < half_columns) {
double x = first_x + width*half_column/half_columns;
// printf("%d ", reflections(x, y));
printf("%c", (unsigned char)reflections(x, y));
half_column += 2;
}
// printf("\n");
half_row += 2;
}
}
void usage(void) {
puts("ARGs!");
}
int main(int argc, char** argv) {
double first_x, last_x, first_y, last_y;
double x, y, zoom;
if ( // center.x center.y zoom
argc == 4 &&
sscanf(argv[1], "%lf", &x) == 1 &&
sscanf(argv[2], "%lf", &y) == 1 &&
sscanf(argv[3], "%lf", &zoom) == 1
) {
x += offset.x;
y += offset.y;
double radius = exp2(-(scale + zoom));
first_x = x - radius;
last_x = x + radius;
first_y = y + radius;
last_y = y - radius;
} else if ( // left right top bottom
argc == 5 &&
sscanf(argv[1], "%lf", &first_x) == 1 &&
sscanf(argv[2], "%lf", &last_x) == 1 &&
sscanf(argv[3], "%lf", &first_y) == 1 &&
sscanf(argv[4], "%lf", &last_y) == 1
) {
} else { // Loser.
usage();
return 1;
}
write_pgm(first_x, last_x, first_y, last_y, 256, 256);
// write_pgm(0, 1, 1, 0, 256, 256);
return 0;
}