function ddh=func_ddh(Volt,Init_Y,Busnum,Y,UAngel,r,c,Angle,Loadi,ContrlCount)
mat_AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum);
mat_INV_AngleIJ=mat_AngleIJ';
yP=Init_Y(1:length(Init_Y)/2);
yQ=Init_Y(length(Init_Y)/2+1:end);
t1=-sparse(1:Busnum,1:Busnum,Y.*cos(mat_INV_AngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP,Busnum,Busnum);
t2=sparse(1:Busnum,1:Busnum,sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP,Busnum,Busnum)*Y.*cos(mat_AngleIJ);
t3=(t1+t2)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum);
t4=-(sparse(1:Busnum,1:Busnum,Y.*cos(mat_AngleIJ)*Volt,Busnum,Busnum) -sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*Y.*cos(mat_INV_AngleIJ) )*sparse(1:Busnum,1:Busnum,sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP,Busnum,Busnum);
ddPdTdT=t3+t4;%ok1
t1=(-sparse(1:Busnum,1:Busnum,Y.*sin(mat_AngleIJ)*Volt,Busnum,Busnum)+sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*Y.*sin(mat_INV_AngleIJ) )*sparse(1:Busnum,1:Busnum,yP,Busnum,Busnum);
t2= -sparse(1:Busnum,1:Busnum, sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP ,Busnum,Busnum)*Y.*sin(mat_AngleIJ)+sparse(1:Busnum,1:Busnum,Y.*sin(mat_INV_AngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP,Busnum,Busnum);
ddPdVdT=t1+t2;%ok1
t1=diag( Y.*sin(mat_INV_AngleIJ)*diag(Volt)*yP);
t2=diag(yP)*Y.*sin(mat_AngleIJ)*diag(Volt);
t3=-diag(yP)*diag(Y.*sin(mat_AngleIJ)*Volt);
t4=-Y.*sin(mat_INV_AngleIJ)*diag( diag(Volt)*yP );
ddPdTdV=t1+t2+t3+t4;
t1=Y.*cos(mat_INV_AngleIJ)*diag(yP);
t2=diag(yP)*Y.*cos(mat_AngleIJ);
ddPdVdV=t1+t2;
t1=-diag(Y.*sin(mat_AngleIJ)*Volt);
t2=diag(Volt)*Y.*sin(mat_INV_AngleIJ);
t3=(t1+t2)*diag( diag(Volt)*yQ );
t4=-diag( diag(Volt)*yQ )*Y.*sin(mat_AngleIJ);
t5=diag(Y.*sin(mat_INV_AngleIJ)*diag(Volt)*yQ);
t6=-(t4+t5)*diag(Volt);
ddQdTdT=t3+t6;%ok1
t1=(diag(Y.*cos(mat_AngleIJ)*Volt)-diag(Volt)*Y.*cos(mat_INV_AngleIJ) )*diag(yQ);
t2=+diag( diag(Volt)*yQ )*Y.*cos(mat_AngleIJ)-diag(Y.*cos(mat_INV_AngleIJ)*diag(Volt)*yQ);
ddQdVdT=t1+t2;
t1=Y.*cos(mat_INV_AngleIJ)*diag(diag(Volt)*yQ);
t2=diag(yQ)*diag(Y.*cos(mat_AngleIJ)*Volt);
t3=-diag(Y.*cos(mat_INV_AngleIJ)*diag(Volt)*yQ);
t4=-diag(yQ)*Y.*cos(mat_AngleIJ)*diag(Volt);
ddQdTdV=t1+t2+t3+t4;
t1=Y.*sin(mat_INV_AngleIJ)*diag(yQ);
t2=diag(yQ)*Y.*sin(mat_AngleIJ);
ddQdVdV=t1+t2;
t=[ddPdVdV+ddQdVdV,ddPdVdT+ddQdVdT;
   ddPdTdV+ddQdTdV,ddPdTdT+ddQdTdT;
];
sizeLoadi=size(Loadi,1)*2;
ddh=[
    sparse(sizeLoadi,ContrlCount);
    sparse(2*Busnum,sizeLoadi),-t;
    ];
end