function ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y)
%决定用循环重写
ContrlCount=size(PVi,1)+size(PGi,1)+Busnum*2;
%deltaPi/deltaThytai_deltaThytaj 非对角元素
dPidTidTj=zeros(Busnum);
t1=-Volt'*Volt;
t2=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
t3=t1.*t2;
for I=1:Busnum
    for J=1:Busnum
        dPidTidTj(I,J)=dPidTidTj(I,J)+Init_Y(2*I-1)*t3(I,J);
    end
end
t2=real(GB).*cos(AngleIJMat)-imag(GB).*sin(AngleIJMat);
t3=t1.*t2;
dPjdTidTj=zeros(Busnum);

for I=1:Busnum
    for J=1:Busnum
        dPjdTidTj(I,J)=dPjdTidTj(I,J)+Init_Y(2*J-1)*t3(I,J);
    end
end
dPdTidTj=dPidTidTj+dPjdTidTj;%最终非对角元素



 %% deltaP/deltaThyta_deltaThyta 对角元素
t1=Volt'*Volt;
t2=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
t3=t1.*t2;
t4=t3-diag(diag(t3));
t5=sum(t4,2);
t6=t5'.*Init_Y(1:2:size(Init_Y,2));%乘y的系数
dPidTidTi=diag(t6);
t1=Volt'*Volt;
t2=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
t2=t1.*t2;
dPidTjdTj=zeros(Busnum);
for J=1:Busnum
    for I=1:Busnum
        
            if I==J
                continue;
            end
            dPidTjdTj(J,J)=dPidTjdTj(J,J)+Init_Y(2*I-1)*t2(I,J);
       
    end
end

dPdTidTi=dPidTidTi+dPidTjdTj;
%%%%%%%%%
dPdTidTj=dPdTidTj-diag(diag(dPdTidTj));
hh=dPdTidTj+dPdTidTi; %#ok<NASGU>
 %% deltaP/deltaThytai_dVi 对角元素
 t1=ones(Busnum,1)*Volt.*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat));
t2=diag(t1);
t3=t1-diag(t2);%去掉了对角元素的
t4=sum(t3,2);
t4=t4'.*Init_Y(1:2:size(Init_Y,2));%乘y的系数
dPidTidVi=diag(t4);

t1=-Volt'*ones(1,Busnum).*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat));
dPjdTidVi=zeros(Busnum);
for J=1:Busnum
    for I=1:Busnum
        
            if I==J
                continue;
            end
            dPjdTidVi(J,J)=dPjdTidVi(J,J)+Init_Y(2*I-1)*t1(I,J);
       
    end
end
dPdTidVi=dPidTidVi+dPjdTidVi;
%% deltaP/deltaThytai_dVj 非对角元素

t1=Volt'*ones(1,Busnum).*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat));%20111224 存疑与学姐给的公式不一致
%t1=ones(1,Busnum)'*Volt.*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat));
dPidTidVj=zeros(Busnum);
for I=1:Busnum
    for J=1:Busnum
        dPidTidVj(I,J)=dPidTidVj(I,J)+Init_Y(2*I-1)*t1(I,J);
    end
end

dPjdTidVj=zeros(Busnum);
t1=Volt'*ones(1,Busnum).*(real(GB).*sin(AngleIJMat)+imag(GB).*cos(AngleIJMat));
for I=1:Busnum
    for J=1:Busnum
        dPjdTidVj(I,J)=dPjdTidVj(I,J)+Init_Y(2*J-1)*t1(I,J);
    end
end

dPdTidVj=dPidTidVj+dPjdTidVj;
%%%%%%%%%%%%%%%%%%%%%%
dPdTidVj=dPdTidVj-diag(diag(dPdTidVj));
hh=dPdTidVj+dPdTidVi; %#ok<NASGU>
%% deltaP/dVi_deltaThytaj 非对角元素
t1=-ones(Busnum,1)*Volt.*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat));
dPidVidTj=zeros(Busnum);
for I=1:Busnum
    for J=1:Busnum
        dPidVidTj(I,J)=dPidVidTj(I,J)+Init_Y(2*I-1)*t1(I,J);
    end
