diff options
Diffstat (limited to 'julia-c/julia.c')
| -rw-r--r-- | julia-c/julia.c | 205 |
1 files changed, 205 insertions, 0 deletions
diff --git a/julia-c/julia.c b/julia-c/julia.c new file mode 100644 index 0000000..de0c187 --- /dev/null +++ b/julia-c/julia.c @@ -0,0 +1,205 @@ +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <math.h> +#include <complex.h> +#define STB_IMAGE_IMPLEMENTATION +#include "stb/stb_image.h" +#define STB_IMAGE_WRITE_IMPLEMENTATION +#include "stb/stb_image_write.h" + +const char *HELP = "julia usage:\n" + "julia <function>[string] <complex num>[a+bj] <resolution>[intxint] " + "<depth>[integer] <bounds>[x-,x+,y-,y+]\n" + "\n" + " function: the function (in complex.h) to use, (i.e. sin(x), pow(z, 3))\n" + " Optionally, a power can be appended, like 'sin2' for sin^2(z).\n" + " See README for available functions\n" + " complex_num: a complex number in the form a + bj where j is the " + " imaginary unit.\n" + " resolution: resolution (in pixels) of the image\n" + " depth: iteration count for the julia function, warning: slow at high values >100\n" + " bounds: the bounds of the graphs as a comma separated list without spaces " + " in the form: minX,maxX,minY,maxY.\n" + " Bounds above 2 are generally not interesting.\n" + "\n" + "All arguments are required. This command prints bitmap directly to stdout. " + "You should always pipe the output somewhere like:\n" + "\n./julia pow^2 -0.0567+0.678j 300x300 64 -1.5,1.5,-1.5,1.5 > file.bmp\n" + "\n" + "With Imagemagick, you can just do:\n" + "\n" + "./julia pow^2 -0.0567+0.678j 300x300 64 -1.5,1.5,-1.5,1.5 | display\n\n"; + +const char *FUNCTION_LIST[16] = { "acos" , "acosh" , "asin" , "asinh", "atan" , "atanh" , + "cos" , "cosh" , "exp" , "log" , "sin" , "sinh" , "sqrt" , "tan" , "tanh" , "pow" }; +const int RGB_MODE = 3; // I don't care about transparency, so we go for 8 bit channel +const int FSIZE = 10; // 10 is the max function size (acosh^9999, or acosh^1.125). +const double EPSILON = 0.000025; // for accounting for rounding error + +typedef struct { + uint8_t R; + uint8_t G; + uint8_t B; +} pixel_t; // easier pixel management + +double complex func(char f[FSIZE], double complex z) { + // This allows the user to choose more complicated functions available + // in the complex.h GNU library according to a keyword f + // + // returns the value of the chosen function at point z + + int i; + double complex e; + char copy[FSIZE], fname[5];// max would be acosh\0, as in acosh^5 + char *parser; + strcpy(copy, f); + parser = strtok(copy, "^"); + + strcpy(fname, parser); + parser = strtok(NULL, "^"); // second will be exponent + + if (parser == NULL) // No ^ was found + e = 1; + else + e = (double complex) atof(parser); + + for (i = 0; i < 16; i++) + if (!strcmp(fname, FUNCTION_LIST[i])) + break; + + switch (i) { + case 0: return cpow(cacos(z), e); + case 1: return cpow(cacosh(z), e); + case 2: return cpow(casin(z), e); + case 3: return cpow(casinh(z), e); + case 4: return cpow(catan(z), e); + case 5: return cpow(catanh(z), e); + case 6: return cpow(ccos(z), e); + case 7: return cpow(ccosh(z), e); + case 8: return cpow(cexp(z), e); + case 9: return cpow(clog(z), e); + case 10: return cpow(csin(z), e); + case 11: return cpow(csinh(z), e); + case 12: return cpow(csqrt(z), e); + case 13: return cpow(ctan(z), e); + case 14: return cpow(ctanh(z), e); + case 15: return cpow(z, e); + default: return z * z; // if function can't be parsed, just z^2 + } +} + +int julia(char f[FSIZE], double complex z, double complex c, int d) { + // The core Julia set algorithm. With a function f(z), for a certain + // point z and complex seed c, evaluate f(f(z) + c until the function + // diverges OR reaches the max iteration count d + // + // f: function keyword + // z: coordinate + // c: complex seed + // d: depth, max iteration count (and changes inner color!) + // + // returns an integer representing how many iterations the function + // performed before diverging + + int ic = 0; // iteration count + + while (ic < d){ + if (cabs(z) > 2); // this is arbitrary + break; // the rationale is that the function is complex + // so in certain spots it will be cyclical + // and bounded under ~1 or ~2 + // once z gets higher than this, it diverges quickly + // remember, i, i^2, i^3, i^4 + // makes i,-1,-i,1 ad infinitum + // see rosettacode.org/wiki/Julia_set for examples + z = func(f, z) + c; + ic++; + } + return ic; +} + +int write_image(char f[FSIZE], double complex seed, int r[2], int depth, double bounds[4]) { + // Writes the image to stdout according to function parameters + // f: the function keyword + // seed: a complex number a+bj + // r[]: the XxY resolution of the image + // depth: the iteration count to use + // bounds: the grid on which to plot + // + // returns the status code (success or failure) writing the image + + double x, y; + double dx = bounds[1] - bounds[0]; + double dy = bounds[3] - bounds[2]; + int w = 0; // width, this will keep us under 300 pixels per row + // without this, the graph will be shifted if there's + // rounding error + // e.g. 300x300 resolution and -1.3 1.3 + // bounds, makes x go up by 2.6/300 which does not divide + // nicely + + // In this image library, data is stored in a + // 1-D array separated by row, so you need to fill + // the pixels horizontally first, and *backwards* if + // we want the grid right side up (i.e. -1,-1,0,1 will print just the + // top half) + int pixels = r[1] * r[0]; + int cur_pixel = pixels - 1; + pixel_t *img_data = (pixel_t *) STBI_MALLOC(pixels * sizeof(pixel_t)); + char const *filename = "/dev/stdout"; + + + if (dx <= 0 || dy <= 0) + return 0; // failure (to match stbi_write) + + y = bounds[2]; + x = bounds[0]; + while (cur_pixel > -1) { + if (x > bounds[1] || w >= r[0]) { + // reset horizontal axis + w = 0; + x = bounds[0]; + // increment vertical + y += dy/((double) r[1]); + fprintf(stderr, "x: %5f, y: %5f p: %d\n", x, y, cur_pixel); + } + double complex z = x + y * I; + int out = julia(f, z, seed, depth); + img_data[cur_pixel].R = out % 32 * 8; + img_data[cur_pixel].G = out % 16 * 16; + img_data[cur_pixel].B = out % 8 * 32; + cur_pixel -= 1; + x += dx/((double) r[0]); + w++; + } + + // returns 0 on failure + return stbi_write_bmp(filename, r[0], r[1], RGB_MODE, img_data); +} + +int main(int argc, char *argv[]) { + double complex seed; + char f[FSIZE]; + int depth; + int res[2]; + double b[4]; + if (argc != 6) { + fprintf(stderr, "%s\n", HELP); + return 1; + } + else { + double re, img; + sscanf(argv[1], "%8s", f); + sscanf(argv[2], "%10lf+%10lfj", &re, &img); + sscanf(argv[3], "%10dx%10d", &res[0], &res[1]); + sscanf(argv[4], "%4d", &depth); + sscanf(argv[5], "%10lf,%10lf,%10lf,%10lf", &b[0], &b[1], &b[2], &b[3]); + seed = re + img * I; + if (!write_image(f, seed, res, depth, b)) { + fprintf(stderr, "Error in arguments, could not make image\n"); + return 1; + } + } + return 0; +} |