12 using namespace cimg_library;
19 #define ITERATIONS 2000
25 float grid[2][WIDTH][HEIGHT];
26 bool fixed[WIDTH][HEIGHT];
30 int main(int argc, char ** argv)
32 float beta2 = pow(DELTA_X / DELTA_Y, 2);
33 //printf("Beta2 is %f\n", beta2);
36 for (int i = 0; i < WIDTH; ++i)
38 for (int j = 0; j < HEIGHT; ++j)
41 fixed[i][j] = (i == 0 || j == 0 || i == WIDTH-1 || j == HEIGHT-1);
46 CImg<unsigned char> src(argv[1]);
49 for (int i = 0; i < w; ++i)
51 for (int j = 0; j < h; ++j)
54 int ii = i * (WIDTH / w);
55 int jj = j * (HEIGHT / h);
56 //printf("At %d %d vs %d %d\n", i, j, ii, jj);
59 for (int d = 0; d < 3; ++d)
61 //printf("rgb[%d] is %d\n", d, rgb[d]);
62 rgb[d] = (int)(src(i, j, 0, d));
65 grid[0][ii][jj] = float(rgb[0] - rgb[2]) / 255.0;
67 if (rgb[0] == rgb[1] && rgb[0] == rgb[2] && rgb[0] == 255)
68 fixed[ii][jj] = false;
71 //printf("Fix potential %d %d at %f\n", ii, jj, grid[0][ii][jj]);
79 float maxChange = 0.0;
80 for (k = 0; k < ITERATIONS; ++k)
84 #pragma omp parallel for
85 for (int i = 0; i < WIDTH; ++i)
87 //#pragma omp parallel for
88 for (int j = 0; j < WIDTH; ++j)
90 if (fixed[i][j] == true)
92 grid[n][i][j] = grid[p][i][j];
96 grid[n][i][j] = (1.0 - OMEGA) * grid[p][i][j];
99 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]));
101 grid[n][i][j] = grid[n][i][j] / (2.0 * (1.0 + beta2));
103 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];
104 if (grid[p][i][j] != 0.00)
105 change = change / grid[p][i][j];
108 if (change >= maxChange)
116 //fprintf(stdout, "# Complete iteration %d; max change was %f\n", k, maxChange);
120 n = (n == 1) ? 0 : 1;
121 p = (p == 1) ? 0 : 1;
123 if (maxChange < EPSILON)
128 sprintf(file, "main.dat");
129 FILE * out = fopen(file, "w");
131 fprintf(out, "# Grid dimensions %d %d\n", WIDTH, HEIGHT);
132 fprintf(out, "# dx = %f, dy = %f, beta2 = %f\n", DELTA_X, DELTA_Y, beta2);
133 fprintf(out, "# Ran for %d iterations\n", k);
134 fprintf(out, "# Largest change is %f\n", maxChange);
135 for (int i = 0; i < WIDTH; ++i)
137 for (int j = 0; j < HEIGHT; ++j)
139 fprintf(out, "%f\t%f\t%f\n", i*DELTA_X, j*DELTA_Y, grid[p][i][j]);
143 fprintf(out, "# DONE\n");
144 if (out != stdout && out != stderr) fclose(out);
148 out = fopen("electron.dat", "w");
152 float x[2] = {80, 195};
155 v[0] = 1.0 * cos(theta);
156 v[1] = 1.0 * sin(theta);
157 for (float t = 0.0; ; t += dt)
159 if (x[0] <= 0 || x[0] >= WIDTH-1 || x[1] <= 0 || x[1] >= HEIGHT-1)
162 fprintf(out, "%f\t%f\t%f\t%f\t%f\n", t, x[0], x[1], v[0], v[1]);
164 float ddx = 0.5 * (grid[p][(int)(x[0])+1][(int)(x[1])] - grid[p][(int)(x[0])-1][(int)x[1]]);
165 float ddy = 0.5 * (grid[p][(int)(x[0])][(int)(x[1])+1] - grid[p][(int)(x[0])][(int)(x[1])-1]);
168 x[0] += v[0]; x[1] += v[1];
171 if (out != stdout && out != stderr) fclose(out);