Need help expanding particle system spread / divergence from 2 to 3 dimensions

I need help. I've been struggling with this for a week now and getting nowhere. I am building a 3D particle system mainly for learning and I am currently working on particle spread / divergence. In specific, introducing random direction to the particle direction so as to create something that looks more like a fountain as opposed to a solid stream.

I have been successful in getting this to work in one axis but no matter what I do, I cannot get it to work in 3 dimensions.

Here is what I am doing:

    // Compute a random angle between -180 to +180 for velocity angle x, y and z. spreadAmount is a float from 0.0 to 1.0 to control degree of spread.
    float velangrndx = spreadAmount * ((((double)(rand() % RAND_MAX) / (RAND_MAX)) - 0.5) * 360.0 * 3.14159265359 / 180.0);
    float velangrndy = spreadAmount * ((((double)(rand() % RAND_MAX) / (RAND_MAX)) - 0.5) * 360.0 * 3.14159265359 / 180.0);
    float velangrndz = spreadAmount * ((((double)(rand() % RAND_MAX) / (RAND_MAX)) - 0.5) * 360.0 * 3.14159265359 / 180.0);

// Compute Angles
float vsin_anglex_dir = -PF_SIN(velangrndx);
float vcos_anglex_dir = -PF_COS(velangrndx);
float vsin_angley_dir = -PF_SIN(velangrndy);
float vcos_angley_dir = -PF_COS(velangrndy);
float vsin_anglez_dir = -PF_SIN(velangrndz);
float vcos_anglez_dir = -PF_COS(velangrndz);

And finally, we get to the particle spread / divergence function below. If I use only the first X axis (comment out the Y and Z) it works as it should (see images), but if I use the Y and Z axis, it works totally incorrectly. px0, py0, and pz0 are temporary float variables so as to preserve the velocity variables.

// X Divergence
px0 = (velx * vsin_anglex_dir);
py0 = (velx * vcos_anglex_dir);
pz0 = velz;
velx = px0; vely = py0; velz = pz0;

// Y Divergence
py0 = (vely * vsin_angley_dir);
pz0 = (vely * vcos_angley_dir);
px0 = velx;
velx = px0; vely = py0; velz = pz0;

// Z Divergence
pz0 = (velz * vsin_anglez_dir);
px0 = (velz * vcos_anglez_dir);
py0 = vely;
velx = px0; vely = py0; velz = pz0;

The velx, vely, and velz are then used to calculate for particle screen position.

This is what the particle spread looks like at 25%, 75% and 100% for the X axis only (if I comment out the Y and Z code). This works as it should and in theory, if the rest of my code was working correctly, I should get this same result for the Y and Z axis. But I don't.

I could really use some help here. Any suggestions on what I am doing wrong and how to correctly expand the currently working spread function from 2 dimensions to 3?

Thanks, -Richard

Likely it is because the values of velx, vely and velz are getting overwritten on subsequent calculations. See whether the below works the way you are expecting.

// X Divergence
float velxXD = (velx * vsin_anglex_dir);
float velyXD = (velx * vcos_anglex_dir);
float velzXD = velz;

// Y Divergence
float velxYD = velx;
float velyYD = (vely * vsin_angley_dir);
float velzYD = (vely * vcos_angley_dir);

// Z Divergence
float velxZD = (velz * vcos_anglez_dir); 
float velyZD = vely;
float velzZD = (velz * vsin_anglez_dir);

velx=velxXD+velxYD+velxZD;
vely=velyXD+velyYD+velyZD;
velz=velzXD+velzYD+velzZD;