12 using namespace cimg_library;
20 #define ITERATIONS 2000
25 //#define debug printf
26 #define debug //printf
29 float U = -17.0; // Primary energy
30 float Va = 380.0; // Accel.
31 float Vf = -20.0;//8.18; // Focus
32 float Vd = 0.0;//0.605; // Deflection
33 float Vw = 10.0; // Wenhault
35 float grid[2][WIDTH][HEIGHT];
36 bool fixed[WIDTH][HEIGHT];
38 void SimulateElectron(float x, float y, float theta, float vMax);
41 int main(int argc, char ** argv)
43 float beta2 = pow(DELTA_X / DELTA_Y, 2);
44 //printf("Beta2 is %f\n", beta2);
47 for (int i = 0; i < WIDTH; ++i)
49 for (int j = 0; j < HEIGHT; ++j)
52 fixed[i][j] = (i <= 1 || j <= 1 || i >= WIDTH-2 || j >= HEIGHT-2);
56 debug("1. fixed[0][0] is %d\n", (int)fixed[0][0]);
60 CImg<unsigned char> src(argv[1]);
63 for (int i = 0; i < w; ++i)
65 for (int j = 0; j < h; ++j)
68 int ii = i * (WIDTH / w);
69 int jj = j * (HEIGHT / h);
70 //printf("At %d %d vs %d %d\n", i, j, ii, jj);
73 for (int d = 0; d < 3; ++d)
75 //printf("rgb[%d] is %d\n", d, rgb[d]);
76 rgb[d] = (int)(src(i, j, 0, d));
81 // Determine the electrode based on colour
82 if (rgb[0] == 255 && rgb[1] == 0 && rgb[2] == 0)
83 grid[0][ii][jj] = Va; //red
84 else if (rgb[0] == 0 && rgb[1] == 0 && rgb[2] == 255)
85 grid[0][ii][jj] = Vf; // blue
86 else if (rgb[0] == 0 && rgb[1] == 255 && rgb[2] == 0)
87 grid[0][ii][jj] = Vw; //green
88 else if (rgb[0] == 255 && rgb[1] == 255 && rgb[2] == 0)
89 grid[0][ii][jj] = Va - Vd; //yellow
90 else if (rgb[0] == 0 && rgb[1] == 255 && rgb[2] == 255)
91 grid[0][ii][jj] = Va + Vd; //light blue
92 else if (rgb[0] == 255 && rgb[1] == 0 && rgb[2] == 255)
93 grid[0][ii][jj] = U; // purple
95 grid[0][ii][jj] = 0.0;
98 // printf("Fixed at %d %d\n", ii, jj);
100 if (fixed[ii][jj] == false)
102 if (rgb[0] == 255 && rgb[1] == 255 && rgb[2] == 255)
103 fixed[ii][jj] = false;
105 fixed[ii][jj] = true;
110 debug("2. fixed[0][0] is %d\n", (int)fixed[0][0]);
112 for (int i = 0; i < WIDTH; ++i)
114 for (int j = 0; j < HEIGHT; ++j)
116 grid[1][i][j] = grid[0][i][j];
122 int n = 1; int p = 0;
124 float maxChange = 0.0;
125 debug("3. fixed[0][0] is %d\n", (int)fixed[0][0]);
127 for (k = 0; k < ITERATIONS; ++k)
131 //#pragma omp parallel for
132 for (int i = 0; i < WIDTH; ++i)
134 //#pragma omp parallel for
135 for (int j = 0; j < HEIGHT; ++j)
137 //printf("4. fixed[0][0] is %d\n", (int)fixed[0][0]);
139 if (fixed[i][j] == true)
142 grid[n][i][j] = grid[p][i][j];
148 if (i <= 1 || j <= 1 || i >= WIDTH-2 || j >= HEIGHT-2)
150 debug("WTF at %d %d\n", i, j);
157 grid[n][i][j] = (1.0 - OMEGA) * grid[p][i][j];
158 if (isnan(grid[n][i][j]))
160 debug("1. NaN at %d %d (n = %d), fixed is %d\n", i, j, n, (int)fixed[i][j]);
164 //printf("4. fixed[0][0] is %d\n", (int)fixed[0][0]);
165 grid[n][i][j] += OMEGA * (grid[n][i-1][j] + grid[p][i+1][j] + beta2 * (grid[n][i][j-1] + grid[p][i][j+1]));
166 if (isnan(grid[n][i][j]))
168 debug("5. fixed[0][0] is %d\n", (int)fixed[0][0]);
169 debug("2. NaN at %d %d (n = %d), fixed is %d\n", i, j, n, (int)fixed[i][j]);
170 debug("6. fixed[0][0] is %d\n", (int)fixed[0][0]);
174 grid[n][i][j] = grid[n][i][j] / (2.0 * (1.0 + beta2));
175 if (isnan(grid[n][i][j]))
177 debug("3. NaN at %d %d (n = %d), fixed is %d\n", i, j, n, (int)fixed[i][j]);
181 float change = (grid[n][i][j] > grid[p][i][j]) ? grid[p][i][j] - grid[n][i][j] : grid[p][i][j] - grid[n][i][j];
182 if (grid[p][i][j] != 0.00)
183 change = change / grid[p][i][j];
186 if (change >= maxChange)
193 //fprintf(stdout, "# Complete iteration %d; max change was %f\n", k, maxChange);
197 n = (n == 1) ? 0 : 1;
198 p = (p == 1) ? 0 : 1;
200 if (maxChange < EPSILON)
205 sprintf(file, "%s.dat", "main");
206 FILE * out = fopen(file, "w");
208 fprintf(out, "# Grid dimensions %d %d\n", WIDTH, HEIGHT);
209 fprintf(out, "# dx = %f, dy = %f, beta2 = %f\n", DELTA_X, DELTA_Y, beta2);
210 fprintf(out, "# Ran for %d iterations\n", k);
211 fprintf(out, "# Largest change is %f\n", maxChange);
212 for (int i = 0; i < WIDTH; ++i)
214 for (int j = 0; j < HEIGHT; ++j)
216 fprintf(out, "%f\t%f\t%f\t", i*DELTA_X, j*DELTA_Y, grid[n][i][j]);
224 fprintf(out, "# DONE\n");
225 if (out != stdout && out != stderr) fclose(out);
230 //#pragma omp parallel for
233 for (float v = 1.0; v <= 1.0; v += 0.01)
237 SimulateElectron(ox, oy, -M_PI/4.0, v);
238 SimulateElectron(ox, oy, +M_PI/4.0, v);
239 SimulateElectron(ox, oy, 0.0, v);
240 SimulateElectron(ox, oy, -M_PI/2.0, v);
241 SimulateElectron(ox, oy, +M_PI/2.0, v);
242 SimulateElectron(ox, oy, -M_PI/8.0, v);
243 SimulateElectron(ox, oy, +M_PI/8.0, v);
244 SimulateElectron(ox, oy, -M_PI/16.0, v);
245 SimulateElectron(ox, oy, +M_PI/16.0, v);
246 SimulateElectron(ox, oy, -M_PI/32.0, v);
247 SimulateElectron(ox, oy, +M_PI/32.0, v);
248 SimulateElectron(ox, oy, -M_PI/64.0, v);
249 SimulateElectron(ox, oy, +M_PI/64.0, v);
250 SimulateElectron(ox, oy, -M_PI/128.0,v);
251 SimulateElectron(ox, oy, +M_PI/128.0, v);
259 void SimulateElectron(float x, float y, float theta, float vMax)
262 sprintf(buffer, "electron%f_%f.dat", theta, vMax);
263 FILE * out = fopen(buffer, "w");
271 p[0] = x; p[1] = y; prev[0] = x; prev[1] = y;
272 v[0] = vMax * cos(theta);
273 v[1] = vMax * sin(theta);
274 for (float t = 0.0; t < 2000.0; t += dt)
277 if ((int)p[0] <= 0 || (int)p[0] >= WIDTH-1 || (int)p[1] <= 0 || (int)p[1] >= HEIGHT-1)
279 if (fixed[(int)p[0]][(int)p[1]] == true && grid[0][(int)p[0]][(int)p[1]] != U)
282 fprintf(out, "%f\t%f\t%f\t%f\t%f\n", t, p[0], p[1], v[0], v[1]);
284 float ddx = 0.5 * (grid[0][(int)(p[0])+1][(int)(p[1])] - grid[0][(int)(p[0])-1][(int)p[1]]);
285 float ddy = 0.5 * (grid[0][(int)(p[0])][(int)(p[1])+1] - grid[0][(int)(p[0])][(int)(p[1])-1]);
288 prev[0] = p[0]; prev[1] = p[1];
289 p[0] += v[0] * dt; p[1] += v[1] * dt;
292 if (out != stdout && out != stderr) fclose(out);