float sample_wavelength_from_absorption_triangle(const float *cf, const float rand, float *lambda)
{
  // cf[0]: emission peak wavelength in nanometer (nm)
  // cf[1]: concentration parameter c, not necessary for sampling the absorption spectrum
  // cf[2]: stokes shift in nm
  if (cf[0] == 0.f) return 0; // avoid sampling if emission peak cf[0] is invalid

  // compute absorption support and peak wavelength
  const float empeak_wn = 1/cf[0];
  const float abspeak_wn = 1/(cf[0] - cf[2]);
  const float abspeak_wl = 1.f / abspeak_wn; // peak wavelength
  const float stokes_wn = abspeak_wn - empeak_wn;
  const float abs_end_wn = abspeak_wn + 0.5*stokes_wn;
  const float abs_beg_wn = abspeak_wn - 0.5*stokes_wn;
  const float abs_end_wl = 1/abs_beg_wn; // end of absorption support
  const float abs_beg_wl = 1/abs_end_wn; // begin of absorption support

  assert(abs_end_wl > abs_beg_wl);

  // the support could also be used for sampling a b-spline (see paper)
  // This section details sampling a triangle function
  // with the given support and peak:

  const float a_p_minus_min = abspeak_wl - abs_beg_wl;
  const float a_support_width = abs_end_wl - abs_beg_wl;
  const float a_ratio_ascending = a_p_minus_min / a_support_width;

  if (rand < a_ratio_ascending)
  {
    // sample ascending triangle
    const float rand_scaled = rand / a_ratio_ascending;
    const float a2 = a_p_minus_min * a_p_minus_min;
    const float wl_offset = sqrt( a2 * rand_scaled);
    *lambda = abs_beg_wl + wl_offset;
    return (2.f * (wl_offset) / a2) * a_ratio_ascending;
  }
  else
  {
    // sample descending triangle
    const float rand_scaled = (rand - a_ratio_ascending) / (1.f - a_ratio_ascending);
    const float a = (abs_end_wl - abspeak_wl);
    const float wl_offset = a - a * sqrt(1.f - rand_scaled);
    *lambda = abspeak_wl + wl_offset;
    return (2.f / a - 2.f * wl_offset / (a*a) ) * (1.f - a_ratio_ascending); // return pdf
  }
}
