!Input: p_lab(3), q_lab(3) - real, momentum and polarization 3-vectors
subroutine transform_muon(p_lab, q_lab, p_out, q_out)

use precision_def
implicit none
real(rp) p_lab(3), q_lab(3), p_out(3), q_out(3)
real(rp) p_1(3), p_2(3), q_1(3), q_2(3)
real(rp) alpha, gamma, beta

!rotate about z by alpha to get p into the x-z plane
! the azimuthal angle is tan(alpha) = p_y/p_x

   alpha = atan2(p_lab(2),p_lab(1))
   call rotate(p_lab, 'z',alpha,p_1)
   call rotate(q_lab, 'z',alpha,q_1)

! then about y to get p along z
   gamma = atan2(p_1(1),p_1(3))
   call rotate(p_1, 'y',gamma,p_out)
   call rotate(q_1, 'y',gamma,q_2)

! then about gamm to get q in x-z plane
   beta = atan2(q_2(2),q_2(1))
   call rotate(q_2, 'z',beta, q_out)

return
end