Skip to content

Graph Labels - Charles Mallalieu 100883136 #8

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
CharlesMallalieu wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Graph Labels - Charles Mallalieu 100883136 #8

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions CondCode/CurrentCode.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,8 +16,8 @@
C.mu_0 = 1.2566370614e-6; % vacuum permeability
C.c = 299792458; % speed of light

nx = 100;
ny = 50;
nx = 100; %nx and ny are the bounds for the graphs made
ny = 50;

Acond = 1;
Bcond = 10;
Expand Down
192 changes: 192 additions & 0 deletions CondCode/GetCurrents.asv
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,192 @@
function [ Curr ] = GetCurrents(ngeo, Max, nx, ny, Acond, Bcond, plot, SimType)
global im fig fc map

if SimType == 'c'
pgeo = rand(ngeo, 2);
pgeo(:, 1) = pgeo(:, 1) * nx;
pgeo(:, 2) = pgeo(:, 2) * ny;
rgeo = rand(ngeo,1)*Max;
rcircs2 = rgeo.^2;

cMap = zeros(nx,ny);

for i = 1:nx
for j = 1:ny
cMap(i,j) = Acond;
for p = 1:ngeo
dx = (pgeo(p, 1) - i);
dy = (pgeo(p, 2) - j);
if dx * dx + dy * dy < rcircs2(p)
cMap(i, j) = Bcond;
end
end
end
end
else %if SimType == 'e'
pgeo = rand(ngeo, 3);
pgeo(:, 1) = pgeo(:, 1) * nx;
pgeo(:, 2) = pgeo(:, 2) * ny;
pgeo(:, 3) = pgeo(:, 2) * 180;
rgeo = ones(ngeo, 1) * Max;

cMap = zeros(nx, ny);

for i = 1:nx
for j = 1:ny
cMap(i, j) = Acond;
for p = 1:ngeo
point = [i, j];
el = [pgeo(p, 1), pgeo(p, 2),rgeo(p), rgeo(p) / 4, pgeo(p, 3)];
if isPointInEllipse(point, el)
cMap(i, j)= Bcond;
end
end
end
end
end


G = sparse(nx*ny);
B = zeros(1,nx*ny);

for i = 1:nx
for j = 1:ny
n = j + (i - 1) * ny;

if i == 1
G(n, :) = 0;
G(n, n) = 1;
B(n) = 1;
elseif i == nx
G(n, :) = 0;
G(n, n) = 1;
elseif j == 1
nxm = j + (i - 2) * ny;
nxp = j + (i) * ny;
nyp = j + 1 + (i - 1) * ny;

rxm = (cMap(i, j) + cMap(i - 1, j)) / 2.0;
rxp = (cMap(i, j) + cMap(i + 1, j)) / 2.0;
ryp = (cMap(i, j) + cMap(i, j + 1)) / 2.0;

G(n, n) = -(rxm+rxp+ryp);
G(n, nxm) = rxm;
G(n, nxp) = rxp;
G(n, nyp) = ryp;

elseif j == ny
nxm = j + (i - 2) * ny;
nxp = j + (i) * ny;
nym = j - 1 + (i - 1) * ny;

rxm = (cMap(i, j) + cMap(i - 1, j)) / 2.0;
rxp = (cMap(i, j) + cMap(i + 1, j)) / 2.0;
rym = (cMap(i, j) + cMap(i, j - 1)) / 2.0;

G(n, n) = -(rxm + rxp + rym);
G(n, nxm) = rxm;
G(n, nxp) = rxp;
G(n, nym) = rym;
else
nxm = j + (i-2)*ny;
nxp = j + (i)*ny;
nym = j-1 + (i-1)*ny;
nyp = j+1 + (i-1)*ny;

rxm = (cMap(i,j) + cMap(i-1,j))/2.0;
rxp = (cMap(i,j) + cMap(i+1,j))/2.0;
rym = (cMap(i,j) + cMap(i,j-1))/2.0;
ryp = (cMap(i,j) + cMap(i,j+1))/2.0;

