/* dala2 --- plot an HPGL mandala 2011-09-10 */

/* Copyright (c) 2011 John Honniball, Froods Software Development */

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#include <math.h>

#include <unistd.h>

#include "hpgllib.h"

void catcher(int sig);

/* Central disc patterns */

double lissajous(const double x0, const double y0, double side, const double f1, const double f2, const double theta, const int npts);

double spiraldisc(const double x0, const double y0, const double radius, const int n, const int turns);

double sunflowerspiral(const double x0, const double y0, const double radius);

double lotus(const double x0, const double y0, double radius, const int nleaves, const int nlayers);

double diamondsquare(const double x0, const double y0, double side);

/* Intermediate ring patterns */

double zigzagring(const double x0, const double y0, const double r1, const double r2, const int npts, const int incr, const int flag);

double radials(const double x0, const double y0, const double inner, const double length, const int n);

double ringoboxes(const double x0, const double y0, const double radius, const int nboxes, const int ninner);

double circlediamonds(const double x0, const double y0, double radius, const int ndiam);

double ringofcircles(const double x0, const double y0, double radius, const int ncirc);

double sqwavering(const double x0, const double y0, const double radius, const double len, int nwaves);

double ellipsering(const double x0, const double y0, const double a, const double b, const int nell);

double ringofpoly(const double x0, const double y0, const double radius, const int nside, const int npoly);

double circlering(const double x0, const double y0, const double r1, const double r2, const int n);

/* Outer circle-to-square transitions */

double circletosuper(const double x0, const double y0, const double r1, const double r2, const int n);

double sqradials(const double x0, const double y0, const double r1, const double r2, const int n);

void curvestitchtrans(const double x0, const double y0, const double r1, const double r2, const int n);

/* Square border drawing functions */

void saltireborder(const double x0, const double y0, const int nx, const int ny, const double w, const double h);

void roundrectborder(const double x0, const double y0, const double r1, const double r2, const int n);

void circleborder(const double x0, const double y0, const int nx, const int ny, const double w, const double h);

void celticstepborder(const double x0, const double y0, const int nx, const int ny, const double w, const double h);

/* Utility drawing functions */

void radialfan(const double x0, const double y0, const double r1, const double startangle, const double dx, const double dy, const double xvec, const double yvec, const int n);

void spiral(const double cx, const double cy, const double radius, const double ang, const int n);

void bordergrid(const double x0, const double y0, const double w, const double h, const int nx, const int ny, const int diag);

void superellipse(const double x0, const double y0, const double a, const double b, const double theta, const double d);

double lotuszigzag(const double x0, const double y0, const double r1, const double r2, const int npts, const int flag);

void drawstep(const double x0, const double y0, const double w, const double h, const int flipx, const int flipy);

int main(int argc, char * const argv[])

