research/simulations/egun/main.cpp~

   1 #include <stdlib.h>
   2 #include <stdio.h>
   3 #include <math.h>
   4 #include <string.h>
   5 #include <stdbool.h>
   6
   7 #include <CImg.h>
   8
   9 #include <omp.h>
  10
  11 using namespace std;
  12 using namespace cimg_library;
  13
  14 #define DELTA_X 1.0
  15 #define DELTA_Y 1.0
  16
  17 #define WIDTH 680
  18 #define HEIGHT 377
  19
  20 #define ITERATIONS 2000
  21 #define EPSILON 0.00
  22
  23 #define OMEGA 1.0
  24
  25 //#define debug printf
  26 #define debug //printf
  27
  28
  29 float U = -17.0; // Primary energy
  30 float Va = 100.0; // Accel.
  31 float Vf = -50.0;//8.18; // Focus
  32 float Vd = 0.0;//0.605; // Deflection
  33 float Vw = 10.0; // Wenhault
  34
  35 float grid[2][WIDTH][HEIGHT];
  36 bool fixed[WIDTH][HEIGHT];
  37
  38 void SimulateElectron(float x, float y, float theta, float vMax);
  39
  40
  41 int main(int argc, char ** argv)
  42 {
  43         float beta2 = pow(DELTA_X / DELTA_Y, 2);
  44         //printf("Beta2 is %f\n", beta2);
  45         omp_set_nested(1);
  46
  47         for (int i = 0; i < WIDTH; ++i)
  48         {
  49                 for (int j = 0; j < HEIGHT; ++j)
  50                 {
  51                         grid[0][i][j] = 0;
  52                         fixed[i][j] = (i <= 1 || j <= 1 || i >= WIDTH-2 || j >= HEIGHT-2);
  53                 }
  54         }
  55
  56         debug("1. fixed[0][0] is %d\n", (int)fixed[0][0]);
  57
  58
  59         // Add the electrodes
  60         CImg<unsigned char> src(argv[1]);
  61         int w = src.width();
  62         int h = src.height();
  63         for (int i = 0; i < w; ++i)
  64         {
  65                 for (int j = 0; j < h; ++j)
  66                 {
  67
  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);
  71
  72                         int rgb[3];
  73                         for (int d = 0; d < 3; ++d)
  74                         {
  75                                 //printf("rgb[%d] is %d\n", d, rgb[d]);
  76                                 rgb[d] = (int)(src(i, j, 0, d));
  77                         }
  78
  79
  80
  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
  94                         else
  95                                 grid[0][ii][jj] = 0.0;
  96
  97                         //if (fixed[ii][jj])
  98                         //      printf("Fixed at %d %d\n", ii, jj);
  99
 100                         if (fixed[ii][jj] == false)
 101                         {
 102                                 if (rgb[0] == 255 && rgb[1] == 255 && rgb[2] == 255)
 103                                         fixed[ii][jj] = false;
 104                                 else
 105                                         fixed[ii][jj] = true;
 106                         }
 107
 108                 }
 109         }
 110         debug("2. fixed[0][0] is %d\n", (int)fixed[0][0]);
 111
 112         for (int i = 0; i < WIDTH; ++i)
 113         {
 114                 for (int j = 0; j < HEIGHT; ++j)
 115                 {
 116                         grid[1][i][j] = grid[0][i][j];
 117                 }
 118         }
 119
 120
 121
 122         int n = 1; int p = 0;
 123         int k = 0;
 124         float maxChange = 0.0;
 125         debug("3. fixed[0][0] is %d\n", (int)fixed[0][0]);
 126
 127         for (k = 0; k < ITERATIONS; ++k)
 128         {
 129                 maxChange = 0.0;
 130
 131                 //#pragma omp parallel for
 132                 for (int i = 0; i < WIDTH; ++i)
 133                 {
 134                         //#pragma omp parallel for
 135                         for (int j = 0; j < HEIGHT; ++j)
 136                         {
 137                                 //printf("4. fixed[0][0] is %d\n", (int)fixed[0][0]);
 138
 139                                 if (fixed[i][j] == true)
 140                                 {
 141
 142                                         grid[n][i][j] = grid[p][i][j];
 143
 144                                         continue;
 145                                 }
 146                                 else
 147                                 {
 148                                         if (i <= 1 || j <= 1 || i >= WIDTH-2 || j >= HEIGHT-2)
 149                                         {
 150                                                 debug("WTF at %d %d\n", i, j);
 151                                                 fixed[i][j] = true;
 152                                                 continue;
 153                                         }
 154
 155                                 }
 156
 157                                 grid[n][i][j] = (1.0 - OMEGA) * grid[p][i][j];
 158                                 if (isnan(grid[n][i][j]))
 159                                 {
 160                                         debug("1. NaN at %d %d (n = %d), fixed is %d\n", i, j, n, (int)fixed[i][j]);
 161
 162                                         exit(EXIT_FAILURE);
 163                                 }
 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]))
 167                                 {
 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]);
 171                                         exit(EXIT_FAILURE);
 172                                 }
 173
 174                                 grid[n][i][j] = grid[n][i][j] / (2.0 * (1.0 + beta2));
 175                                 if (isnan(grid[n][i][j]))
 176                                 {
 177                                         debug("3. NaN at %d %d (n = %d), fixed is %d\n", i, j, n, (int)fixed[i][j]);
 178                                         exit(EXIT_FAILURE);
 179                                 }
 180
 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];
 184                                 else
 185                                         change = 1.0;
 186                                 if (change >= maxChange)
 187                                         maxChange = change;
 188
 189
 190                         }
 191                 }
 192
 193                 //fprintf(stdout, "# Complete iteration %d; max change was %f\n", k, maxChange);
 194
 195
 196
 197                 n = (n == 1) ? 0 : 1;
 198                 p = (p == 1) ? 0 : 1;
 199
 200                 if (maxChange < EPSILON)
 201                         break;
 202
 203         }
 204                 char file[100];
 205                 sprintf(file, "%s.dat", "main");
 206                 FILE * out = fopen(file, "w");
 207
 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)
 213                 {
 214                         for (int j = 0; j < HEIGHT; ++j)
 215                         {
 216                                 fprintf(out, "%f\t%f\t%f\t", i*DELTA_X, j*DELTA_Y, grid[n][i][j]);
 217                                 if (fixed[i][j])
 218                                         fprintf(out,"1\n");
 219                                 else
 220                                         fprintf(out,"0\n");
 221                         }
 222                         fprintf(out, "\n");
 223                 }
 224                 fprintf(out, "# DONE\n");
 225                 if (out != stdout && out != stderr) fclose(out);
 226
 227
 228                 //return 0;
 229
 230                 //#pragma omp parallel for
 231                 int ox = 90;
 232                 int oy = 190;
 233                 for (float v = 1.0; v <= 1.0; v += 0.01)
 234                 {
 235
 236
 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);
 252
 253                 }
 254
 255
 256
 257 }
 258
 259 void SimulateElectron(float x, float y, float theta, float vMax)
 260 {
 261         char buffer[256];
 262         sprintf(buffer, "electron%f_%f.dat", theta, vMax);
 263         FILE * out = fopen(buffer, "w");
 264
 265         float dt = 0.01;
 266         float p[2];
 267         float prev[2];
 268
 269         float v[2];
 270         float a[2];
 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)
 275         {
 276
 277                 if ((int)p[0] <= 0 || (int)p[0] >= WIDTH-1 || (int)p[1] <= 0 || (int)p[1] >= HEIGHT-1)
 278                         break;
 279                 if (fixed[(int)p[0]][(int)p[1]] == true && grid[0][(int)p[0]][(int)p[1]] != U)
 280                         break;
 281
 282                 fprintf(out, "%f\t%f\t%f\t%f\t%f\n", t, p[0], p[1], v[0], v[1]);
 283
 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]);
 286                 v[0] += ddx;
 287                 v[1] += ddy;
 288                 prev[0] = p[0]; prev[1] = p[1];
 289                 p[0] += v[0] * dt; p[1] += v[1] * dt;
 290         }
 291
 292         if (out != stdout && out != stderr) fclose(out);
 293
 294 }
 295
 296 //EOF