How to expand particle system angle spread to 3 dimensions

I am trying to understand how to implement angle spread for a 3D particle system, to achieve an effect similar to a fountain. I can get it to work for a 2D system, but not a 3D. I would really appreciate any help as I've tried just about everything.

Here is what I'm doing:

Compute an initial random angle between -180 to +180. spreadAmount is a float from 0.0 to 1.0 to control degree of spread.

float velangrnd = spreadAmount * ((((double)(rand() % RAND_MAX) / (RAND_MAX)) - 0.5) * 360.0 * 3.14159265359 / 180.0);

Compute angles:

float vsin_anglex_dir = -SIN(velangrnd);
float vcos_anglex_dir = -COS(velangrnd);

And finally, calculate the angle spread. Vel is the speed from 0-1:

// XY Spread
float px0 = (vel * vsin_anglex_dir);
float py0 = (vel * vcos_anglex_dir);
float pz0 = 0;

After that, I simply compute the screen position. x, y, z are the emitter coordinates:

px0 = x + px0 * time;
py0 = y + py0 * time;
pz0 = z + pz0 * time;

This creates a perfect circle of particles in the XY axis when spreadAmount is 1.0. In other words, particles will shoot out in 360 degrees and according to the velocity (vel). At a lower value, it will create a 2D fountain effect.

However, not matter what I try, I cannot expand this to a second axis. I am trying to create a 3D fountain. So that one axis shoots particles outwards, another adds random spread angle in one direction, and the third a random spread angle in the other direction. Thus, a fountain.

Are there any suggestions on how to do this before I pull out the remainder of my hair?

Thank you!

in 3D you need main direction unit vector t and your angle a from that you create 2 basis vectors u,v that are perpendicular to each and to t and from that you construct your random directions... Something like this (if I did not do any silly mistake):

float a=?,b,r;
vec3 u,v,t=vec3(?,?,?),d;

// create u,v from t
u=vec3(1,0,0); // u is any non zero vector
if (fabs(dot(u,t))>0.75) u=vec3(0,1,0); // but not (anti)parallel to t
u=normalize(cross(u,t)); // make it perpendicular and unit
v=normalize(cross(u,t)); // make it perpendicular and unit

// compute random direction d
b=2.0*M_PI*Random();   // random angle
r=tan(0.5*a)*Random(); // random radius
d=normalize((r*u*cos(b))+(r*v*sin(b))+t);

So basically its random vector within cone where t is central axis and cone angle is a. I used my GLSL_math.h however you can use any vector math like GLM ...

[Edit1] lib less and component wise version

#include <math.h>
void normalize(float &x,float &y,float &z)
    {
    float l=sqrt((x*x)+(y*y)+(z*z));
    if (l>1e-6) l=1.0/l; else l=0.0;
    x*=l; y*=l; z*=l;
    }
void cross(float &x,float &y,float &z, float ax,float ay,float az, float bx,float by,float bz)
    {
    x=(ay*bz)-(az*by);
    y=(az*bx)-(ax*bz);
    z=(ax*by)-(ay*bx);
    }
void cone(float &dx,float &dy,float &dz, float tx,float ty,float tz, float a)
    {
    // (dx,dy,dz) <- random direction inside cone
    // (tx,ty,tz) -> cone axis direction
    // a          -> cone angle in [rad]
    float b,c,s,r;
    float ux,uy,uz;
    float vx,vy,vz;
    // create u,v from t
    ux=1.0; uy=0.0; uz=0.0;                 // u is any non zero vector
    if (fabs((ux*tx)+(uy*ty)+(uz*tz))>0.75) // but not (anti)parallel to t
        { ux=0.0; uy=1.0; uz=0.0; }
    cross(ux,uy,uz, ux,uy,uz, tx,ty,tz);    // make it perpendicular
    normalize(ux,uy,uz);                    // make it unit
    cross(vx,vy,vz, ux,uy,uz, tx,ty,tz);    // make it perpendicular
    normalize(vx,vy,vz);                    // make it unit
    // compute random direction d
    b=2.0*M_PI*Random();                    // random angle
    r=tan(0.5*a)*Random();                  // random radius
    c=r*cos(b); s=r*sin(b);
    dx=(ux*c)+(vx*s)+tx;                    // random direction inside cone
    dy=(uy*c)+(vy*s)+ty;
    dz=(uz*c)+(vz*s)+tz;
    normalize(dx,dy,dz);                    // make it unit
    }

