potential/SelfAdaptSimulation.m

clc
clear
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ӧģ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɷ<EFBFBD>
% [1].	<EFBFBD><EFBFBD>ΡΡ, <EFBFBD><EFBFBD>.A.<EFBFBD><EFBFBD>.A., һ<EFBFBD>ֽ<EFBFBD>׼ȷ<EFBFBD>ķ<EFBFBD><EFBFBD>ѵ<EFBFBD><EFBFBD>߱<EFBFBD><EFBFBD>泡ǿ<EFBFBD><EFBFBD><EFBFBD>㷽<EFBFBD><EFBFBD>. <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>, 2006(04): <EFBFBD><EFBFBD>92-96ҳ.
% [2].	<EFBFBD><EFBFBD>ϰ<EFBFBD><EFBFBD>, <EFBFBD>ظ<EFBFBD>ѹֱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>·<EFBFBD><EFBFBD><EFBFBD>Ż<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>о<EFBFBD>, 2012, <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ͨ<EFBFBD><EFBFBD>ѧ.
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>׵Ľ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>	<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>, <EFBFBD><EFBFBD>.A.<EFBFBD><EFBFBD>.A., Analysis of Conductors<EFBFBD><EFBFBD> Surface Electric Field of UHVDC Transmission Lines Based on Optimized Charge Simulation Method. <EFBFBD>ߵ<EFBFBD>ѹ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>, 2008(12): p. 2547-2551.
%%
%<EFBFBD><EFBFBD><EFBFBD>ü<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
semi_lineDistance=450;%<EFBFBD><EFBFBD><EFBFBD>Ѽ<EFBFBD><EFBFBD><EFBFBD>
semi_lineCount=6;%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
ConductorX=[-11000,11000];%<EFBFBD><EFBFBD><EFBFBD>߼<EFBFBD><EFBFBD><EFBFBD>
ConductorY=[22000,22000,];%<EFBFBD><EFBFBD><EFBFBD>߾<EFBFBD><EFBFBD>ظ߶<EFBFBD>
GroundX=[-11000,11000];%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
GroundY=[30000,30000];
CSM_N=50;%ÿһ<EFBFBD><EFBFBD><EFBFBD>ӵ<EFBFBD><EFBFBD>ߵ<EFBFBD>ģ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
subconductorR=16.8;%<EFBFBD>ӵ<EFBFBD><EFBFBD>߰뾶
phaseN=2;%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
groundN=2;%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%%
%<EFBFBD><EFBFBD><EFBFBD>õ<EFBFBD>ѹ
% Volt_=[1100/sqrt(3);1100/sqrt(3)*exp(1j*4/3*pi);1100/sqrt(3)*exp(1j*2/3*pi);];
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
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ͷ<EFBFBD><EFBFBD>ѵ<EFBFBD><EFBFBD>߼<EFBFBD><EFBFBD>಼<EFBFBD>õ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>·<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%<EFBFBD>ü<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
arc=2*pi/semi_lineCount;
CSM_arc=2*pi/CSM_N;
%<EFBFBD>ӵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĵľ<EFBFBD><EFBFBD><EFBFBD>
R=semi_lineDistance/2/sin(arc/2);
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ģ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɵ<EFBFBD>λ<EFBFBD><EFBFBD>
r1=8;
error=10000;
step=1/10;
maxLoop=round((subconductorR-r1)/step);
for Loop=1:maxLoop;
    simulationChargePos=ones(CSM_N,1);
    simulationChargeABCPos=[];
    matchPos=[];
    for I=1:CSM_N
        simulationChargePos(I)=exp(1j*((I-1)*CSM_arc+CSM_arc/2))*r1;%<EFBFBD><EFBFBD>ʱ<EFBFBD><EFBFBD>תһ<EFBFBD><EFBFBD><EFBFBD>Ƕ<EFBFBD>
    end
    for phaseLoop=1:phaseN
        for sC=1:semi_lineCount
            simulationChargeABCPos=[simulationChargeABCPos;simulationChargePos+ConductorX(phaseLoop)+1j*ConductorY(phaseLoop)+exp(1j*((sC-1)*arc+arc/2))*R];%<EFBFBD>ƶ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ӵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
            %ͬʱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƥ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>λ<EFBFBD><EFBFBD>
            matchPos=[matchPos;simulationChargePos/r1*subconductorR+ConductorX(phaseLoop)+1j*ConductorY(phaseLoop)+exp(1j*((sC-1)*arc+arc/2))*R];
        end
    end
    %<EFBFBD><EFBFBD><EFBFBD>ߵ<EFBFBD>
    for groundLoop=1:groundN
        %simulationChargeABCPos=[simulationChargeABCPos;simulationChargePos+GroundX(phaseLoop)+1j*GroundY(phaseLoop)];%<EFBFBD>ƶ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ӵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
        simulationChargeABCPos=[simulationChargeABCPos;simulationChargePos+GroundX(groundLoop)+1j*GroundY(groundLoop)];%<EFBFBD>ƶ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ӵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
            %ͬʱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƥ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>λ<EFBFBD><EFBFBD>
            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;
    %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>λϵ<EFBFBD><EFBFBD>
