thesis/simulation/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 #define ITERATIONS 2000
  20 #define EPSILON 0.00
  21
  22 #define OMEGA 1.0
  23
  24
  25 float grid[2][WIDTH][HEIGHT];
  26 bool fixed[WIDTH][HEIGHT];
  27
  28
  29
  30 int main(int argc, char ** argv)
  31 {
  32         float beta2 = pow(DELTA_X / DELTA_Y, 2);
  33         //printf("Beta2 is %f\n", beta2);
  34         omp_set_nested(1);
  35
  36         for (int i = 0; i < WIDTH; ++i)
  37         {
  38                 for (int j = 0; j < HEIGHT; ++j)
  39                 {
  40                         grid[0][i][j] = 0;
  41                         fixed[i][j] = (i == 0 || j == 0 || i == WIDTH-1 || j == HEIGHT-1);
  42                 }
  43         }
  44
  45         // Add the electrodes
  46         CImg<unsigned char> src(argv[1]);
  47         int w = src.width();
  48         int h = src.height();
  49         for (int i = 0; i < w; ++i)
  50         {
  51                 for (int j = 0; j < h; ++j)
  52                 {
  53
  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);
  57
  58                         int rgb[3];
  59                         for (int d = 0; d < 3; ++d)
  60                         {
  61                                 //printf("rgb[%d] is %d\n", d, rgb[d]);
  62                                 rgb[d] = (int)(src(i, j, 0, d));
  63                         }
  64
  65                         grid[0][ii][jj] = float(rgb[0] - rgb[2]) / 255.0;
  66
  67                         if (rgb[0] == rgb[1] && rgb[0] == rgb[2] && rgb[0] == 255)
  68                                 fixed[ii][jj] = false;
  69                         else
  70                         {
  71                                 //printf("Fix potential %d %d at %f\n", ii, jj, grid[0][ii][jj]);
  72                                 fixed[ii][jj] = true;
  73                         }
  74                 }
  75         }
  76
  77         int n = 1; int p = 0;
  78         int k = 0;
  79         float maxChange = 0.0;
  80         for (k = 0; k < ITERATIONS; ++k)
  81         {
  82                 maxChange = 0.0;
  83
  84                 #pragma omp parallel for
  85                 for (int i = 0; i < WIDTH; ++i)
  86                 {
  87                         //#pragma omp parallel for
  88                         for (int j = 0; j < WIDTH; ++j)
  89                         {
  90                                 if (fixed[i][j] == true)
  91                                 {
  92                                         grid[n][i][j] = grid[p][i][j];
  93                                         continue;
  94                                 }
  95
  96                                 grid[n][i][j] = (1.0 - OMEGA) * grid[p][i][j];
  97
  98
  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]));
 100
 101                                 grid[n][i][j] = grid[n][i][j] / (2.0 * (1.0 + beta2));
 102
 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];
 106                                 else
 107                                         change = 1.0;
 108                                 if (change >= maxChange)
 109                                         maxChange = change;
 110
 111
 112
 113                         }
 114                 }
 115
 116                 //fprintf(stdout, "# Complete iteration %d; max change was %f\n", k, maxChange);
 117
 118
 119
 120                 n = (n == 1) ? 0 : 1;
 121                 p = (p == 1) ? 0 : 1;
 122
 123                 if (maxChange < EPSILON)
 124                         break;
 125
 126         }
 127                 char file[100];
 128                 sprintf(file, "main.dat");
 129                 FILE * out = fopen(file, "w");
 130
 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)
 136                 {
 137                         for (int j = 0; j < HEIGHT; ++j)
 138                         {
 139                                 fprintf(out, "%f\t%f\t%f\n", i*DELTA_X, j*DELTA_Y, grid[p][i][j]);
 140                         }
 141                         fprintf(out, "\n");
 142                 }
 143                 fprintf(out, "# DONE\n");
 144                 if (out != stdout && out != stderr) fclose(out);
 145
 146
 147
 148                 out = fopen("electron.dat", "w");
 149
 150                 float theta = 0.0;
 151                 float dt = 1.0;
 152                 float x[2] = {80, 195};
 153                 float v[2];
 154                 float a[2];
 155                 v[0] = 1.0 * cos(theta);
 156                 v[1] = 1.0 * sin(theta);
 157                 for (float t = 0.0; ; t += dt)
 158                 {
 159                         if (x[0] <= 0 || x[0] >= WIDTH-1 || x[1] <= 0 || x[1] >= HEIGHT-1)
 160                                 break;
 161
 162                         fprintf(out, "%f\t%f\t%f\t%f\t%f\n", t, x[0], x[1], v[0], v[1]);
 163
 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]);
 166                         v[0] -= ddx;
 167                         v[1] -= ddy;
 168                         x[0] += v[0]; x[1] += v[1];
 169                 }
 170
 171                 if (out != stdout && out != stderr) fclose(out);
 172
 173 }
 174
 175 //EOF