等效电荷法计算电场

这是一段模拟电荷法（scm)的matlab源程序，求解了超高压输电线路下方工频电场，仿真出了离地面1.5m处的二维电场模型，该程序经过验证，准确无误，谢谢。

clear
grid on
x=-40:0.5:40;
a=13.72;b=12.19;
y=1.5;
c=8.85e-12;
QR=2*pi*c*[54705 -73523 24287];
Qi=2*pi*c*[46378 19701 -67517];
Q=QR+Qi*j;
A=[（y-b）./sqrt（（a-x）.^2+（b-1.5）.^2）-（y-b）./sqrt（（a-x）.^2+（-b-1.5）.^2）;（y-b）./sqrt（x.^2+（b-1.5）.^2）-（y-b）./sqrt（x.^2+（-b-1.5）.^2）;
（y-b）./sqrt（（-a-x）.^2+（b-1.5）.^2）-（y-b）./sqrt（（-a-x）.^2+（-b-1.5）.^2）];
Ey=abs（Q*A./2*pi*c）*1e21;
subplot（331）;
y=1.5;
plot（xy）
hold on;
plot（[0 a -a -a 0 a][b b b -b -b -b]‘o‘）
title（‘算例示意图‘）;
ylabel（‘离地面高度‘）;
xlabel（‘水平距离‘）;
subplot（332）;
plot（xEy）;
title（‘算例一的电场强度曲线‘）;
xlabel（‘x（m）‘）;
ylabel（‘E（V/m）‘）;
axis（[-40 40 0 12000]）;
a1=20;b1=8;c1=25;
QR=2*pi*c*[75864 -40753 -31074];
Qi=2*pi*c*[-7926 70575 -69662];
Q=QR+Qi*j;
A=[（y-b1）./sqrt（（a1-x）.^2+（b1-1.5）.^2）-（y-b1）./sqrt（（a1-x）.^2+（-b1-1.5）.^2）;（y-c1）./sqrt（x.^2+（b1-1.5）.^2）-（y-c1）./sqrt（x.^2+（-b1-1.5）.^2）;
（y-b1）./sqrt（（-a1-x）.^2+（b1-1.5）.^2）-（y-b1）./sqrt（（-a1-x）.^2+（-b1-1.5）.^2）];
B=[（x-a1）./sqrt（（a1-x）.^2+（b1-1.5）.^2）-（x-a1）./sqrt（（a1-x）.^2+（-b1-1.5）.^2）;（x）./sqrt（x.^2+（b1-1.5）.^2）-（x）./sqrt（x.^2+（-b1-1.5）.^2）;
（x+a1）./sqrt（（-a1-x）.^2+（b1-1.5）.^2）-（x+a1）./sqrt（（-a1-x）.^2+（-b1-1.5）.^2）];
Ey1=abs（Q*A./2*pi*c）*1e20;
Ey2=abs（Q*B./2*pi*c）*1e20;
Ey3=sqrt（Ey1.^2+Ey2.^2）;
subplot（333）;
y=1.5;
plot（xy）
hold on;
plot（[0 a1 -a1 -a1 0 a1][c1 b1 b1 -b1 -c1 -b1]‘o‘）
title（‘算例示意图‘）;
ylabel（‘离地面高度‘）;
xlabel（‘水平距离‘）;
subplot（334）;
plot（xEy1）;
title（‘算例二的电场强度曲线‘）;
xlabel（‘x（m）‘）;
ylabel（‘E（V/m）‘）;
subplot（335）;
plot（xEy2）
subplot（336）;
plot（xEy3）

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       1592  2008-05-16 21:10  Untitled8.m

 -rw-r--r--       221  2013-07-17 15:17  readme_verysource.com.txt

----------- ---------  ---------- -----  ----

                 1813                    2