clear;
%设模拟电荷点坐标
point_1=[0.007,18];point_2=[0,18.007];point_3=[-0.007,18];point_4=[0,17.993];
S_point=[point_1;point_2;point_3;point_4];
%x=zeros(1,4);y=zeros(1,4);
x=S_point(1:4,1)';
y=S_point(1:4,2)';
%设镜像模拟电荷点坐标
Ix=x;Iy=-y;
%设匹配点坐标
Mx=[0.01,0,-0.01,0];
My=[18,18.01,18,17,99];
%设镜像匹配点坐标IMx=Mx;IMy=-My;
%电压等级
U1=500000*1.05/sqrt(3);ep=8.85e-12;
%开始计算
num=4;%总的点电荷数
P=zeros(num,num);%原始系数矩阵
R=zeros(1,num);
for i=1:num
for j=1:num
R(j)=(Mx(j)-x(i))^2+(My(j)-y(i))^2;
P(j,i)=P(j,i)+(1/((4*pi*ep)*((R(j))^(1/2))));
end
end
%镜像系数矩阵
I=zeros(num,num);%原始系数矩阵
for i=1:num
for j=1:num
R(j)=(Mx(j)-Ix(i))^2+(My(j)-Iy(i))^2;
I(j,i)=I(j,i)+(1/((4*pi*ep)*((R(j))^(1/2))));
end
end
A=P+I;%合成总系数矩阵
U_P=zeros(4,1);
U2=repmat(U1,num,1);
[UHm,SHm,VHm]=svd(A);
Tau=VHm*pinv(SHm)*UHm'*U2;%模拟电荷数
%开始计算
Fx=-40:1:40;Fy=1.5;
EX=zeros(length(Fx),1);EY=zeros(length(Fx),1);
R=zeros(length(Fx),1);
%原始网格点场值
%模拟电荷在场点产生的电场
for i=1:num
B1=Tau(i);
for k=1:length(Fx)
R(k)=(Fx(k)-x(i))^2+(Fy-y(i))^2;
EX(k,1)=EX(k,1)+B1.*(Fx(k)-x(i))./((4*pi*ep).*(R(k).^(3/2)));
EY(k,1)=EX(k,1)+B1.*(Fy-y(i))./((4*pi*ep).*(R(k).^(3/2)));
end
end
%镜像电荷对计算场点的作业
for i=1:num
B1=-Tau(i);
for k=1:length(Fx)
R(k)=(Fx(k)-Ix(i))^2+(Fy-Iy(i))^2;
EX(k,1)=EX(k,1)+B1.*(Fx(k)-Ix(i))./((4*pi*ep).*(R(k).^(3/2)));
EY(k,1)=EX(k,1)+B1.*(Fy-Iy(i))./((4*pi*ep).*(R(k).^(3/2)));
end
end
%计算总合成电场及画图
E=(abs(EX).^2+abs(EY).^2).^0.5;
plot(Fx,E,'*')
hold on
plot(Fx,EX,'r.')
hold on
plot(Fx,EY,'gs')