I have a set of points that I would like to express with respect to one of them. This is the set of points:

I would like to express all the points with respect to cornerD. So I thought I would translate the points in order to bring cornerD to the origin and then rotate by the angle that the line (cornerA, cornerD) forms with the vertical axis. The points are actually in 3D expressed in homogeneous coordinates but they all lie in a plane. This is the matlab code:
%represent the points according to the 3D reference system in cornerD
%these are the points
cornerA3d = [cornerAI(1) cornerAI(2) 0 1]';
cornerB3d = [cornerBI(1) cornerBI(2) 0 1]';
cornerC3d = [cornerCI(1) cornerCI(2) 0 1]';
cornerD3d = [cornerDI(1) cornerDI(2) 0 1]';
peschiera3d = [peschieraI(1) peschieraI(2) 0 1]';
riva3d = [rivaI(1) rivaI(2) 0 1]';
salo3d = [saloI(1) saloI(2) 0 1]';
garda3d = [gardaI(1) gardaI(2) 0 1]';
iseo3d = [iseoI(1) iseoI(2) 0 1]';
lovere3d = [lovereI(1) lovereI(2) 0 1]';
%%%%%
%find angle theta to rotate the points
x = cornerA3d(1) - cornerD3d(1);
y = cornerA3d(2) - cornerD3d(2);
theta = atan2(x, y);
%rotate wrt z axis
Hrot = [
cos(theta) -sin(theta) 0 0;
sin(theta) cos(theta) 0 0;
0 0 1 0;
0 0 0 1
];
%translate to bring cornerD in the origin
Htrasl = [
1 0 0 -cornerD3d(1);
0 1 0 -cornerD3d(2);
0 0 1 -cornerD3d(3);
0 0 0 1
];
H = Hrot*Htrasl;
%new set of points
cornerA3d = H * cornerA3d;
cornerB3d = H * cornerB3d;
cornerC3d = H * cornerC3d;
cornerD3d = H * cornerD3d;
peschiera3d = H * peschiera3d;
riva3d = H * riva3d;
salo3d = H * salo3d;
garda3d = H * garda3d;
iseo3d = H * iseo3d;
lovere3d = H * lovere3d;
However this is the result:

I don't see where I went wrong... Any ideas? Thank you. :)