{

int opt;

double xc, yc;

double maxx, maxy;

double height;

double radius;

double inner_r;

double scale = 40.0;

while ((opt = getopt(argc, argv, "no:p:s:t:v:")) != -1) {

switch (opt) {

case 's':

if (strchr(optarg, '1'))

scale = 80.0;

else if (strchr(optarg, '2'))

scale = 56.57;

else if (strchr(optarg, '4'))

scale = 28.28;

else if (strchr(optarg, '5'))

scale = 20.0;

case 'n':

case 'o':

case 'p':

case 't':

case 'v':

plotopt (opt, optarg);

break;

default: /* '?' */

fprintf(stderr, "Usage: %s [-p pen] [-s <size>] [-t title]\n", argv[0]);

fprintf(stderr, " <size> ::= A1 | A2 | A3 | A4 | A5\n");

exit(EXIT_FAILURE);

}

}

if (plotbegin(0) < 0) {

fputs("Failed to initialise HPGL library\n", stderr);

exit(EXIT_FAILURE);

}

getplotsize(&maxx, &maxy);

xc = maxx / 2.0;

yc = maxy / 2.0;

height = maxy;

/* Draw square border */

rectangle(xc - (height / 2.0), 0.0, xc + (height / 2.0), maxy);

radius = 20.0 * scale;

/* Draw central spirals */

// radius = spiraldisc(xc, yc, radius, 3, 2);

radius = lotus(xc, yc, radius, 16, 7);

// radius = lissajous(xc, yc, 38.0 * scale, 3.0, 4.0, 0.0, 144);

/* Plot ring of ellipses around spirals */

circle(xc, yc, radius);

radius = ellipsering(xc, yc, radius + (1.0 * scale), radius + (9.0 * scale), 3);

radius += 1.0 * scale;

/* Pre-draw circles around zigzags */

circle(xc, yc, radius);

circle(xc, yc, radius + (50.0 * scale));

/* Plot radial lines */

radius = radials(xc, yc, radius, 5.0 * scale, 63);

/* Plot two rings of zigzags */

radius = zigzagring(xc, yc, radius, radius + (20.0 * scale), 63, 4, 0);

radius = zigzagring(xc, yc, radius, radius + (25.0 * scale), 63, 8, 1);

/* Plot steampunk cog */

// radius = sqwavering(xc, yc, radius + (2.0 * scale), 12.0 * scale, 18);

/* Plot ring of concentric boxes */

radius = ringoboxes(xc, yc, radius + (2.0 * scale), 48, 3);

/* Plot ring of polygons */

// radius = ringofpoly(xc, yc, radius + (12.0 * scale), 7, 5);

/* Plot circle with diamonds */

circle(xc, yc, radius + (2.0 * scale));

radius = circlediamonds(xc, yc, radius + (2.0 * scale), 36);

circle(xc, yc, radius);

radius += 2.0 * scale;

/* Plot superellipses with increasingly sharp corners */

// circletosuper(xc, yc, radius, height / 2.0, 6);

// sqradials(xc, yc, radius, height / 2.0, 160);

// circlering(xc, yc, radius, height / 2.0, 4);

// ringofcircles(xc, yc, radius, 96);

inner_r = (height / 2.0) - (height / 12.0);

// circle(xc, yc, inner_r);

curvestitchtrans(xc, yc, radius, inner_r, 20);

/* Draw a square border pattern */

// roundrectborder(xc, yc, inner_r, height / 2.0, 5);

// saltireborder(xc - (height / 2.0), 0.0, 12, 12, height / 12.0, height / 12.0);

// circleborder(xc - (height / 2.0), 0.0, 12, 12, height / 12.0, height / 12.0);

celticstepborder(xc - (height / 2.0), 0.0, 12, 12, height / 12.0, height / 12.0);

plotend();

return (0);

}

/* lissajous --- draw a lissajous pattern */

double lissajous(const double x0, const double y0, double side, const double f1, const double f2, const double theta, const int npts)

{

const double delta = (2.0 * M_PI) / (double)npts;

const double sintheta = sin(theta);

const double costheta = cos(theta);

double r;

int i;

side /= 2.0;

r = side;

for (i = 0; i <= npts; i++) {

const double t = (double)i * delta;

const double t1 = t * f1;

const double t2 = t * f2;

const double x = (r * cos(t1) * costheta) - (r * sin(t2) * sintheta);

const double y = (r * cos(t1) * sintheta) + (r * sin(t2) * costheta);

if (i == 0)

moveto(x0 + x, y0 + y);

else

lineto(x0 + x, y0 + y);

}

/* Return radius of circumscribing circle */

return (sqrt((side * side) * 2.0));

}

/* spiraldisc --- draw multiple spirals from a common centre point */

double spiraldisc(const double x0, const double y0, const double radius, const int n, const int turns)

{

int i;

const double delta = (M_PI * 2.0) / (double)n;

for (i = 0; i < n; i++) {

const double theta = (double)i * delta;

spiral(x0, y0, radius, theta, turns);

}

return (radius);

}

/* diamondsquare --- draw concentric squares at 45-degree offsets */

double diamondsquare(const double x0, const double y0, double side)

{

side /= 2.0;

/* Plot outer square */

rectangle(x0 - side, y0 - side, x0 + side, y0 + side);

/* Plot inner diagonal square */

moveto(x0, y0 - side);

lineto(x0 + side, y0);

lineto(x0, y0 + side);

lineto(x0 - side, y0);

lineto(x0, y0 - side);

/* Return radius of circumscribing circle */

return (sqrt((side * side) * 2.0));

}