end
t1=-ones(Busnum,1)*Volt.*(real(GB).*sin(AngleIJMat)+imag(GB).*cos(AngleIJMat));
dPjdVidTj=zeros(Busnum);
for I=1:Busnum
    for J=1:Busnum
        dPjdVidTj(I,J)=dPjdVidTj(I,J)+Init_Y(2*J-1)*t1(I,J);
    end
end
dPdVidTj=dPidVidTj+dPjdVidTj;
%% deltaPi/dVi_deltaThyta 对角元素

dPdVidTi=dPdTidVi;
%%%%%%%%%%%%%%%%
dPdVidTj=dPdVidTj-diag(diag(dPdVidTj));
hh=dPdVidTj+dPdVidTi;
%% deltaP/dVi_dVj 非对角元素
t1=-(real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat));
dPidVidVj=zeros(Busnum);
for I=1:Busnum
    for J=1:Busnum
        dPidVidVj(I,J)=dPidVidVj(I,J)+Init_Y(2*I-1)*t1(I,J);
    end
end


t1=-(real(GB).*cos(AngleIJMat)-imag(GB).*sin(AngleIJMat));
dPjdVidVj=zeros(Busnum);
for I=1:Busnum
    for J=1:Busnum
        dPjdVidVj(I,J)=dPjdVidVj(I,J)+Init_Y(2*J-1)*t1(I,J);
    end
end
dPdVidVj=dPidVidVj+dPjdVidVj;
 %% deltaP/dVi_dVi 对角元素
t0=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
t1=diag(t0);
t2=t1'.*Init_Y(1:2:size(Init_Y,2));%乘y的系数
dPidVidVi=-2*diag(t2);
dPidVjdVj=0;
dPdVidVi=dPidVidVi+dPidVjdVj;
%%%%%%%%%%%
dPdVidVj=dPdVidVj-diag(diag(dPdVidVj));
hh=dPdVidVj+dPdVidVi;
%% 生成APi
APi=zeros(2*Busnum,2*Busnum);
dPdTidTj=dPdTidTj-diag(diag(dPdTidTj));
dPdTidVj=dPdTidVj-diag(diag(dPdTidVj));
dPdVidTj=dPdVidTj-diag(diag(dPdVidTj));
dPdVidVj=dPdVidVj-diag(diag(dPdVidVj));

APi(1:2:2*Busnum,1:2:2*Busnum)=dPdTidTj;%%非对角 TT
APi(1:2:2*Busnum,2:2:2*Busnum)=dPdTidVj;%%非对角 TV
APi(2:2:2*Busnum,1:2:2*Busnum)=dPdVidTj;%%非对角 VT
APi(2:2:2*Busnum,2:2:2*Busnum)=dPdVidVj;%%非对角 VV

APi(1:2:2*Busnum,1:2:2*Busnum)=APi(1:2:2*Busnum,1:2:2*Busnum)+dPdTidTi;%%对角
APi(1:2:2*Busnum,2:2:2*Busnum)=APi(1:2:2*Busnum,2:2:2*Busnum)+dPdTidVi;%%对角
APi(2:2:2*Busnum,1:2:2*Busnum)=APi(2:2:2*Busnum,1:2:2*Busnum)+dPdVidTi;%%对角
APi(2:2:2*Busnum,2:2:2*Busnum)=APi(2:2:2*Busnum,2:2:2*Busnum)+dPdVidVi;%%对角
%% deltaQ/deltaThyta_deltaThyta 非对角元素
t1=-Volt'*Volt;
% %t1=Volt'*Volt;
t2=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat);
t3=t1.*t2;
dQidTidTj=zeros(Busnum);
for I=1:Busnum
    for J=1:Busnum
        dQidTidTj(I,J)=dQidTidTj(I,J)+Init_Y(2*I)*t3(I,J);
    end
end

t2=real(GB).*sin(AngleIJMat)+imag(GB).*cos(AngleIJMat);
t3=-t1.*t2;

