X-Git-Url: https://git.ucc.asn.au/?a=blobdiff_plain;f=thesis%2Fsimulation%2Fegun%2Fmain.cpp~;fp=thesis%2Fsimulation%2Fegun%2Fmain.cpp~;h=9ee96953d36c2b1c796b07defad86a7d26f800be;hb=66ef7590cb47f237f793079dd4689c8a2941498c;hp=0000000000000000000000000000000000000000;hpb=a0fee95dc5a9513c530afb6eeab956145d0380f3;p=matches%2Fhonours.git

diff --git a/thesis/simulation/egun/main.cpp~ b/thesis/simulation/egun/main.cpp~
new file mode 100644
index 00000000..9ee96953
--- /dev/null
+++ b/thesis/simulation/egun/main.cpp~
@@ -0,0 +1,175 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <stdbool.h>
+
+#include <CImg.h>
+
+#include <omp.h>
+
+using namespace std;
+using namespace cimg_library;
+
+#define DELTA_X 1.0
+#define DELTA_Y 1.0
+
+#define WIDTH 680
+#define HEIGHT 377
+#define ITERATIONS 2000
+#define EPSILON 0.00
+
+#define OMEGA 1.0
+
+
+float grid[2][WIDTH][HEIGHT];
+bool fixed[WIDTH][HEIGHT];
+
+
+
+int main(int argc, char ** argv)
+{
+	float beta2 = pow(DELTA_X / DELTA_Y, 2);
+	//printf("Beta2 is %f\n", beta2);
+	omp_set_nested(1);
+	
+	for (int i = 0; i < WIDTH; ++i)
+	{
+		for (int j = 0; j < HEIGHT; ++j)
+		{
+			grid[0][i][j] = 0;
+			fixed[i][j] = (i == 0 || j == 0 || i == WIDTH-1 || j == HEIGHT-1);
+		}
+	}
+
+	// Add the electrodes
+	CImg<unsigned char> src(argv[1]);
+	int w = src.width();
+	int h = src.height();
+	for (int i = 0; i < w; ++i)
+	{
+		for (int j = 0; j < h; ++j)
+		{
+
+			int ii = i * (WIDTH / w);
+			int jj = j * (HEIGHT / h);
+			//printf("At %d %d vs %d %d\n", i, j, ii, jj);
+			
+			int rgb[3];
+			for (int d = 0; d < 3; ++d)
+			{
+				//printf("rgb[%d] is %d\n", d, rgb[d]);
+				rgb[d] = (int)(src(i, j, 0, d));
+			}
+
+			grid[0][ii][jj] = float(rgb[0] - rgb[2]) / 255.0;
+			
+			if (rgb[0] == rgb[1] && rgb[0] == rgb[2] && rgb[0] == 255)
+				fixed[ii][jj] = false;
+			else
+			{
+				//printf("Fix potential %d %d at %f\n", ii, jj, grid[0][ii][jj]);
+				fixed[ii][jj] = true;
+			}
+		}
+	}
+
+	int n = 1; int p = 0;
+	int k = 0;
+	float maxChange = 0.0;
+	for (k = 0; k < ITERATIONS; ++k)
+	{
+		maxChange = 0.0;
+		
+		#pragma omp parallel for	
+		for (int i = 0; i < WIDTH; ++i)
+		{
+			//#pragma omp parallel for
+			for (int j = 0; j < WIDTH; ++j)
+			{
+				if (fixed[i][j] == true)
+				{
+					grid[n][i][j] = grid[p][i][j];
+					continue;
+				}
+				
+				grid[n][i][j] = (1.0 - OMEGA) * grid[p][i][j];
+			
+
+				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]));
+				
+				grid[n][i][j] = grid[n][i][j] / (2.0 * (1.0 + beta2));
+
+				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];
+				if (grid[p][i][j] != 0.00)
+					change = change / grid[p][i][j];
+				else
+					change = 1.0;
+				if (change >= maxChange)
+					maxChange = change;
+				
+				
+
+			}
+		}
+		
+		//fprintf(stdout, "# Complete iteration %d; max change was %f\n", k, maxChange);
+
+		
+		
+		n = (n == 1) ? 0 : 1;
+		p = (p == 1) ? 0 : 1;
+
+		if (maxChange < EPSILON)
+			break;
+
+	}
+		char file[100];
+		sprintf(file, "main.dat");
+		FILE * out = fopen(file, "w");
+
+		fprintf(out, "# Grid dimensions %d %d\n", WIDTH, HEIGHT);
+		fprintf(out, "# dx = %f, dy = %f, beta2 = %f\n", DELTA_X, DELTA_Y, beta2);
+		fprintf(out, "# Ran for %d iterations\n", k);
+		fprintf(out, "# Largest change is %f\n", maxChange);
+		for (int i = 0; i < WIDTH; ++i)
+		{
+			for (int j = 0; j < HEIGHT; ++j)
+			{
+				fprintf(out, "%f\t%f\t%f\n", i*DELTA_X, j*DELTA_Y, grid[p][i][j]);
+			}
+			fprintf(out, "\n");
+		}
+		fprintf(out, "# DONE\n");
+		if (out != stdout && out != stderr) fclose(out);
+
+
+
+		out = fopen("electron.dat", "w");
+
+		float theta = 0.0;
+		float dt = 1.0;
+		float x[2] = {80, 195};
+		float v[2];
+		float a[2];
+		v[0] = 1.0 * cos(theta);
+		v[1] = 1.0 * sin(theta);		
+		for (float t = 0.0; ; t += dt)
+		{
+			if (x[0] <= 0 || x[0] >= WIDTH-1 || x[1] <= 0 || x[1] >= HEIGHT-1)
+				break;
+
+			fprintf(out, "%f\t%f\t%f\t%f\t%f\n", t, x[0], x[1], v[0], v[1]);
+		
+			float ddx = 0.5 * (grid[p][(int)(x[0])+1][(int)(x[1])] - grid[p][(int)(x[0])-1][(int)x[1]]);
+			float ddy = 0.5 * (grid[p][(int)(x[0])][(int)(x[1])+1] - grid[p][(int)(x[0])][(int)(x[1])-1]);
+			v[0] -= ddx;
+			v[1] -= ddy;
+			x[0] += v[0]; x[1] += v[1];
+		}
+
+		if (out != stdout && out != stderr) fclose(out);
+
+}
+
+//EOF