G(n,n) = -(rxm+rxp+rym+ryp);
G(n,nxm) = rxm;
G(n,nxp) = rxp;
G(n,nym) = rym;
G(n,nyp) = ryp;
end

end
end

V = G\B';

Vmap = zeros(nx,ny);
for i = 1:nx
for j = 1:ny
n = j + (i - 1) * ny;

Vmap(i, j) = V(n);
end
end

for i = 1:nx
for j = 1:ny
if i == 1
Ex(i, j) = (Vmap(i + 1, j) - Vmap(i, j));
elseif i == nx
Ex(i, j) = (Vmap(i, j) - Vmap(i - 1, j));
else
Ex(i, j) = (Vmap(i + 1, j) - Vmap(i - 1, j)) * 0.5;
end
if j == 1
Ey(i, j) = (Vmap(i, j + 1) - Vmap(i, j));
elseif j == ny
Ey(i, j) = (Vmap(i, j) - Vmap(i, j - 1));
else
Ey(i, j) = (Vmap(i, j + 1) - Vmap(i, j - 1)) * 0.5;
end
end
end

Ex = -Ex;
Ey = -Ey;

eFlowx = cMap .* Ex;
eFlowy = cMap .* Ey;

if plot
if fc == 1
fig = figure();
end

subplot(2, 2, 1), H = surf(cMap');
set(H, 'linestyle', 'none');
view(0, 90)
title('Conductivity');
xlabel('');
ylabel('');

subplot(2, 2, 2), H = surf(Vmap');
set(H, 'linestyle', 'none');
view(0, 90)
title('Potential');
xlabel('');
ylabel('');

subplot(2, 2, 3), quiver(Ex', Ey');
axis([0 nx 0 ny]);
title('Electric Field');
xlabel('');
ylabel('');

subplot(2, 2, 4), quiver(eFlowx', eFlowy');
axis([0 nx 0 ny]);
title('Electron Flow');
xlabel('');
ylabel('');

if fc == 1
f = getframe(fig);
[im, map] = rgb2ind(f.cdata, 256, 'nodither');
im(1, 1, 1, 2) = 0;
else
f = getframe(fig);
im(:, :, 1, fc) = rgb2ind(f.cdata, map, 'nodither');
end
end

C0 = sum(eFlowx(1, :));
Cnx = sum(eFlowx(nx, :));
Curr = (C0 + Cnx) * 0.5;

end
17 changes: 17 additions & 0 deletions CondCode/GetCurrents.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -152,13 +152,30 @@
subplot(2, 2, 1), H = surf(cMap');
set(H, 'linestyle', 'none');
view(0, 90)
title('Conductivity');
xlabel('X Position (m)'); %The graph plots the index position
ylabel('Y Position (m)');
zlabel('Conductivity (S)')

subplot(2, 2, 2), H = surf(Vmap');
set(H, 'linestyle', 'none');
view(0, 90)
title('Potential');
xlabel('X Position (m)'); %The graph plots the index position
ylabel('Y Position (m)');
zlabel('Voltage (V)')

subplot(2, 2, 3), quiver(Ex', Ey');
axis([0 nx 0 ny]);
title('Electric Field');
xlabel('X Position (m)'); %The graph plots the index position
ylabel('Y Position (m)');

subplot(2, 2, 4), quiver(eFlowx', eFlowy');
axis([0 nx 0 ny]);
title('Electron Flow');
xlabel('X Position (m)'); %The graph plots the index position
ylabel('Y Position (m)');

if fc == 1
f = getframe(fig);
Expand Down
2 changes: 1 addition & 1 deletion MCTraj/MCtraj.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@
C.eps_0 = 8.854187817e-12; % vacuum permittivity
C.mu_0 = 1.2566370614e-6; % vacuum permeability
C.c = 299792458; % speed of light
C.g = 9.80665 %metres (32.1740 ft) per s�
C.g = 9.80665; %metres (32.1740 ft) per s�

nTime = 100;
nTraj = 10;
Expand Down