function [k]=Q4_2Dk(X, Y, t, E, Nu)
%clear all;                                   


D=E/(1-Nu^2)*[1 Nu 0;Nu 1 0; 0 0 (1-Nu)/2];

syms Xi Eta real;



XiC=[-1 1 1 -1];
EtaC=[-1 -1 1 1];

Ni=[ ( (XiC(2)-Xi)*(EtaC(4)-Eta) )/( (XiC(2)-XiC(1))*(EtaC(4)-EtaC(1))),    ((XiC(1)-Xi)*(EtaC(3)-Eta))/( (XiC(1)-XiC(2))*(EtaC(3)-EtaC(2))),      ((XiC(4)-Xi)*(EtaC(2)-Eta))/( (XiC(4)-XiC(3))*(EtaC(2)-EtaC(3))),      ((XiC(3)-Xi)*(EtaC(1)-Eta))/( (XiC(3)-XiC(4))*(EtaC(1)-EtaC(4))) ];

dNibydXi=diff(Ni,Xi);
dNibydEta=diff(Ni,Eta);




x=Ni*X';
y=Ni*Y';

J=[ diff(x,Xi), diff(y,Xi); diff(x,Eta), diff(y,Eta)];
Jinv=inv(J);

dNibydx=simplify(Jinv(1,1)*dNibydXi+Jinv(1,2)*dNibydEta);
dNibydy=simplify(Jinv(2,1)*dNibydXi+Jinv(2,2)*dNibydEta);

B=[dNibydx(1)      0      dNibydx(2)      0      dNibydx(3)      0        dNibydx(4)      0      ; ...
         0      dNibydy(1)      0      dNibydy(2)      0      dNibydy(3)      0      dNibydy(4);   ...
   dNibydy(1)   dNibydx(1) dNibydy(2)  dNibydx(2)  dNibydy(3) dNibydx(3)  dNibydy(4) dNibydx(4)];    

BTDBJ1=subs(B'*D*B*det(J)*t,{Xi, Eta},{ -1/sqrt(3), -1/sqrt(3)});
BTDBJ2=subs(B'*D*B*det(J)*t,{Xi, Eta},{ 1/sqrt(3), -1/sqrt(3)});
BTDBJ3=subs(B'*D*B*det(J)*t,{Xi, Eta},{ 1/sqrt(3), 1/sqrt(3)});
BTDBJ4=subs(B'*D*B*det(J)*t,{Xi, Eta},{ -1/sqrt(3), 1/sqrt(3)});
k=BTDBJ1+BTDBJ2+BTDBJ3+BTDBJ4;





   