dQjdTidTj=zeros(Busnum);
for I=1:Busnum
    for J=1:Busnum
        dQjdTidTj(I,J)=dQjdTidTj(I,J)+Init_Y(2*J)*t3(I,J);
    end
end

dQdTidTj=dQidTidTj+dQjdTidTj;
 %% deltaQ/deltaThyta_deltaThyta 对角元素
t1=Volt'*Volt;
t2=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat);
t3=t1.*t2;
t4=t3-diag(diag(t3));
t5=sum(t4,2);
t6=t5'.*Init_Y(2:2:size(Init_Y,2));
dQidTidTi=diag(t6);
t1=Volt'*Volt;
 t2=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat);
 t3=t1.*t2;
 dQjdTidTi=zeros(Busnum);
 for J=1:Busnum
    for I=1:Busnum
        
            if I==J
                continue;
            end
            dQjdTidTi(J,J)=dQjdTidTi(J,J)+Init_Y(2*I)*t3(I,J);
       
    end
 end

dQdTidTi=dQjdTidTi+dQidTidTi;
%%%%%%%%%%%%
dQdTidTj=dQdTidTj-diag(diag(dQdTidTj));
hh=dQdTidTj+dQdTidTi;
%% deltaQ/deltaThyta_deltaV 非对角元素
t1=-Volt;
t2=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
t3=t1'*ones(1,Busnum).*t2;
dQidTidVj=zeros(Busnum);
for I=1:Busnum
    for J=1:Busnum
        dQidTidVj(I,J)=dQidTidVj(I,J)+Init_Y(2*I)*t3(I,J);
    end
end
t2=real(GB).*cos(AngleIJMat)-imag(GB).*sin(AngleIJMat);
t3=Volt'*ones(1,Busnum).*t2;
dQjdTidVj=zeros(Busnum);
for I=1:Busnum
    for J=1:Busnum
        %dQjdTidVj(I,J)=dQidTidVj(I,J)+Init_Y(2*J)*t3(I,J); 20111225
        dQjdTidVj(I,J)=dQjdTidVj(I,J)+Init_Y(2*J)*t3(I,J);
    end
end

dQdTidVj=dQidTidVj+dQjdTidVj;
 %% deltaQ/deltaThyta_deltaV 对角元素

t1=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
t2=-ones(Busnum,1)*Volt.*t1;
t2=t2-diag(diag(t2));
t3=sum(t2,2);
t3=t3'.*Init_Y(2:2:size(Init_Y,2));
dQidTidVi=diag(t3);
 


t1=Volt'*ones(1,Busnum).*(real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat));
t2=t1-diag(diag(t1));
dQjdTidVi=zeros(Busnum);
for J=1:Busnum
    for I=1:Busnum
        
            if I==J
                continue;
            end
            dQjdTidVi(J,J)=dQjdTidVi(J,J)+Init_Y(2*I)*t2(I,J);
       
    end
end

dQdTidVi=dQidTidVi+dQjdTidVi;
dQdTidVj=dQdTidVj-diag(diag(dQdTidVj));
hh=dQdTidVj+dQdTidVi;
%% deltaQ/deltaV_deltaV 非对角元素
t1=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat);
t2=-t1;
dQidVidVj=zeros(Busnum);
for I=1:Busnum
    for J=1:Busnum
        dQidVidVj(I,J)=dQidVidVj(I,J)+Init_Y(2*I)*t2(I,J);
    end
end

 t1=real(GB).*sin(AngleIJMat)+imag(GB).*cos(AngleIJMat);
dQjdVidVj=zeros(Busnum);
for I=1:Busnum
    for J=1:Busnum
        dQjdVidVj(I,J)=dQjdVidVj(I,J)+Init_Y(2*J)*t1(I,J);
    end
end

dQdVidVj=dQidVidVj+dQjdVidVj;
%% deltaQ/deltaV_deltaV 对角元素
t1=-2*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat));
t2=diag(t1);
t3=t2'.*Init_Y(2:2:size(Init_Y,2));%乘y的系数
dQidVidVi=diag(t3);

