-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathcircle.m
More file actions
93 lines (81 loc) · 2.19 KB
/
circle.m
File metadata and controls
93 lines (81 loc) · 2.19 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
clc;clear all;
%save different electrical field
row_num = 320;
colum_num = 400;
for z = 52:100
% background + edges
for i = 1:row_num
for j = 1:colum_num
a(i,j,1) = 0;
a(i,j,2) = 0;
end
end
%sphere
center = [160 200];
radius = 50;
for i = 1:row_num
for j = 1:colum_num
if norm(center-[i j])<radius
a(i,j,2) = 100-z;
a(i,j,1) = 1;
end
end
end
%%{
%iterative finite difference method
complete = 1; % 1 all correct 0 have false
sprintf('Start')
for k=1:500000
complete = 1;
for i=2:row_num-1
for j=2:colum_num-1
if a(i,j,1)==0
if complete== 1
if abs((a(i,j,2)-differential_value(a,i,j))/a(i,j,2)) > 0.0001
complete = 0;
end
end
a(i,j,2) = differential_value(a,i,j);
end
end
end
if complete == 1
break
end
end
sprintf('End')
disp(k)
%}
%caculate electric field
divident = 1;
e_x = zeros(row_num/divident,colum_num/divident); %
e_y = zeros(row_num/divident,colum_num/divident); %
for i=2:row_num/divident-1
for j = 2:colum_num/divident-1
[e_x(i,j), e_y(i,j)] = electric_vector(a,i*divident,j*divident);
end
end
%{
%quiver for corridinate
for i = 1:row_num/divident
for j = 1:colum_num/divident
x_cordinate(i,j) = j*divident;
end
end
for i = 1:row_num/divident
for j = 1:colum_num/divident
y_cordinate(i,j) = i*divident;
end
end
quiver(x_cordinate,y_cordinate,e_x, e_y,'linewidth',2,'MaxHeadSize',2000);
%}
save(strcat('potential/save_potential_ex',int2str(z),'.mat'),'e_x')
save(strcat('potential/save_potential_ey',int2str(z),'.mat'),'e_y')
end
function y = differential_value(a, i, j)
y = (a(i+1,j,2) + a(i-1,j,2) + a(i,j+1,2) +a(i, j-1,2))/4;
end
function [x,y] = electric_vector(a,i,j)
x= -(a(i,j+1,2)-a(i,j-1,2))/(2);
y = -(a(i+1,j,2)-a(i-1,j,2))/(2);
end