diff --git a/src/CUDA/CudaCollimatorPhysics.cu b/src/CUDA/CudaCollimatorPhysics.cu
index b44c360..4aca093 100644
--- a/src/CUDA/CudaCollimatorPhysics.cu
+++ b/src/CUDA/CudaCollimatorPhysics.cu
@@ -42,6 +42,23 @@ __device__ inline double3 cross(double3 &lhs, double3 &rhs) {
   return tmp;
 }
 
+__device__ inline double3 ArbitraryRotation(double3 &W, double3 &Rorg, double Theta) {
+  double c=cos(Theta);
+  double s=sin(Theta);
+  double dotW = sqrt(dot(W,W));
+  W.x = W.x / dotW;
+  W.y = W.y / dotW;
+  W.z = W.z / dotW;
+
+  double dotWR = dot(W, Rorg) * (1.0 - c);
+  double3 crossW = cross(W, Rorg);
+  double3 tmp;
+  tmp.x = Rorg.x * c + crossW.x * s + W.x * dotWR;
+  tmp.y = Rorg.y * c + crossW.y * s + W.y * dotWR;
+  tmp.z = Rorg.z * c + crossW.z * s + W.z * dotWR;
+  return tmp;
+} 
+
 __device__ inline bool checkHit(double &z, double *par) {
 
   /* check if particle is in the degrader material */
@@ -101,6 +118,7 @@ __device__ inline void Rot(double &px, double &pz, double &x, double &z, double
   double Psixz;
   double pxz;
 
+  /*
   if (px>=0 && pz>=0)
     Psixz = atan(px/pz);
   else if (px>0 && pz<0)
@@ -109,7 +127,8 @@ __device__ inline void Rot(double &px, double &pz, double &x, double &z, double
     Psixz = atan(px/pz) + 2*PI;
   else
     Psixz = atan(px/pz) + PI;
-
+  */
+  Psixz = atan2(px, pz);
   pxz = sqrt(px*px + pz*pz);
 
   if(coord==1) {
@@ -151,20 +170,9 @@ __device__ inline void coulombScat(double3 &R, double3 &P, curandState &state, d
   }
 
   //__syncthreads();  
-
   double xplane = z1 * deltas * theta0 / sqrt(12.0) + z2 * deltas * theta0 / 2.0;
   Rot(P.x, P.z, R.x, R.z, xplane, normP, thetacou, deltas, 1, par);
 
-  double P2 = curand_uniform_double(&state);//gsl_rng_uniform(rGen_m);
-  if( (P2 < 0.0047) && enableRutherfordScattering) {
-    double P3 = curand_uniform_double(&state);//gsl_rng_uniform(rGen_m);
-    double thetaru = 2.5 * sqrt(1 / P3) * sqrt(2.0) * theta0;
-    double P4 = curand_uniform_double(&state);//gsl_rng_uniform(rGen_m);
-    if(P4 > 0.5)
-      thetaru = -thetaru;
-    Rot(P.x,P.z,R.x,R.z, xplane, normP, thetaru, deltas, 0, par);
-  }
-
   // y-direction: See Physical Review, "Multiple Scattering"
   z1 = curand_normal_double(&state);//gsl_ran_gaussian(rGen_m,1.0);
   z2 = curand_normal_double(&state);//gsl_ran_gaussian(rGen_m,1.0);
@@ -181,14 +189,23 @@ __device__ inline void coulombScat(double3 &R, double3 &P, curandState &state, d
   double yplane = z1 * deltas * theta0 / sqrt(12.0) + z2 * deltas * theta0 / 2.0;
   Rot(P.y,P.z,R.y,R.z, yplane, normP, thetacou, deltas, 2, par);
 
-  P2 = curand_uniform_double(&state);//gsl_rng_uniform(rGen_m);
+  double P2 = curand_uniform_double(&state);//gsl_rng_uniform(rGen_m);
   if( (P2 < 0.0047) && enableRutherfordScattering) {
     double P3 = curand_uniform_double(&state);//gsl_rng_uniform(rGen_m);
-    double thetaru = 2.5 * sqrt(1 / P3) * sqrt(2.0) * theta0;
-    double P4 = curand_uniform_double(&state);//gsl_rng_uniform(rGen_m);
-    if(P4 > 0.5)
-      thetaru = -thetaru;
-    Rot(P.y,P.z,R.y,R.z, yplane, normP, thetaru, deltas, 0, par);
+    //double thetaru = 2.5 * sqrt(1 / P3) * sqrt(2.0) * theta0;
+    double thetaru = 2.5 * sqrt(1 / P3) * 2.0 * theta0;
+    double phiru = 2.0 * M_PI * curand_uniform_double(&state);
+    double th0=atan2(sqrt(P.x*P.x+P.y*P.y),fabs(P.z));
+    double3 W,X;
+    
+    double dotP = sqrt(dot(P,P));
+    X.x = cos(phiru)*sin(thetaru) * dotP;
+    X.y = sin(phiru)*sin(thetaru) * dotP;
+    X.z = cos(thetaru) * dotP;
+    W.x = -P.y;
+    W.y = P.x;
+    W.z = 0.0;
+    P = ArbitraryRotation(W, X, th0);
   }
 
 }