diff options
Diffstat (limited to 'julia-c/julia.c')
| -rw-r--r-- | julia-c/julia.c | 119 |
1 files changed, 78 insertions, 41 deletions
diff --git a/julia-c/julia.c b/julia-c/julia.c index adc80a6..de0c187 100644 --- a/julia-c/julia.c +++ b/julia-c/julia.c @@ -12,23 +12,24 @@ const char *HELP = "julia usage:\n" "julia <function>[string] <complex num>[a+bj] <resolution>[intxint] " "<depth>[integer] <bounds>[x-,x+,y-,y+]\n" "\n" - "\tfunction: the function (in complex.h) to use, (i.e. sin(x), pow(z, 3))\n" - "\tOptionally, a power can be appended, like 'sin2' for sin^2(z).\n" - "\tSee README for available functions\n" - "\tcomplex_num: a complex number in the form a + bj where j is the " - "imaginary unit.\n" - "\tresolution: resolution (in pixels) of the image\n" - "\tdepth: iteration count for the julia function, warning: slow at high values >100\n" - "\tbounds: the bounds of the graphs as a comma separated list without spaces " - "in the form: minX,maxX,minY,maxY.\n" - "\tBounds above 2 are generally not interesting.\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"; + "./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" }; @@ -40,10 +41,15 @@ typedef struct { uint8_t R; uint8_t G; uint8_t B; -} pixel_t; +} pixel_t; // easier pixel management double complex func(char f[FSIZE], double complex z) { - int i = 0; + // 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; @@ -58,7 +64,7 @@ double complex func(char f[FSIZE], double complex z) { else e = (double complex) atof(parser); - for ( ; i < 16; i++) + for (i = 0; i < 16; i++) if (!strcmp(fname, FUNCTION_LIST[i])) break; @@ -84,15 +90,29 @@ double complex func(char f[FSIZE], double complex z) { } 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, iteration count + // 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) - break; + 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++; } @@ -100,54 +120,71 @@ int julia(char f[FSIZE], double complex z, double complex c, int d) { } 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 + // 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; - int pixels = r[1] * r[0]; - int cur_pixel = pixels; - char const *filename = "/dev/stdout"; - pixel_t *img_data = (pixel_t *) STBI_MALLOC(pixels * sizeof(pixel_t)); + 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 + // 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"; - double dx = bounds[1] - bounds[0]; - double dy = bounds[3] - bounds[2]; + if (dx <= 0 || dy <= 0) - return 0; // failure (to match stbi_write) + return 0; // failure (to match stbi_write) - for(y = bounds[2]; y <= bounds[3]; y += dy/((double) r[1])) { - for(x = bounds[0]; x <= bounds[1]; x += dx/((double) r[0])) { - 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; - fprintf(stderr, "x: %5f, y: %5f p: %d\n", x, y, cur_pixel); + 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); } - if (cur_pixel < 300) - break; - + 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]; // 10 is the max function size (acosh^9999) + char f[FSIZE]; int depth; int res[2]; double b[4]; - if(argc != 6) { + if (argc != 6) { fprintf(stderr, "%s\n", HELP); return 1; } |