If you interested the vector math equations (and even array version implementation) is in here (Edit2 near bottom):

• Understanding 4x4 homogenous transform matrices

Using array or even classes is much more convenient (less amount and much clearer code as you can see when you compare with previous version).

[Edit2] better version using sort of spherical coordinates allowing spread <0 , 2*PI> [rad] with "uniform" distribution

I was thinking about merging my cone approach with spherical coordinate system to avoid the rectangular cone cap ... Here the result:

vec3 rnd_dir(vec3 t,float spread)
    {
    vec3 b,n;
    float a,r,x,y,z;
    static const float  pih=0.5*M_PI;
    static const float _pih=2.0/M_PI;
    // (x,y,z) = unit spread where +x is main direction (central axis)
    x=Random();                 // random angle within spread scaled to <0,1>
    a=x*spread*0.5;             // random angle within spread [rad]
    if (a>pih)                  // linearize the point distribution (avoid high density on poles)
        {
        a=M_PI-a;
        a=sqrt(a*_pih)*pih;     // sqrt looks good probably because surface area is scaled with ^2
        a=M_PI-a;
        }
    else{
        a=sqrt(a*_pih)*pih;     // sqrt looks good probably because surface area is scaled with ^2
        }
    x=cos(a);                   // convert angle to main axis coordinate
    r=sqrt(1.0-(x*x));          // max radius of cone cap still inside unit sphere
    a=Random()*2.0*M_PI;        // random polar angle inside the cone cap [rad]
    y=r*cos(a);                 
    z=r*sin(a);
    // create normal n, binormal b from tangent t ... TBN matrix
//  t=normalize(t);             // t should be unit if it is not uncomment this line
    n=vec3(1,0,0);              // n is any non zero vector
    if (fabs(dot(n,t))>0.75) n=vec3(0,1,0); // but not (anti)parallel to t
    n=normalize(cross(n,t));    // make it perpendicular and unit
    b=normalize(cross(n,t));    // make it perpendicular and unit
    // convert (x,y,z) so t is main direction
    return (t*x)+(b*y)+(n*z);
    }

Its using vec3 again (I am too lazy to code in x,y,z) however the axis aligned spread generation itself is not using vectors at all so you can directly use that and use the conversion to t,b,n from previous version (t,u,v) its the same...

Here preview of all spreads with 5 deg step (1000 points):

rendered in OpenGL like this:

void gl_draw()
    {
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

    float aspect=float(xs)/float(ys);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    gluPerspective(60.0/aspect,aspect,0.1,100.0);
    glMatrixMode(GL_TEXTURE);
    glLoadIdentity();
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    glTranslatef(0.0,0.0,-5.5);
    static float a=0.0; a+=5.5; if (a>360.0) a-=360.0;
    a=-75.0;
    glRotatef(a,0.0,1.0,0.0);

    glEnable(GL_DEPTH_TEST);
    glDisable(GL_TEXTURE_2D);

    int i;
    vec3 d,t=normalize(vec3(1.0,0.5,-0.2));

    Randomize();
    RandSeed=0x1234567;
    glPointSize(1);

    glBegin(GL_POINTS);
    glColor3f(0.1,0.8,1.0);
    for (i=0;i<1000;i++)
        {
        d=rnd_dir(t,spread);
        glVertex3fv(d.dat);
        }
    glEnd();

    glPointSize(1);

    glFlush();
    SwapBuffers(hdc);
    }

Where float spread=M_PI; is global variable changed with mouse wheel ...