Fix CPU Particles spread

Fixes https://github.com/godotengine/godot/issues/51162
This is the same change as for the GPU Particles in https://github.com/godotengine/godot/pull/47310, including the recent fix for zero direction vectors.

scene/3d/cpu_particles.cpp

@@ -694,15 +694,28 @@ void CPUParticles::_particles_process(float p_delta) {
 				p.velocity = rot * parameters[PARAM_INITIAL_LINEAR_VELOCITY] * Math::lerp(1.0f, float(Math::randf()), randomness[PARAM_INITIAL_LINEAR_VELOCITY]);
 			} else {
 				//initiate velocity spread in 3D
-				float angle1_rad = Math::atan2(direction.x, direction.z) + (Math::randf() * 2.0 - 1.0) * Math_PI * spread / 180.0;
-				float angle2_rad = Math::atan2(direction.y, Math::abs(direction.z)) + (Math::randf() * 2.0 - 1.0) * (1.0 - flatness) * Math_PI * spread / 180.0;
+				float angle1_rad = (Math::randf() * 2.0 - 1.0) * Math_PI * spread / 180.0;
+				float angle2_rad = (Math::randf() * 2.0 - 1.0) * (1.0 - flatness) * Math_PI * spread / 180.0;
 
 				Vector3 direction_xz = Vector3(Math::sin(angle1_rad), 0, Math::cos(angle1_rad));
 				Vector3 direction_yz = Vector3(0, Math::sin(angle2_rad), Math::cos(angle2_rad));
-				direction_yz.z = direction_yz.z / MAX(0.0001, Math::sqrt(ABS(direction_yz.z))); //better uniform distribution
-				Vector3 direction = Vector3(direction_xz.x * direction_yz.z, direction_yz.y, direction_xz.z * direction_yz.z);
-				direction.normalize();
-				p.velocity = direction * parameters[PARAM_INITIAL_LINEAR_VELOCITY] * Math::lerp(1.0f, float(Math::randf()), randomness[PARAM_INITIAL_LINEAR_VELOCITY]);
+				Vector3 spread_direction = Vector3(direction_xz.x * direction_yz.z, direction_yz.y, direction_xz.z * direction_yz.z);
+				Vector3 direction_nrm = direction;
+				if (direction_nrm.length_squared() > 0) {
+					direction_nrm.normalize();
+				} else {
+					direction_nrm = Vector3(0, 0, 1);
+				}
+				// rotate spread to direction
+				Vector3 binormal = Vector3(0.0, 1.0, 0.0).cross(direction_nrm);
+				if (binormal.length_squared() < 0.00000001) {
+					// direction is parallel to Y. Choose Z as the binormal.
+					binormal = Vector3(0.0, 0.0, 1.0);
+				}
+				binormal.normalize();
+				Vector3 normal = binormal.cross(direction_nrm);
+				spread_direction = binormal * spread_direction.x + normal * spread_direction.y + direction_nrm * spread_direction.z;
+				p.velocity = spread_direction * parameters[PARAM_INITIAL_LINEAR_VELOCITY] * Math::lerp(1.0f, float(Math::randf()), randomness[PARAM_INITIAL_LINEAR_VELOCITY]);
 			}
 
 			float base_angle = (parameters[PARAM_ANGLE] + tex_angle) * Math::lerp(1.0f, p.angle_rand, randomness[PARAM_ANGLE]);