clc
clear
close all
a=exp(1j*2*pi/3);
Tp2f=1/3*[1 1 1;
    1 a a^2;
    1 a^2 a];
sigma=0.03;
rLoadP=[10;10;10];
rLoadQ=[9;9;9];
eLoadP=normrnd(0,sigma,length(rLoadP),1);
mLoadP=rLoadP.*(1+eLoadP);
eLoadQ=normrnd(0,sigma,length(rLoadQ),1);
mLoadQ=rLoadQ.*(1+eLoadQ);
rVoltP=[1;1*exp(1j*-120/180*pi);1*exp(1j*+120/180*pi)];
mVoltP=rVoltP.*(1+normrnd(0,sigma,length(rVoltP),1));
rI=conj((rLoadP+1j*rLoadQ)./rVoltP);
mI=conj((mLoadP+1j*mLoadQ)./mVoltP);
rFI=Tp2f*rI;
mFI=Tp2f*mI;
%¼ÆË㲻ƽºâ¶È
LmP=(sum(mLoadP.^4))/(sum(mLoadP.^2)*sum(mLoadP.^2));
LmQ=(sum(mLoadQ.^4))/(sum(mLoadQ.^2)*sum(mLoadQ.^2));
ubmP=sqrt( (1-sqrt(3-6*LmP))/(1+sqrt(3-6*LmP) ) )*100;
ubmQ=sqrt( (1-sqrt(3-6*LmQ))/(1+sqrt(3-6*LmQ) ) )*100;
LrP=(sum(rLoadP.^4))/(sum(rLoadP.^2)*sum(rLoadP.^2));
LrQ=(sum(rLoadQ.^4))/(sum(rLoadQ.^2)*sum(rLoadQ.^2));
ubrP=sqrt( (1-sqrt(3-6*LrP))/(1+sqrt(3-6*LrP) ) )*100;
ubrQ=sqrt( (1-sqrt(3-6*LrQ))/(1+sqrt(3-6*LrQ) ) )*100;
ramI=conj((rLoadP+1j*rLoadQ)./mVoltP)+conj((rLoadP.*eLoadP+1j*rLoadQ.*eLoadQ)./mVoltP);
ramFI=Tp2f*conj((rLoadP+1j*rLoadQ)./mVoltP)+Tp2f*conj((rLoadP.*eLoadP+1j*rLoadQ.*eLoadQ)./mVoltP);
count=9;
plotF0=zeros(count,1);
plotF1=zeros(count,1);
plotF2=zeros(count,1);
for I=1:count
    usedLoadP=rLoadP;
    usedLoadQ=rLoadQ;
%     usedLoadP(3)=rLoadP(3)*(1-I/10);
%     usedLoadQ(2)=rLoadQ(2)*(1-I/10);
    usedVoltP=rVoltP;
    usedVoltP(2)=rVoltP(2)*(1-I/10);
    usedVoltP(3)=rVoltP(3)*(1-I/10);
    t=Tp2f*conj((usedLoadP+1j*usedLoadQ)./usedVoltP);
    f0=t(1);
    f1=t(2);
    f2=t(3);
    plotF0(I)=f0;
    plotF1(I)=f1;
    plotF2(I)=f2;
end
subplot(6,1,1);
plot(1:count,real(plotF0));
subplot(6,1,2);
plot(1:count,imag(plotF0));
subplot(6,1,3);
plot(1:count,real(plotF1));
subplot(6,1,4);
plot(1:count,imag(plotF1));
subplot(6,1,5);
plot(1:count,real(plotF2));
subplot(6,1,6);
plot(1:count,imag(plotF2));
% abs(t)