As Daniel points out, you will always get a good speed-up by compiling. I want to take a different approach and consider how you can get a faster execution by writing your code in a more "functional" style.

My machine is slower than yours and ran your code in 17.2 seconds, according to AbsoluteTiming.

The first thing I noticed is that you are creating and overwriting the same variable definitions a lot in the body of your loop. This often slows things down.

The first thing to consider is creating all the random numbers in a single tensor instead of as two vectors that are re-defined each time.

starter = 
   Transpose[{RandomReal[{0., Pi}, {numsteps, numparticles}], 
     RandomReal[{0, 2.*Pi}, {numsteps, numparticles}]}, {2, 1, 
     3}]; // Timing

{0.885516, Null}

This creates a numsteps by 2 by numparticles tensor.

I'm going to set your loop up as a Fold that starts with the starting values particlesx etc and "folds" in the random numbers in each row of starter to update the values determined in the previous step.

The other thing I've done is change the function to determine if a crossing has taken place to use Sign rather than If in a Map. There is in any case a function Boole which does the If[something, 1, 0] construct, and it's Listable. Unfortunately Greater and Less aren't Listable so that's not an option.

Fold[With[{newx = #1[[1]] + 0.01*Sin[#2[[1]]]*Cos[#2[[2]]], 
      newy = #1[[2]] + 0.01*Sin[#2[[1]]]*Sin[#2[[2]]], 
      newz = #1[[3]] + 0.01*Cos[#2[[1]]], 
      rold = #1[[1]]^2 + #1[[2]]^2 + #1[[3]]^2}, 
  {newx, newy, newz, 
      0.5 - 0.5 Sign[(rold - radius)*(newx^2 + newy^2 + newz^2 - 
            radius)] + #1[[4]]}] &, 
 {particlesx, particlesy, particlesz, numcrossings}, starter]; // AbsoluteTiming

{11.741202, Null}

Just this change in programming style gives a 30% speed-up, without compiling.

An alternative would be to use UnitStep in place of Sign, like this, but on my machine this gives only a small speedup:

Fold[With[{newx = #1[[1]] + 0.01*Sin[#2[[1]]]*Cos[#2[[2]]], 
      newy = #1[[2]] + 0.01*Sin[#2[[1]]]*Sin[#2[[2]]], 
      newz = #1[[3]] + 0.01*Cos[#2[[1]]], 
      rold = #1[[1]]^2 + #1[[2]]^2 + #1[[3]]^2}, {newx, newy, newz, 
      UnitStep[-(rold - radius)*(newx^2 + newy^2 + newz^2 - 
           radius)] + #1[[4]]}] &, {particlesx, particlesy, 
    particlesz, numcrossings}, starter]; // AbsoluteTiming

{15.007808, Null}

You might also be able to speed up the loop by constructing starter as the pre-calculated trig values

trigstarter = 
   0.01*Transpose[
     With[{theta = RandomReal[{0., Pi}, {numsteps, numparticles}], 
       phi = RandomReal[{0, 2.*Pi}, {numsteps, numparticles}]}, 
   {Sin[theta]*Cos[phi], Sin[theta]*Sin[phi], Cos[theta]}], {2, 1, 3}]; // AbsoluteTiming

{3.010354, Null}

And then the Fold is even faster:

result=Fold[With[{newx = #1[[1]] + #2[[1]], newy = #1[[2]] + #2[[2]], 
      newz = #1[[3]] + #2[[3]], 
      rold = #1[[1]]^2 + #1[[2]]^2 + #1[[3]]^2}, {newx, newy, newz, 
      0.5 - 0.5 Sign[(rold - radius)*(newx^2 + newy^2 + newz^2 - 
            radius)] + #1[[4]]}] &, 
 {particlesx, particlesy, particlesz, numcrossings}, trigstarter]; // AbsoluteTiming

{2.495747, Null}

This is a major speedup. By taking the trigonometric functions out of the loop, you are only pre-compiling them internally once.

You can then do things like:

Histogram[Last@result]