/* zigzagring --- draw a ring of zigzags */

double zigzagring(const double x0, const double y0, const double r1, const double r2, const int npts, const int incr, const int flag)

{

int i;

double x1[128], y1[128];

double x2[128], y2[128];

double theta1, theta2;

const double delta = (2.0 * M_PI) / (double)npts;

int n1, n2;

for (i = 0; i < npts; i++) {

if (flag) {

theta2 = delta * (double)i;

theta1 = (delta * (double)i) + (delta / 2.0);

}

else {

theta1 = delta * (double)i;

theta2 = (delta * (double)i) + (delta / 2.0);

}

x1[i] = (cos(theta1) * r1) + x0;

y1[i] = (sin(theta1) * r1) + y0;

x2[i] = (cos(theta2) * r2) + x0;

y2[i] = (sin(theta2) * r2) + y0;

}

moveto(x1[0], y1[0]);

n1 = 0;

n2 = incr / 2;

for (i = 0; i < npts; i++) {

lineto(x1[n1], y1[n1]);

lineto(x2[n2], y2[n2]);

n1 = (n1 + incr) % npts;

n2 = (n2 + incr) % npts;

}

lineto(x1[0], y1[0]);

return (r2);

}

/* radials --- draw simple radial lines */

double radials(const double x0, const double y0, const double inner, const double length, const int n)

{

int i;

const double delta = (2.0 * M_PI) / (double)n;

for (i = 0; i < n; i++) {

const double theta = (double)i * delta;

const double xvec = cos(theta);

const double yvec = sin(theta);

const double x1 = xvec * inner;

const double y1 = yvec * inner;

const double x2 = xvec * (inner + length);

const double y2 = yvec * (inner + length);

if (i & 1) {

moveto(x0 + x1, y0 + y1);

lineto(x0 + x2, y0 + y2);

}

else {

moveto(x0 + x2, y0 + y2);

lineto(x0 + x1, y0 + y1);

}

}

return (inner + length);

}

double ringoboxes(const double x0, const double y0, const double radius, const int nboxes, const int ninner)

{

int i, j, k;

double side, s2;

const double delta = (2.0 * M_PI) / (double)nboxes;

double x[4], y[4];

double rx[4], ry[4];

double inc;

side = (2.0 * M_PI * radius) / (double)nboxes;

side *= 0.8;

s2 = side / 2.0;

inc = s2 / (double)ninner;

for (i = 0; i < nboxes; i++) {

const double theta = (double)i * delta;

const double s = sin(theta);

const double c = cos(theta);

for (k = 0; k < ninner; k++) {

/* Set up a square */

x[0] = -s2;

y[0] = -s2;

x[1] = s2;

y[1] = -s2;

x[2] = s2;

y[2] = s2;

x[3] = -s2;

y[3] = s2;

/* Shrink, rotate and translate square */

for (j = 0; j < 4; j++) {

if (x[j] < 0)

x[j] += k * inc;

else

x[j] -= k * inc;

if (y[j] < 0)

y[j] += k * inc;

else

y[j] -= k * inc;

rx[j] = (x[j] * c) - (y[j] * s);

ry[j] = (x[j] * s) + (y[j] * c);

rx[j] += x0 + (c * (radius + s2));

ry[j] += y0 + (s * (radius + s2));

}

/* Draw the rotated square */

moveto(rx[0], ry[0]);

lineto(rx[1], ry[1]);

lineto(rx[2], ry[2]);

lineto(rx[3], ry[3]);

lineto(rx[0], ry[0]);

}

}

return (radius + side);

}

double circlediamonds(const double x0, const double y0, double radius, const int ndiam)

