Calculate 3D unit vectors in a cone for MC-integration

cone_vector.h

struct vec2f {float x, y;};
struct vec3f {float x, y, z;};

//============================================================================
// cone_uniform_vector
//============================================================================
// Returns uniformly distributed unit vector on a [0, 0, 1] oriented cone of
// given apex angle and uniform random vector xi ([x, y] in range [0, 1]).
// e.g. cos_half_apex_angle = 0 returns samples on a hemisphere (cos(pi/2)=0),
// while cos_half_apex_angle = -1 returns samples on a sphere (cos(pi)=-1)
vec3f cone_uniform_vector(const vec2f &xi_, float cos_half_apex_angle_)
{
  float z=1.0f-(1.0f-cos_half_apex_angle_)*xi_.x;
  float r=sqrt(1.0f-z*z);
  float a=xi_.y*2.0f*3.1415926536f;
  float sphi=sin(a);
  float cphi=cos(a);
  vec3f res={r*cphi, r*sphi, z};
  return res;
}
//----------------------------------------------------------------------------


//============================================================================
// cone_strata_vector
//============================================================================
// Returns ith strata unit vector for a cone with given apex angle using
// "golden angle spiral method" described here: http://blog.marmakoide.org/?p=1
vec3f cone_strata_vector(unsigned sample_idx_, unsigned num_samples_, float cos_half_apex_angle_)
{
  float z=1.0f-(1.0f-cos_half_apex_angle_)*(sample_idx_+0.5f)/num_samples_;
  float r=sqrt(1.0f-z*z);
  float a=sample_idx_*(3.1415926536f*(3.0f-sqrt(5.0f)));
  float sphi=sin(a);
  float cphi=cos(a);
  vec3f res={r*cphi, r*sphi, z};
  return res;
}
//----------------------------------------------------------------------------


//============================================================================
// cone_solid_angle
//============================================================================
// Returns solid angle for a cone with given apex angle: https://en.wikipedia.org/wiki/Solid_angle
float cone_solid_angle(float cos_half_apex_angle_)
{
  return 2.0f*3.1415926536f*(1.0f-cos_half_apex_angle_);
}
//----------------------------------------------------------------------------