function R = rotationmatrix(phi)
R = [cos(phi), -sin(phi);
     sin(phi),  cos(phi)];