{

int i;

const double delta = (M_PI * 2.0) / (double)ndiam;

const double diag = (2.0 * M_PI * radius) / (double)ndiam;

double r1, r2;

double x1[128], y1[128];

double x2[128], y2[128];

radius += diag / 2.0;

circle (x0, y0, radius);

r1 = radius - (diag / 2.0);

r2 = radius + (diag / 2.0);

for (i = 0; i < ndiam; i++) {

const double theta = (double)i * delta;

x1[i] = (r1 * cos(theta)) + x0;

y1[i] = (r1 * sin(theta)) + y0;

x2[i] = (r2 * cos(theta)) + x0;

y2[i] = (r2 * sin(theta)) + y0;

}

moveto(x1[0], y1[0]);

for (i = 0; i < ndiam; i++) {

if (i & 1)

lineto(x2[i], y2[i]);

else

lineto(x1[i], y1[i]);

}

lineto(x1[0], y1[0]);

moveto(x2[0], y2[0]);

for (i = 0; i < ndiam; i++) {

if (i & 1)

lineto(x1[i], y1[i]);

else

lineto(x2[i], y2[i]);

}

lineto(x2[0], y2[0]);

return (r2);

}

/* ringofpoly --- draw overlapping polygons to form a ring */

double ringofpoly(const double x0, const double y0, const double radius, const int nside, const int npoly)

{

int i, j;

const double delta = (M_PI * 2.0) / (double)npoly;

const double angperside = (M_PI * 2.0) / (double)nside;

for (i = 0; i < npoly; i++) {

const double theta = (double)i * delta;

double x = radius * cos(theta);

double y = radius * sin(theta);

moveto(x0 + x, y0 + y);

for (j = 0; j < nside; j++) {

const double intang = (double)(j + 1) * angperside;

x = radius * cos(theta + intang);

y = radius * sin(theta + intang);

lineto(x0 + x, y0 + y);

}

}

return (radius);

}

double ringofcircles(const double x0, const double y0, double radius, const int ncirc)

{

const double delta = (M_PI * 2.0) / (double)ncirc;

double r2;

int i;

/* Compute radius as if centres are on circumference of 'radius' */

r2 = radius * sin(delta);

/* Increase radius by radius of small circles */

radius += r2;

/* Recompute radius of smaller circles */

r2 = radius * sin(delta);

for (i = 0; i < ncirc; i++) {

const double theta = (double)i * delta;

const double x = radius * cos(theta);

const double y = radius * sin(theta);

circle2(x0 + x, y0 + y, r2, 10.0);

}

return (radius + r2);

}

/* sqwavering --- draw a cog-like ring of square waves */

double sqwavering(const double x0, const double y0, const double radius, const double len, int nwaves)

{

double delta;

double degrees;

int i;

/* Number of half-waves; must be even */

nwaves *= 2;

delta = (2.0 * M_PI) / (double)nwaves;

degrees = 360.0 / (double)nwaves;

for (i = 0; i < nwaves; i++) {

const double theta = (double)i * delta;

const double xvec = cos(theta);

const double yvec = sin(theta);

const double x1 = xvec * radius;

const double y1 = yvec * radius;

const double x2 = xvec * (radius + len);

const double y2 = yvec * (radius + len);

if (i == 0)

moveto(x0 + x2, y0 + y2);

if (i & 1)

lineto(x0 + x2, y0 + y2);

else

lineto(x0 + x1, y0 + y1);

arc(x0, y0, degrees);

}

return (radius + len);

}

double ellipsering(const double x0, const double y0, const double a, const double b, const int nell)

{

const double delta = M_PI / (double)nell;

int i;

for (i = 0; i < nell; i++)

ellipse(x0, y0, a, b, delta * i);

if (a > b)

return (a);

else

return (b);

}

/* circletosuper --- draw transition from circle to superellipses */

double circletosuper(const double x0, const double y0, const double r1, const double r2, const int n)

{

static double pow2[6] = {

2.8, 4.0, 5.65, 8.0, 11.3, 16.0

};

const double step = (r2 - r1) / (double)n;

int i;

for (i = 0; i < n; i++) {

const double a = r1 + (i * step);

const double b = r1 + (i * step);

superellipse(x0, y0, a, b, 0.0, pow2[i]);

}

return (r2);

}

double circlering(const double x0, const double y0, const double r1, const double r2, const int n)

