subroutine signed_phi(svec,pvec, phi)

use sim_utils
implicit none
real(rp) svec(3), pvec(3),khat(3)
real(rp) phi, cosphi
real(rp) svec_pp(3), svec_p(3), svec_ppp(3), pvec_p(3), pvec_pp(3)
real(rp) theta_z, theta_y
real(rp) cross(3)
real(rp) R_z(3,3), R_y(3,3), R_zp(3,3)
real(rp) svec_mag, pvec_mag
real(rp) phi0


svec_mag = sqrt(dot_product(svec,svec))
pvec_mag = sqrt(dot_product(pvec,pvec))
cosphi = dot_product(svec,pvec)/svec_mag/pvec_mag

!rotate about z-axis so that momentum is in the x-z plane
theta_z = atan2(pvec(2),pvec(1))
call rot_matrix(theta_z, 'z', R_z)
pvec_p = matmul(R_z,pvec)
svec_p = matmul(R_z,svec)
!rotate about y so that momentum is in z-direction
theta_y=atan2(-pvec_p(1),pvec_p(3))
call rot_matrix(theta_y,'y',R_y)
pvec_pp = matmul(R_y,pvec_p)
svec_pp = matmul(R_y,matmul(R_z,svec) )

cross = cross_product(svec_pp,pvec_pp)

!rotate about z so that svec_pp is in x-z plane
theta_z = atan2(svec_pp(2),svec_pp(1))
if(theta_z <0)theta_z = theta_z+twopi/2
call rot_matrix(theta_z,'z',R_zp)
svec_ppp=matmul(R_zp,svec_pp)
!now pvec_pp is in the z direction, svec_ppp is in the x-z plane so that svec_ppp X pvec_pp is in the y direction
khat = cross_product(svec_ppp,pvec_pp)

phi = atan2(khat(2),dot_product(svec,pvec))*sign(1.0_rp,svec_ppp(1))
phi = acos(cosphi)*sign(1.0_rp, svec(1)) 
if(phi < 0) phi = phi+twopi
return
end subroutine signed_phi



!call cross_product(A,B,C)
!subroutine cross_product(A,B,C)
!  use precision_def
!  implicit none
!  real(rp) A(3),B(3),C(3)

!  C = (/A(2)*B(3)-A(3)*B(2),A(3)*B(1)-A(1)*B(3),A(1)*B(2)-A(2)*B(1)/)
!  return
!  end subroutine cross_product

subroutine rot_matrix(theta, axis, R)
use precision_def
implicit none

real(rp) theta,R(3,3)
character*1 axis
integer i

R(:,:)=0
forall(i=1:3)R(i,i)=1
if(axis == 'x')then
   R(2,2)= cos(theta)
   R(2,3)=sin(theta)
   R(3,2)=-R(2,3)
   R(3,3)=R(2,2)
endif
if(axis == 'y')then
   R(1,1)= cos(theta)
   R(1,3)=sin(theta)
   R(3,1)=-R(1,3)
   R(3,3)=R(1,1)
endif
if(axis == 'z')then
   R(1,1)= cos(theta)
   R(1,2)=sin(theta)
   R(2,1)=-R(1,2)
   R(2,2)=R(1,1)
endif

return
end subroutine rot_matrix