dQjdVidVi=0;
dQdVidVi=dQidVidVi+dQjdVidVi; 
%%%%%%%%%%%%%%
dQdVidVj=dQdVidVj-diag(diag(dQdVidVj));
hh=dQdVidVi+dQdVidVj;
 %% deltaQ/deltaV_deltaThyta 非对角元素

t1=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
t2=ones(Busnum,1)*Volt.*(t1);
dQidVidTj=zeros(Busnum);
for I=1:Busnum
    for J=1:Busnum
        dQidVidTj(I,J)=dQidVidTj(I,J)+Init_Y(2*I)*t2(I,J);
    end
end


t1=real(GB).*cos(AngleIJMat)-imag(GB).*sin(AngleIJMat);
t2=-ones(Busnum,1)*Volt.*(t1);
dQjdVidTj=zeros(Busnum);
for I=1:Busnum
    for J=1:Busnum
        dQjdVidTj(I,J)=dQjdVidTj(I,J)+Init_Y(2*J)*t2(I,J);
    end
end

 dQdVidTj=dQidVidTj+dQjdVidTj;
%% deltaQ/deltaV_deltaThyta 对角元素
dQdVidTi=dQdTidVi;% @
dQdVidTj=dQdVidTj-diag(diag(dQdVidTj));
hhd=dQdVidTi+dQdVidTj;
%% 生成AQi
AQi=zeros(2*Busnum,2*Busnum);
dQdTidTj=dQdTidTj-diag(diag(dQdTidTj));
dQdTidVj=dQdTidVj-diag(diag(dQdTidVj));
dQdVidTj=dQdVidTj-diag(diag(dQdVidTj));
dQdVidVj=dQdVidVj-diag(diag(dQdVidVj));

AQi(1:2:2*Busnum,1:2:2*Busnum)=dQdTidTj;%%非对角 TT
AQi(1:2:2*Busnum,2:2:2*Busnum)=dQdTidVj;%%非对角 TV
AQi(2:2:2*Busnum,1:2:2*Busnum)=dQdVidTj;%%非对角 VT
AQi(2:2:2*Busnum,2:2:2*Busnum)=dQdVidVj;%%非对角 VV

