clc
clear
%% 自适应模拟电荷法
% [1].	任巍巍, 孙.A.宗.A., 一种较准确的分裂导线表面场强计算方法. 电网技术, 2006(04): 第92-96页.
% [2].	陈习文, 特高压直流输电线路电磁环境的研究, 2012, 北京交通大学.
%%
%设置几个参数
semi_lineDistance=257;%分裂间距
semi_lineCount=4;%分裂数
ConductorX=[-14500,14500];%导线距地高度
ConductorY=[16500,16500];%导线间距
CSM_N=80;%每一个子导线的模拟电荷数
subconductorR=30;%子导线半径
%%
%设置电压
Volt=[500;500;500;500;-500;-500;-500;-500];
Volt=[500*ones(CSM_N*semi_lineCount,1);-500*ones(CSM_N*semi_lineCount,1)];
%按分裂数和分裂导线间距布置单相线路导线
%用极坐标
arc=2*pi/semi_lineCount;
CSM_arc=2*pi/CSM_N;
%子导线中心到导线中心的距离
R=semi_lineDistance/2/sin(arc/2);
%计算模拟电荷的位置
r1=20;
error=10000;
step=1/10;
maxLoop=round((subconductorR-r1)/step);
for Loop=1:maxLoop;
    simulationChargePos=ones(CSM_N,1);
    simulationChargeAPos=[];
    simulationChargeBPos=[];
    for sC=1:semi_lineCount
        for I=1:CSM_N
%             simulationChargePos(I)=exp(1j*((I-1)*CSM_arc+CSM_arc/2))*(R+r1);%逆时针转一个角度
              simulationChargePos(I)=exp(1j*((I-1)*CSM_arc+CSM_arc/2))*r1;%逆时针转一个角度
        end
        simulationChargeAPos=[simulationChargeAPos;simulationChargePos+ConductorX(1)+1j*ConductorY(1)+exp(1j*((sC-1)*arc+arc/2))*R];%移动到子导线中心
        simulationChargeBPos=[simulationChargeBPos;simulationChargePos+ConductorX(2)+1j*ConductorY(2)+exp(1j*((sC-1)*arc+arc/2))*R];%移动到子导线中心
    end
%     simulationChargeAPos=simulationChargePos+ConductorX(1)+1j*ConductorY(1);
%     simulationChargeBPos=simulationChargePos+ConductorX(2)+1j*ConductorY(2);
    simulationChargePos=[simulationChargeAPos;simulationChargeBPos];
    %计算电位系数
    H=diag(imag(simulationChargePos));
    r=subconductorR*eye(length(imag(simulationChargePos)));%导线自几何均距
    %导线与导线的距离
    matSimulationChargePos=repmat(simulationChargePos,1,length(simulationChargePos));
    conductor2conductorDistance=matSimulationChargePos-conj(matSimulationChargePos');
    conductor2conductorDistance=abs(conductor2conductorDistance-diag(diag(conductor2conductorDistance)));
    matMirrorChargePos=conj(matSimulationChargePos);%虚部取负号
    conductor2MirrorDistance=matSimulationChargePos-conj(matMirrorChargePos');
    conductor2MirrorDistance=abs(conductor2MirrorDistance-diag(diag(conductor2MirrorDistance)));
    eslong=8.854187817*10;
    P1=1/2/pi/eslong*log(2*H./r);
    P1(isnan(P1))=0;
    P2=1/2/pi/eslong*log(conductor2MirrorDistance./conductor2conductorDistance);
    P2(isnan(P2))=0;
    P=P1+P2;
    %求电荷
    QRI=P\Volt;
    %以下是验证部分
    if error<0.0001
        break;
    end
    %选检验导线上一个角度
    vrfRelA=linspace(0,2*pi,200)';%vrf=verify
    %计算检验点相对于子导线的位置
    vrfRelPos=exp(1j*vrfRelA)*subconductorR;
    %移动坐标，使验证的子导线中心和实际子导线中心重合。
    vrfPos=[];
    for sC=1:semi_lineCount
        vrfPos=[vrfPos;exp(1j*((sC-1)*arc+arc/2))*R+ConductorX(1)+1j*ConductorY(1)+vrfRelPos];
    end
%     vrfPos=ConductorX(1)+1j*ConductorY(1)+vrfRelPos;
    %计算这一点的电位系数
    matVrfPos=repmat(vrfPos,1,length(simulationChargePos));
    vrf2ConductorDistance=abs(matVrfPos-repmat(conj(simulationChargePos'),length(vrfPos),1));
    vrf2MirrorDistance=abs(matVrfPos-repmat(conj(conj(simulationChargePos')),length(vrfPos),1));
    Pij=1/2/pi/eslong*log(vrf2MirrorDistance./vrf2ConductorDistance);
    %计算电压
    V=Pij*QRI;
    error=sum(abs(V-500)./500)/length(V);
    r1=r1+step;
end
display('Finished.');
if Loop<maxLoop
    display('Converged.');
end
display(Loop);
scatter(real(simulationChargeAPos),imag(simulationChargeAPos),[],'r');
hold on;
scatter(real(vrfPos),imag(vrfPos),[],'k');
% scatter(real([simulationChargeAPos(1:10);]),imag([simulationChargeAPos(1:10);]),10,'red');
% hold;
% scatter(real([simulationChargeAPos(11:20);]),imag([simulationChargeAPos(11:20);]),10,'blue');
% scatter(real([simulationChargeAPos(21:30);]),imag([simulationChargeAPos(21:30);]),10,'green');
% scatter(real([simulationChargeAPos(31:40);]),imag([simulationChargeAPos(31:40);]),10,'yellow');