{

int i;

const double delta = (M_PI * 2.0) / (double)n;

const double r = (r2 - r1) / 2.0;

for (i = 0; i < n; i++) {

const double theta = (double)i * delta;

const double x = r * cos(theta);

const double y = r * sin(theta);

circle(x0 + x, y0 + y, r1 + r);

}

return (r2);

}

/* sqradials --- draw radial lines from inner circle to outer square */

double sqradials(const double x0, const double y0, const double r1, const double r2, const int n)

{

const double side = 2.0 * r2; /* Length of side of overall square */

const double r45 = M_PI / 4.0; /* 45 degrees in radians */

radialfan(x0, y0, r1, r45, r2, r2, -side, 0.0, n / 4);

radialfan(x0, y0, r1, r45 * 3.0, -r2, r2, 0.0, -side, n / 4);

radialfan(x0, y0, r1, r45 * 5.0, -r2, -r2, side, 0.0, n / 4);

radialfan(x0, y0, r1, r45 * 7.0, r2, -r2, 0.0, side, n / 4);

return (r2);

}

/* roundrectborder --- draw border as nested rounded rectangles */

void roundrectborder(const double x0, const double y0, const double r1, const double r2, const int n)

{

int i;

const double incr = (r2 - r1) / (double)n;

/* Innermost rectangle does not have rounded corners */

rectangle(x0 - r1, y0 - r1, x0 + r1, y0 + r1);

for (i = 1; i < n; i++) {

const double radius = incr * (double)i;

const double xmin = x0 - r1 - radius;

const double ymin = y0 - r1 - radius;

const double xmax = x0 + r1 + radius;

const double ymax = y0 + r1 + radius;

roundrect(xmin, ymin, xmax, ymax, radius);

}

}

/* saltireborder --- draw border squares as saltires */

void saltireborder(const double x0, const double y0, const int nx, const int ny, const double w, const double h)

{

/* nx and ny must be even, but need not be equal */

const double xmax = x0 + (nx * w);

const double ymax = y0 + (ny * h);

int i;

bordergrid (x0, y0, w, h, nx, ny, 1);

/* Zig-zag across bottom */

moveto(x0, y0 + (2.0 * h));

lineto(x0 + (2.0 * w), y0);

for (i = 3; i < (nx - 1); i += 2) {

lineto(x0 + (i * w), y0 + h);

lineto(x0 + ((i + 1) * w), y0);

}

/* Zig-zag up RHS */

lineto(xmax, y0 + (2.0 * h));

for (i = 3; i < (ny - 1); i += 2) {

lineto(xmax - w, y0 + (i * h));

lineto(xmax, y0 + ((i + 1) * h));

}

/* Zig-zag across top, R to L */

lineto(xmax - (2.0 * w), ymax);

for (i = nx - 3; i > 1; i -= 2) {

lineto(x0 + (i * w), ymax - h);

lineto(x0 + ((i - 1) * w), ymax);

}

/* Zig-zag down LHS */

lineto(x0, ymax - (2.0 * h));

for (i = ny - 3; i > 1; i -= 2) {

lineto(x0 + w, y0 + (i * h));

lineto(x0, y0 + ((i - 1) * h));

}

/* Zig-zag across bottom again, L to R */

moveto(x0, y0 + h);

lineto(x0 + w, y0);

for (i = 2; i < nx; i += 2) {

lineto(x0 + (i * w), y0 + h);

lineto(x0 + ((i + 1) * w), y0);

}

/* Zig-zag up RHS */

lineto(xmax, y0 + h);

for (i = 2; i < ny; i += 2) {

lineto(xmax - w, y0 + (i * h));

lineto(xmax, y0 + ((i + 1) * h));

}

/* Zig-zag across top, R to L */

lineto(xmax - w, ymax);

for (i = nx - 2; i > 1; i -= 2) {

lineto(x0 + (i * w), ymax - h);

lineto(x0 + ((i - 1) * w), ymax);

}

/* Zig-zag down LHS */

lineto(x0, ymax - h);

for (i = ny - 2; i > 1; i -= 2) {

lineto(x0 + w, y0 + (i * h));

lineto(x0, y0 + ((i - 1) * h));

}

}

/* circleborder --- draw border squares as circles */

void circleborder(const double x0, const double y0, const int nx, const int ny, const double w, const double h)