AQi(1:2:2*Busnum,1:2:2*Busnum)=AQi(1:2:2*Busnum,1:2:2*Busnum)+dQdTidTi;%%对角
AQi(1:2:2*Busnum,2:2:2*Busnum)=AQi(1:2:2*Busnum,2:2:2*Busnum)+dQdTidVi;%%对角
AQi(2:2:2*Busnum,1:2:2*Busnum)=AQi(2:2:2*Busnum,1:2:2*Busnum)+dQdVidTi;%%对角
AQi(2:2:2*Busnum,2:2:2*Busnum)=AQi(2:2:2*Busnum,2:2:2*Busnum)+dQdVidVi;%%对角
%% 生成ddh
t=[zeros(size(PGi,1)+size(PVi,1),ContrlCount);
    zeros(2*Busnum,size(PVi,1)+size(PGi,1)),AQi+APi;
];
ddh=t;
%% 以下是学姐给的公式
AngleIJ=AngleIJMat-angle(GB);
%yP=Init_Y(1:2:size(Init_Y,2));
yP=Init_Y(1:2:size(Init_Y,2));%暂时改这里  20111227
%yQ=Init_Y(2:2:size(Init_Y,2));
t1=-diag(Y.*cos(AngleIJ')*diag(Volt)*yP');
t2=diag(diag(Volt)*yP')*Y.*cos(AngleIJ);
t3=(t1+t2)*diag(Volt);
t4=-(diag(Y.*cos(AngleIJ)*Volt') -diag(Volt)*Y.*cos(AngleIJ') )*diag(diag(Volt)*yP');
ddPdTdT=t3+t4;%ok1
tttt=t2*diag(Volt);%ok1
ttttt=diag(Volt)*Y.*cos(AngleIJ')*diag(diag(Volt)*yP');%ok1
tttttt=tttt+ttttt;%ok1
ttttttt=-diag(Y.*cos(AngleIJ)*Volt')*diag(diag(Volt)*yP')+tttt;%ok1
tttttttt=diag(Volt)*Y.*cos(AngleIJ')*diag(diag(Volt)*yP')+t1*diag(Volt);%ok1
t1=(-diag(Y.*sin(AngleIJ)*Volt')+diag(Volt)*Y.*sin(AngleIJ') )*diag(yP);
t2= -diag( diag(Volt)*yP' )*Y.*sin(AngleIJ)+diag(Y.*sin(AngleIJ')*diag(Volt)*yP');
ddPdVdT=t1+t2;%ok1
tttt=-diag( diag(Volt)*yP' )*Y.*sin(AngleIJ);
t1=diag( Y.*sin(AngleIJ')*diag(Volt)*yP');
t2=diag(yP)*Y.*sin(AngleIJ)*diag(Volt);
t3=-diag(yP)*diag(Y.*sin(AngleIJ)*Volt');
t4=-Y.*sin(AngleIJ')*diag( diag(Volt)*yP' );
ddPdTdV=t1+t2+t3+t4;%存疑与我的不一样
tttt=t2;
ttttt=t4;
t1=Y.*cos(AngleIJ')*diag(yP);
t2=diag(yP)*Y.*cos(AngleIJ);
ddPdVdV=t1+t2;
t1=-diag(Y.*sin(AngleIJ)*Volt');
t2=diag(Volt)*Y.*sin(AngleIJ');
t3=(t1+t2)*diag( diag(Volt)*yQ' );
t4=-diag( diag(Volt)*yQ' )*Y.*sin(AngleIJ);

t5=diag(Y.*sin(AngleIJ')*diag(Volt)*yQ');
t6=-(t4+t5)*diag(Volt);
ddQdTdT=t3+t6;%ok1
tttt=-(t4)*diag(Volt);
ttttt=t2*diag( diag(Volt)*yQ' );
tttttt=t1*diag( diag(Volt)*yQ' )+tttt;
ttttttt=-t5*diag(Volt)+t2*diag( diag(Volt)*yQ' );
t1=(diag(Y.*cos(AngleIJ)*Volt')-diag(Volt)*Y.*cos(AngleIJ') )*diag(yQ);
t2=+diag( diag(Volt)*yQ' )*Y.*cos(AngleIJ)-diag(Y.*cos(AngleIJ')*diag(Volt)*yQ');
tttt=diag( diag(Volt)*yQ' )*Y.*cos(AngleIJ);
ttttt=-diag(Volt)*Y.*cos(AngleIJ') *diag(yQ);
ddQdVdT=t1+t2;
t1=Y.*cos(AngleIJ')*diag(diag(Volt)*yQ');
t2=diag(yQ)*diag(Y.*cos(AngleIJ)*Volt');
t3=-diag(Y.*cos(AngleIJ')*diag(Volt)*yQ');
t4=-diag(yQ)*Y.*cos(AngleIJ)*diag(Volt);
ddQdTdV=t1+t2+t3+t4;
t1=Y.*sin(AngleIJ')*diag(yQ);
t2=diag(yQ)*Y.*sin(AngleIJ);
ddQdVdV=t1+t2;
%%%%
t=zeros(2*Busnum);
% t(1:2:2*Busnum,1:2:2*Busnum)=ddPdTdT+ddQdTdT;
% %t(1:2:2*Busnum,2:2:2*Busnum)=ddPdTdV+ddQdTdV;
% %t(2:2:2*Busnum,1:2:2*Busnum)=ddPdVdT+ddQdVdT;
% t(1:2:2*Busnum,2:2:2*Busnum)=ddPdVdT+ddQdVdT;
% t(2:2:2*Busnum,1:2:2*Busnum)=ddPdTdV+ddQdTdV;
% t(2:2:2*Busnum,2:2:2*Busnum)=ddPdVdV+ddQdVdV;暂时改一下  20111227
t=[ddPdVdV+ddQdVdV,ddPdVdT+ddQdVdT;
    ddPdTdV+ddQdTdV,ddPdTdT+ddQdTdT;
];
t=[zeros(size(PGi,1)+size(PVi,1),ContrlCount);
    zeros(2*Busnum,size(PVi,1)+size(PGi,1)),-t;
];
ddh=t;
end