1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
|
#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"
"\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"
"\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"
"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";
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;
double complex func(char f[FSIZE], double complex z) {
int i = 0;
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 < 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) {
// f: function keyword
// z: coordinate
// c: complex seed
// d: depth, iteration count
int ic = 0; // iteration count
while(ic < d){
if(cabs(z) > 2)
break;
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
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));
// 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
double dx = bounds[1] - bounds[0];
double dy = bounds[3] - bounds[2];
if (dx <= 0 || dy <= 0)
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);
}
if (cur_pixel < 300)
break;
}
// 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)
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;
}
|