%     H=diag(imag(simulationChargeABCPos));
%     r=subconductorR*eye(length(imag(simulationChargeABCPos)));%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Լ<EFBFBD><EFBFBD>ξ<EFBFBD><EFBFBD><EFBFBD>
    %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>뵼<EFBFBD>ߵľ<EFBFBD><EFBFBD><EFBFBD>
    matSimulationChargePos=repmat(simulationChargeABCPos,1,length(simulationChargeABCPos));
    matMatchPos=repmat(conj(matchPos'),length(simulationChargeABCPos),1);
    CMS2MatchPointDistance=abs(matSimulationChargePos-matMatchPos);
%     conductor2conductorDistance=abs(conductor2conductorDistance-diag(diag(conductor2conductorDistance)));
    matMirrorChargePos=conj(matSimulationChargePos);%<EFBFBD>鲿ȡ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
    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);
%     P2(isnan(P2))=0;
%     P=P1+P2;
    P=P2;
    %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
    QRI=P\Volt;
    %ѡ<EFBFBD><EFBFBD><EFBFBD>鵼<EFBFBD><EFBFBD><EFBFBD><EFBFBD>һ<EFBFBD><EFBFBD><EFBFBD>Ƕ<EFBFBD>
    vrfRelA=linspace(0,2*pi,200)';%vrf=verify
    %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ӵ<EFBFBD><EFBFBD>ߵ<EFBFBD>λ<EFBFBD><EFBFBD>
    vrfRelPos=exp(1j*vrfRelA)*subconductorR;
    %<EFBFBD>ƶ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>꣬ʹ<EFBFBD><EFBFBD>֤<EFBFBD><EFBFBD><EFBFBD>ӵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĺ<EFBFBD>ʵ<EFBFBD><EFBFBD><EFBFBD>ӵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>غϡ<EFBFBD>
    vrfPos=[];
    for phaseLoop=1:phaseN
        for sC=1:semi_lineCount
            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
    %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>һ<EFBFBD><EFBFBD><EFBFBD>ĵ<EFBFBD>λϵ<EFBFBD><EFBFBD>
    matVrfPos=repmat(vrfPos,1,length(simulationChargeABCPos));
    vrf2ConductorDistance=abs(matVrfPos-repmat(conj(simulationChargeABCPos'),length(vrfPos),1));
    vrf2MirrorDistance=abs(matVrfPos-repmat(conj(conj(simulationChargeABCPos')),length(vrfPos),1));
    Pij=1/2/pi/eslong*log(vrf2MirrorDistance./vrf2ConductorDistance);
    %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѹ
    V=Pij*QRI;
    Vvalidation=[];
    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)
        %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>֤<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
    if error<0.01
        break;
    end
    r1=r1+1*step;
end
display('Finished.');
if Loop<maxLoop
    display('Converged.');
end
display(Loop);
%<EFBFBD><EFBFBD><EFBFBD>㳡ǿ
ABCy=imag(repmat(simulationChargeABCPos,1,length(vrfPos)));
ABCx=real(repmat(simulationChargeABCPos,1,length(vrfPos)));
y=imag(conj(matVrfPos'));
x=real(conj(matVrfPos'));
ERy=sum( ( (ABCy-y)./( (ABCy-y).^2+(ABCx-x).^2 )-(ABCy+y)./( (ABCy+y).^2+(ABCx-x).^2 ) ).*repmat(real(QRI),1,length(vrfPos))./2/pi/eslong,1 );
EIy=sum( ( (ABCy-y)./( (ABCy-y).^2+(ABCx-x).^2 )-(ABCy+y)./( (ABCy+y).^2+(ABCx-x).^2 ) ).*repmat(imag(QRI),1,length(vrfPos))./2/pi/eslong,1 );  
ERx=sum( ( (ABCx-x)./( (ABCy-y).^2+(ABCx-x).^2 )-(ABCx-x)./( (ABCy+y).^2+(ABCx-x).^2 ) ).*repmat(real(QRI),1,length(vrfPos))./2/pi/eslong,1 );
EIx=sum( ( (ABCx-x)./( (ABCy-y).^2+(ABCx-x).^2 )-(ABCx-x)./( (ABCy+y).^2+(ABCx-x).^2 ) ).*repmat(imag(QRI),1,length(vrfPos))./2/pi/eslong,1 );
E2=sqrt(ERy.^2+EIy.^2+ERx.^2+EIx.^2+((ERy.^2-EIy.^2+ERx.^2-EIx.^2).^2+4*(ERy.*EIy+ERx.*EIx).^2).^.5);
E3=sqrt(ERy.^2+EIy.^2+ERx.^2+EIx.^2);
Emat=1/pi/2./eslong.*repmat(conj(QRI'),length(vrfPos),1)./(vrf2ConductorDistance.^2).*(matVrfPos-repmat(conj(simulationChargeABCPos'),length(vrfPos),1))./vrf2ConductorDistance;
E=sum(Emat,2);



max(sqrt(2)*abs(E));

scatter(real(simulationChargeABCPos(1:length(simulationChargeABCPos)/1)),imag(simulationChargeABCPos(1:length(simulationChargeABCPos)/1)),[],'r');
axis equal
hold on;
scatter(real(vrfPos),imag(vrfPos),[],'k');