{

/* w and h must be the same for the circles to fit in the grid */

const double xmax = x0 + (nx * w);

const double ymax = y0 + (ny * h);

int i;

bordergrid(x0, y0, w, h, nx, ny, 0);

for (i = 1; i < nx; i++) {

circle(x0 + (w / 2.0) + (w * (double)i), h / 2.0, h / 2.0);

}

for (i = 1; i < ny; i++) {

circle(xmax - (w / 2.0), y0 + (h / 2.0) + (h * (double)i), h / 2.0);

}

for (i = nx - 2; i >= 0; i--) {

circle(x0 + (w / 2.0) + (w * (double)i), ymax - (h / 2.0), h / 2.0);

}

for (i = ny - 2; i >= 0; i--) {

circle(x0 + (w / 2.0), y0 + (h / 2.0) + (h * (double)i), h / 2.0);

}

}

/* celticstepborder --- draw a border grid of Celtic step pattern */

void celticstepborder(const double x0, const double y0, const int nx, const int ny, const double w, const double h)

{

const double xmax = x0 + (nx * w);

const double ymax = y0 + (ny * h);

int i;

bordergrid(x0, y0, w, h, nx, ny, 0);

for (i = 1; i < nx; i++) {

drawstep(x0 + (w * (double)i), y0, w, h, i & 1, 0);

}

for (i = 1; i < ny; i++) {

drawstep(xmax - w, y0 + (h * (double)i), w, h, 1, i & 1);

}

for (i = nx - 2; i >= 0; i--) {

drawstep(x0 + (w * (double)i), ymax - h, w, h, i & 1, 1);

}

for (i = ny - 2; i >= 0; i--) {

drawstep(x0, y0 + (h * (double)i), w, h, 0, i & 1);

}

}

/* curvestitchtrans --- transition by curve-stiching (pin-and-cotton) */

void curvestitchtrans(const double x0, const double y0, const double r1, const double r2, const int n)

{

int i;

double xmin, xmax;

double ymin, ymax;

const double xinc = r2 / (double)n;

const double yinc = r2 / (double)n;

/* Draw lower RH corner */

xmin = x0;

ymin = y0 - r2;

xmax = x0 + r2;

ymax = y0;

for (i = 1; i < n; i++) {

if (i & 1) {

moveto(xmin + (xinc * (double)i), ymin);

lineto(xmax, ymin + (yinc * (double)i));

}

else {

moveto(xmax, ymin + (yinc * (double)i));

lineto(xmin + (xinc * (double)i), ymin);

}

}

/* Draw upper RH corner */

xmin = x0;

ymin = y0;

xmax = x0 + r2;

ymax = y0 + r2;

for (i = 1; i < n; i++) {

if (i & 1) {

moveto(xmax, ymin + (yinc * (double)i));

lineto(xmax - (xinc * (double)i), ymax);

}

else {

moveto(xmax - (xinc * (double)i), ymax);

lineto(xmax, ymin + (yinc * (double)i));

}

}

/* Draw upper LH corner */

xmin = x0 - r2;

ymin = y0;

xmax = x0;

ymax = y0 + r2;

for (i = 1; i < n; i++) {

if (i & 1) {

moveto(xmax - (xinc * (double)i), ymax);

lineto(xmin, ymax - (yinc * (double)i));

}

else {

moveto(xmin, ymax - (yinc * (double)i));

lineto(xmax - (xinc * (double)i), ymax);

}

}

/* Draw lower LH corner */

xmin = x0 - r2;

ymin = y0 - r2;

xmax = x0;

ymax = y0;

for (i = 1; i < n; i++) {

if (i & 1) {

moveto(xmin, ymax - (yinc * (double)i));

lineto(xmin + (xinc * (double)i), ymin);

}

else {

moveto(xmin + (xinc * (double)i), ymin);

lineto(xmin, ymax - (yinc * (double)i));

}

}

}

void radialfan(const double x0, const double y0, const double r1, const double startangle, const double dx, const double dy, const double xvec, const double yvec, const int n)

