clc
clear
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;16;12];
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);