1.把模拟电荷法的地线加上了

2.初始半径太小时数值性差，接近奇异

Signed-off-by: dmy@lab <dmy@lab.com>

SelfAdaptSimulation.m

@@ -10,17 +10,23 @@ semi_lineDistance=450;%
 semi_lineCount=6;%分裂数
 ConductorX=[-11000,11000];%导线间距
 ConductorY=[22000,22000,];%导线距地高度
-CSM_N=40;%每一个子导线的模拟电荷数
+GroundX=[-11000,11000];%地线坐标
+GroundY=[30000,30000];
+CSM_N=50;%每一个子导线的模拟电荷数
 subconductorR=16.8;%子导线半径
 phaseN=2;%相数，单回三相
+groundN=2;%地线数量
 %%
 %设置电压
 % Volt_=[1100/sqrt(3);1100/sqrt(3)*exp(1j*4/3*pi);1100/sqrt(3)*exp(1j*2/3*pi);];
-Volt_=[800;-800;];
+Volt_=[800;-800;0;0];
 Volt=[];
 for vLoop=1:phaseN
     Volt=[Volt;Volt_(vLoop)*ones(CSM_N*semi_lineCount,1);];
 end
+for vLoop=1:groundN
+    Volt=[Volt;Volt_(vLoop+phaseN)*ones(CSM_N,1);];
+end
 %按分裂数和分裂导线间距布置单相线路导线
 %用极坐标
 arc=2*pi/semi_lineCount;
@@ -28,7 +34,7 @@ CSM_arc=2*pi/CSM_N;
 %子导线中心到导线中心的距离
 R=semi_lineDistance/2/sin(arc/2);
 %计算模拟电荷的位置
-r1=5;
+r1=8;
 error=10000;
 step=1/10;
 maxLoop=round((subconductorR-r1)/step);
@@ -46,6 +52,14 @@ for Loop=1:maxLoop;
             matchPos=[matchPos;simulationChargePos/r1*subconductorR+ConductorX(phaseLoop)+1j*ConductorY(phaseLoop)+exp(1j*((sC-1)*arc+arc/2))*R];
         end
     end
+    %地线的
+    for groundLoop=1:groundN
+        %simulationChargeABCPos=[simulationChargeABCPos;simulationChargePos+GroundX(phaseLoop)+1j*GroundY(phaseLoop)];%移动到子导线中心
+        simulationChargeABCPos=[simulationChargeABCPos;simulationChargePos+GroundX(groundLoop)+1j*GroundY(groundLoop)];%移动到子导线中心
+            %同时计算匹配点的位置
+            matchPos=[matchPos;simulationChargePos/r1*subconductorR+GroundX(groundLoop)+1j*GroundY(groundLoop)];
+%             matchPos=[matchPos;GroundX(groundLoop)+1j*GroundY(groundLoop)];
+    end
 %     simulationChargeAPos=simulationChargePos+ConductorX(1)+1j*ConductorY(1);
 %     simulationChargeBPos=simulationChargePos+ConductorX(2)+1j*ConductorY(2);
 %     simulationChargePos=simulationChargeABCPos;
@@ -61,6 +75,7 @@ for Loop=1:maxLoop;
     mirrorCharge2MatchPointDistance=abs(matMirrorChargePos-matMatchPos);
 %     conductor2MirrorDistance=abs(conductor2MirrorDistance-diag(diag(conductor2MirrorDistance)));
     eslong=8.854187817*10^-12*1000;
+    eslong=1;
 %     P1=1/2/pi/eslong*log(2*H./r);
 %     P1(isnan(P1))=0;
     P2=1/2/pi/eslong*log(mirrorCharge2MatchPointDistance./CMS2MatchPointDistance);
@@ -80,6 +95,9 @@ for Loop=1:maxLoop;
             vrfPos=[vrfPos;exp(1j*((sC-1)*arc+arc/2))*R+ConductorX(phaseLoop)+1j*ConductorY(phaseLoop)+vrfRelPos];
         end
     end
+    for groundLoop=1:groundN
+            vrfPos=[vrfPos;GroundX(groundLoop)+1j*GroundY(groundLoop)+vrfRelPos];
+    end
     %计算这一点的电位系数
     matVrfPos=repmat(vrfPos,1,length(simulationChargeABCPos));
     vrf2ConductorDistance=abs(matVrfPos-repmat(conj(simulationChargeABCPos'),length(vrfPos),1));
@@ -91,11 +109,15 @@ for Loop=1:maxLoop;
     for phaseLoop=1:phaseN
         Vvalidation=[Vvalidation;Volt_(phaseLoop)*ones(semi_lineCount*200,1);];
     end
+    for groundLoop=1:groundN
+        Vvalidation=[Vvalidation;Volt_(groundLoop+phaseN)*ones(200,1);];
+    end
     error=abs((V-Vvalidation)./Vvalidation);
+    error=[error(1:phaseN*semi_lineCount*200);abs(V(phaseN*semi_lineCount*200+1:end)-Vvalidation(phaseN*semi_lineCount*200+1:end))];
 %     error(isinf(error))=0;
     error=sum(error)/length(Vvalidation)
         %以下是验证部分
-    if error<0.001
+    if error<0.01
         break;
     end
     r1=r1+1*step;