{

int i;

const double delta1 = (M_PI / 2.0) / (double)n;

const double delta2 = 1.0 / (double)n;

for (i = 0; i < n; i++) {

const double theta = startangle + ((double)i * delta1);

const double x1 = x0 + (r1 * cos(theta));

const double y1 = y0 + (r1 * sin(theta));

const double x2 = x0 + dx + (xvec * ((double)i * delta2));

const double y2 = y0 + dy + (yvec * ((double)i * delta2));

if (i & 1) {

moveto(x1, y1);

lineto(x2, y2);

}

else {

moveto(x2, y2);

lineto(x1, y1);

}

}

}

/* bordergrid --- draw grid around page with optional diagonals in corners */

void bordergrid(const double x0, const double y0, const double w, const double h, const int nx, const int ny, const int diag)

{

const double xmax = x0 + (nx * w);

const double ymax = y0 + (ny * h);

int i;

/* Draw four lines in the form of #, in optimum sequence */

moveto(x0 + w, y0);

lineto(x0 + w, ymax);

moveto(x0, ymax - h);

lineto(xmax, ymax - h);

moveto(xmax - w, ymax);

lineto(xmax - w, y0);

moveto(xmax, y0 + h);

lineto(x0, y0 + h);

/* Draw remainder of grid, along the bottom */

for (i = 2; i < (nx - 1); i++) {

moveto(x0 + (i * w), y0);

lineto(x0 + (i * w), y0 + h);

}

/* Draw lower RH corner diagonal */

if (diag) {

moveto(xmax, y0);

lineto(xmax - w, y0 + h);

}

/* Up the right-hand side */

for (i = 2; i < (ny - 1); i++) {

moveto(xmax, y0 + (i * h));

lineto(xmax - w, y0 + (i * h));

}

/* Draw upper RH corner diagonal */

if (diag) {

moveto(xmax, ymax);

lineto(xmax - w, ymax - h);

}

/* Along the top, leftwards */

for (i = nx - 2; i > 1; i--) {

moveto(x0 + (i * w), ymax);

lineto(x0 + (i * w), ymax - h);

}

/* Draw upper LH corner diagonal */

if (diag) {

moveto(x0, ymax);

lineto(x0 + w, ymax - h);

}

/* Left-hand side, downwards */

for (i = ny - 2; i > 1; i--) {

moveto(x0, y0 + (i * h));

lineto(x0 + w, y0 + (i * h));

}

/* Draw lower LH corner diagonal */

if (diag) {

moveto(x0, y0);

lineto(x0 + w, y0 + h);

}

}

/* superellipse --- draw a superellipse */

void superellipse(const double x0, const double y0, const double a, const double b, const double theta, const double d)

{

const int npts = 72;

const double delta = (2.0 * M_PI) / (double)npts;

const double sintheta = sin(theta);

const double costheta = cos(theta);

int i;

for (i = 0; i <= npts; i++) {

double sinpt, cospt;

const double t = (double)i * delta;

const double st = sin(t);

const double ct = cos(t);

if (st < 0.0)

sinpt = -pow(-st, 2.0 / d);

else

sinpt = pow(st, 2.0 / d);

if (ct < 0.0)

cospt = -pow(-ct, 2.0 / d);

else

cospt = pow(ct, 2.0 / d);

const double x = (a * cospt * costheta) - (b * sinpt * sintheta);

const double y = (a * cospt * sintheta) + (b * sinpt * costheta);

if (i == 0)

moveto(x0 + x, y0 + y);

else

lineto(x0 + x, y0 + y);

}

}

/* spiral --- draw a spiral given angle and number of turns */

void spiral(const double cx, const double cy, const double radius, const double ang, const int n)

{

double delta = 2.0 * M_PI / 72.0;

int i;

int npts = 72 * n;

moveto(cx, cy);

for (i = 1; i <= npts; i++) {

const double theta = ang + (delta * (double)i);

const double r = (radius * i) / (double)npts;

const double x = (r * cos(theta)) + cx;

const double y = (r * sin(theta)) + cy;

lineto(x, y);

}

}

/* lotus --- draw a central disc in the form of a lotus pattern */

double lotus(const double x0, const double y0, double radius, const int nleaves, const int nlayers)

{