// ---------------------------------------------------------------------------
// Deterministic sampling of goniometric light sources using
// density stretching scheme in 2 dimensions. Works only for hemisphere.
// Avoids rejecting any input sample.
//
// ---------------------------------------------------------------------------
//
// Written by Vlastimil Havran, August, 2002.
// e-mail: havran@mpi-sb.mpg.de
//
// Max-Planck-Institut fur Informatik, Germany
// ---------------------------------------------------------------------------
// C++ standard headers

// standard headers
#include <iostream>
#include <cstdio>
//using namespace std;

#include "config.h"
#include "hemigon.h"

__BEGIN_BRDFLIB_SOURCE

// The constructor of the hemisphere gonio light
CHemisphereGonio::CHemisphereGonio()
{
  // the number of vertical angles (theta)
  VERTICAL_COUNT = -1;
  // the number of horizontal angles (phi)
  HORIZONTAL_COUNT = -1;

  // the angles in rads in vertical direction (theta)
  VERTICAL_ANGLES = 0;
  // the angles in rads in vertical direction (theta)
  HORIZONTAL_ANGLES = 0;

  // this is the description of the intensities in goniometric diagram
  INTENSITY = 0;

  // ------------------------------------------------------------------------
  // Auxiliary variables used by the algorithm:

  // sin of vertical angles, after remapping from 3D to 2D plane
  SIN_VERTICAL_ANGLES = 0;
  // square of the value above
  SIN_VERTICAL_ANGLES_SQR = 0;
  // the multiplier used to stretch the sample in 2D (radius)
  MULT = 0.0;
  // index of phi angle that has as a sum maximum energy
  MULT_index = -1;
  // precomputed array of values required to get fast search
  // in the patches changing the vertical angle
  GG = 0;
  // for rescaling in phi angle
  ALPHA = 0;
  // multiplier to stretch phi angle
  _maxPhiRescale = 0.0;

  return;
}

void
CHemisphereGonio::CleanUp()
{
  delete [] VERTICAL_ANGLES;
  VERTICAL_ANGLES = 0;

  delete [] HORIZONTAL_ANGLES;
  HORIZONTAL_ANGLES = 0;

  delete [] INTENSITY;
  INTENSITY = 0;

  delete [] SIN_VERTICAL_ANGLES;
  SIN_VERTICAL_ANGLES = 0;

  delete [] SIN_VERTICAL_ANGLES_SQR;
  SIN_VERTICAL_ANGLES_SQR = 0;

  delete [] GG;
  GG = 0;

  delete [] ALPHA;
  ALPHA = 0;

  return;
}


// Allocate arrayes to represent hemisphere intensities
// and also auxiliary variables.
void
CHemisphereGonio::AllocArrays()
{
  assert(VERTICAL_COUNT > 1);
  assert(HORIZONTAL_COUNT > 1);

  // input values;
  VERTICAL_ANGLES = new double[VERTICAL_COUNT];
  assert(VERTICAL_ANGLES);

  HORIZONTAL_ANGLES = new double[HORIZONTAL_COUNT+1];
  assert(HORIZONTAL_ANGLES);
  HORIZONTAL_ANGLES[HORIZONTAL_COUNT] = M_PI * 2.0;

  INTENSITY = new double[VERTICAL_COUNT * HORIZONTAL_COUNT];
  assert(INTENSITY);

  // precomputed values
  SIN_VERTICAL_ANGLES = new double[VERTICAL_COUNT];
  assert(SIN_VERTICAL_ANGLES);

  SIN_VERTICAL_ANGLES_SQR = new double[VERTICAL_COUNT];
  assert(SIN_VERTICAL_ANGLES_SQR);

  GG = new double[VERTICAL_COUNT * HORIZONTAL_COUNT];
  assert(GG);

  ALPHA = new double[HORIZONTAL_COUNT+1];
  assert(ALPHA);

  return;
}

// Initialization of the intensities in keypoints of goniometric web
// Returns 'false' on success, 'true' on failure.
bool
CHemisphereGonio::InitIntensities(
  int phi_num, double *phi,
  int theta_num, double *theta,
  double *diagram)
{
  // the number of vertical angles (theta)
  VERTICAL_COUNT = theta_num;
  assert(VERTICAL_COUNT > 1);

  // the number of horizontal angles (phi)
  HORIZONTAL_COUNT = phi_num;
  assert(HORIZONTAL_COUNT > 1);

  // initialize all arrays
  AllocArrays();

  // loop variables
  int i, j;

  for (i = 0; i < VERTICAL_COUNT; ++i)
    VERTICAL_ANGLES[i] = theta[i];

  // the theta[0] must specify the pole of the hemisphere
  if (VERTICAL_ANGLES[0] != 0)
    return true; // error;

  // setting phi angles
  for (i = 0; i < HORIZONTAL_COUNT; ++i)
    HORIZONTAL_ANGLES[i] = phi[i];
  // This below is not completely correct!
  if (HORIZONTAL_ANGLES[0] < 0)
    HORIZONTAL_ANGLES[0] = 0;
  // set the last angle correctly
  HORIZONTAL_ANGLES[HORIZONTAL_COUNT] =
     HORIZONTAL_ANGLES[0] + 2.0 * M_PI;

  // Two setting below are equal, they just differ
  // by indexing organization
  for (i = 0; i < VERTICAL_COUNT; ++i) { // theta
    for (j = 0; j < HORIZONTAL_COUNT; ++j) { // phi
      INTENSITY[GetIndex(i, j)] = diagram[i + theta_num * j];
      //cout << "Int[" << i << "][" << j << "] = Int(" << GetIndex(i, j) << ") = " << INTENSITY[GetIndex(i, j)] << "\n";
    } // for j
  } // for i

  // initialization of auxiliary variables for bands
  InitTheta();

  // initialization of auxiliary variables for various phi angles.
  InitPhi();

  // Compute the total power given by the goniometric web description
  // as the integral over the intensities.
  ComputeTotalPower();

  // success,, all variables initialized
  return false;
}

// --------------------------------------------------------------------------
// Initialization of variables required for sampling in theta direction
void
CHemisphereGonio::InitTheta()
{
  int i;

#ifdef _DEBUG
  for (i = 0; i < HORIZONTAL_COUNT; i++)
  {
    assert(HORIZONTAL_ANGLES[i] >= 0.0);
    assert(HORIZONTAL_ANGLES[i] <= 2.0 * M_PI);
    if (i > 0) {
      assert(HORIZONTAL_ANGLES[i] > HORIZONTAL_ANGLES[i-1]);
    }
  }
#endif // _DEBUG

  // for each vertical slize of hemisphere, increasing theta gradually
  for (i = 0; i < VERTICAL_COUNT; i++)
  {
    assert(VERTICAL_ANGLES[i] >= 0.0);
    assert(VERTICAL_ANGLES[i] <= M_PI_2);
    if (i > 0) {
      assert(VERTICAL_ANGLES[i] > VERTICAL_ANGLES[i-1]);
    }

    // compute the sinus of the angle that is used actually!

    // The formula below takes the 2D->3D mapping into account
    // so it preserves the angles specified in 3D space!
    // The 2D->3D mapping must be then: theta = acos(1-r*r), where
    // 'r' is the radius on 2D disk.
    SIN_VERTICAL_ANGLES_SQR[i] = 1.0 - cos(VERTICAL_ANGLES[i]);
    SIN_VERTICAL_ANGLES[i] = sqrt(SIN_VERTICAL_ANGLES_SQR[i]);
  } // for i

  // find out the maximum when integrating the energy for all bands!
  MULT = 0.0;
  int phiIndex;
  // for each Horizontal slice of hemisphere, increasing phi gradually
  for (phiIndex = 0; phiIndex < HORIZONTAL_COUNT; phiIndex++) {
    int thetaIndex;
    double multMax = 0.0;
    double multMaxSqr = 0.0;
    // for each angle theta, increasing theta gradually
    for (thetaIndex = 0; thetaIndex < VERTICAL_COUNT; thetaIndex++) {
      // compute the index of relative intensity
      int index = GetIndex(thetaIndex, phiIndex);
      if (INTENSITY[index] < 0.0) {
	cout << "The relative intensity for general goniometric lightsource"
	      << " must be always greater than zero!" << endl;
	abort();
      }

      // create the mapping entry for this segment
      if (thetaIndex == 0) {
	assert(VERTICAL_ANGLES[thetaIndex] == 0.0);
	GG[index] = 0.0; // for initial angle=0.0 it is always so
      }
      else
      {
	int index2 = GetIndex(thetaIndex-1, phiIndex);
 	// perform remmapping because of further 2D->3D remapping
	// work on sin of these angles (projection to 2D plane)

	// the intensities corresponding to the angles
	double intensityMaxAngle = INTENSITY[index];
	double intensityMinAngle = INTENSITY[index2];
	
	if (intensityMinAngle == intensityMaxAngle) {
	  // constant linear mapping
	  double g = multMaxSqr + (SIN_VERTICAL_ANGLES_SQR[thetaIndex] -
				  SIN_VERTICAL_ANGLES_SQR[thetaIndex-1]) * intensityMinAngle;
	  GG[index] = g;
	  multMaxSqr = g;
	  multMax = sqrt(g);
	}
	else {
	  // arbitrary linear mapping, but not constant!
	  double a = SIN_VERTICAL_ANGLES[thetaIndex-1];
	  double aa = SIN_VERTICAL_ANGLES_SQR[thetaIndex-1]; // or sqr(a)
	  double c = SIN_VERTICAL_ANGLES[thetaIndex];
	  double cc = SIN_VERTICAL_ANGLES_SQR[thetaIndex]; // or sqr(c)
	  double tt = (intensityMinAngle - intensityMaxAngle)/(c-a);
	  double g = (c*(1.0/3.0 * cc - aa) + aa * a * 2.0/3.0) * tt +
	             multMaxSqr + intensityMaxAngle * (cc - aa);
	  assert(g >= multMaxSqr);
	  GG[index] = g;
	  multMaxSqr = g;
	  multMax = sqrt(g);
	}
#if 0
	cout << "ThetaIndex=" << thetaIndex
	     << " phiIndex=" << phiIndex << " index = " << index
	     << " ANGLEorigMin = " << VERTICAL_ANGLES[thetaIndex-1]
	     << " ANGLEorigMax = " << VERTICAL_ANGLES[thetaIndex]
	     << " f = " << GG[index2] << " g = " << GG[index]
	     << endl;
#endif
      } // if
    } // for thetaIndex

#if 0
    cout << "phiIndex = " << phiIndex << " multMax = " << multMax << endl;
#endif
    // update the maximum
    if (multMax > MULT) {
      MULT = multMax;
      MULT_index = phiIndex;
    }
  } // for phiIndex

#if 1
  // to guarantee the same stretching for all bands
  cout << "MULT = " << setprecision(10) << MULT << endl;
#endif

  return;
}

// --------------------------------------------------------------------------
// Initialization of variables required for sampling in phi direction
void
CHemisphereGonio::InitPhi()
{
  int phiIndex;

  // cumulative sum of angles
  ALPHA[0] = 0.0;
  
  // for each Horizontal slice of hemisphere, increasing phi gradually
  for (phiIndex = 0; phiIndex < HORIZONTAL_COUNT; phiIndex++) {
    double D0 = GG[GetIndex(VERTICAL_COUNT-1, phiIndex % HORIZONTAL_COUNT)];
    double D1 = GG[GetIndex(VERTICAL_COUNT-1, (phiIndex+1) % HORIZONTAL_COUNT)];
    // base angle of the result phi angle
    double alpha = ALPHA[phiIndex % HORIZONTAL_COUNT ];
    // increase in phi angle
    // double beta = M_PI * 2.0 / (double)HORIZONTAL_COUNT;

    // angle between two vertices of goniometric web in in phi-direction
    double PHI_SEGMENT =
      HORIZONTAL_ANGLES[phiIndex+1] - HORIZONTAL_ANGLES[phiIndex];

    // now compute gamma for this particular range of phi
    ALPHA[phiIndex + 1] = alpha + PHI_SEGMENT / sqr(MULT) * (D1 + D0) / 2.0;

#if 0
    cout << "phiIndex = " << phiIndex << " D0 = " << D0 << " D1 = " << D1
	 << " ALPHA[phiIndex]=" << ALPHA[phiIndex + 1]
	 << " origalpha[phiIndex] = "
	 << (double)(phiIndex + 1) * M_PI * 2.0 / HORIZONTAL_COUNT
	 << endl;
#endif
  }
  cout << "----------" << endl;


  // phi is rescaled by this parameter
  _maxPhiRescale = ALPHA[ HORIZONTAL_COUNT] / (2.0 * M_PI);
  cout << "_maxPhiRescale = " << _maxPhiRescale << endl;

  return;
}

// ---------------------------------------------------------
// Max
//   Returns maximum of two values
// ---------------------------------------------------------
/*template<class T>
inline const T&
Max(const T &x, const T &y)
{
  return (x < y) ? y : x;
}*/

// -----------------------------------------------------------------------
// Computes total power (flux) over the hemisphere using the intesities
// description via the goniometric web
void
CHemisphereGonio::ComputeTotalPower()
{
  int iphi, itheta, ii;
  double a1, a2, b1, b2, dp, dt, dt0, dt1, aa, bb;

  // the information to be computed is initialized
  maximumIntensity = -1.0;
  totalPower = 0.0;
  
  for (iphi = 0; iphi < HORIZONTAL_COUNT - 1; ++iphi) {
    for (itheta = 0; itheta < VERTICAL_COUNT - 1; ++itheta) {
      // the four values for given goniometric web diagram patch
      a1 = INTENSITY[ii = GetIndex(itheta, iphi)]; // INTENSITY[ii = iphi * theta_num + itheta];
      b1 = INTENSITY[GetIndex(itheta, iphi+1)];    // INTENSITY[ii + 1];
      a2 = INTENSITY[GetIndex(itheta+1, iphi)];    // diagram[ii += theta_num];
      b2 = INTENSITY[GetIndex(itheta+1, iphi+1)];  // diagram[ii + 1];

      if ((dt = Max(Max(a1, a2), Max(b1, b2))) > maximumIntensity)
        maximumIntensity = dt;

      dt = (dt1 = VERTICAL_ANGLES [itheta + 1]) - (dt0 = VERTICAL_ANGLES[itheta]);
      dp = HORIZONTAL_ANGLES[iphi + 1] - HORIZONTAL_ANGLES[iphi];

      // A * dp * 0.5
      aa = 0.5 * dp / dt * (-sin(dt1) + sin(dt0) + dt * cos(dt0));
      // B * dp * 0.5
      bb = 0.5 * dp / dt * ( sin(dt1) - sin(dt0) - dt * cos(dt1));
      //aa = 0.25; bb = 0.25;

      totalPower += (a1 + a2) * aa + (b1 + b2) * bb;
    } // for itheta
  } // for iphi

  return;
}

// -------------------------------------------------------------------------
// Solve a special case of equation: R^3 + R^2 * B + C = 0
// Returns always real root found (root is always >= 0.0)
double
CHemisphereGonio::SolveCubicRoot(bool ed, double B, double C)
{
  double B3 = B * B * B;
  double result;
  
  if (ed) { // e < d
    // evaluation of the first root using Cardan formulas
    double V = C * (81.0 * C + 12.0 * B3);

    // we must use complex numbers, however as a result we need only
    // real component
    complex<double> D(V, 0);
    D = sqrt(D);
    complex<double> T1 = -108.0 * C - 8.0 * B3 + 12.0 * D;      
    T1 = pow(T1, (double)(1.0/3.0));
    complex<double> W =
      -1.0/12.0*T1*(complex<double>(1.0, sqrt(3.0)))
      -1.0/3.0*B*(1.0 + B/T1*(complex<double>(1.0, - sqrt(3.0) )));

#ifdef _DEBUG    
    // note W.imag() == 0!
    if (fabs(W.imag()) >= 1e-3) {
      cout << "W.imag() is not zero = " << W.imag() << endl;
    }
#endif // _DEBUG
    result = W.real();
  }
  else { // e > d
    // evaluation the second root using Cardan formulas
    double V = C * (81.0 * C + 12.0 * B3);

    if (V > 0.0) {
      V = sqrt(V);
      double T2 = -108.0 * C - 8.0 * B3 + 12.0 * V;
      if (T2 > 0.0) {
	double T1 = pow(T2, 1.0/3.0);
	result = 1.0/6.0 * T1 + 2.0/3.0 * B * B/T1 - 1.0/3.0 * B;
      }
      else { // T2 < 0.0
	complex<double> T1 = pow(T2, 1.0/3.0);
	complex<double> W = 1.0/6.0 * T1 + 2.0/3.0 * B * B/T1 - 1.0/3.0 * B;
	result = W.real();
#ifdef _DEBUG
	// note W.imag() == 0! for real root
	if (fabs(W.imag()) >= 1e-3) {
	  cout << "W.imag() is not zero = " << W.imag() << endl;
	}
#endif // _DEBUG
      }
    }
    else { // V < 0.0
      complex<double> D(V, 0);
      D = sqrt(D);
      complex<double> T1 = -108.0 * C - 8.0 * B3 + 12.0 * D;
      T1 = pow(T1, (double)(1.0/3.0));
      complex<double> W = 1.0/6.0 * T1 + 2.0/3.0 * B * B/T1 - 1.0/3.0 * B;
      result = W.real();
#ifdef _DEBUG    
      // note W.imag() == 0! for real root
      if (fabs(W.imag()) >= 1e-3) {
	cout << "W.imag() is not zero = " << W.imag() << endl;
      }
#endif // _DEBUG
    }
  }
  return result;
}

// -------------------------------------------------------------------------------
// Searches a digon with remapping of 'phi' variable
void
CHemisphereGonio::LocateDigonStretch(double &phi, double &radius,
				     int &left, int &right,
				     double &alpha, double &effectiveDensity)
{
  // remap phi angle (phi)
  phi *= _maxPhiRescale;
  // locate the slice for which fits phi angle
  right = -1;
  // how 's' is remapped
  effectiveDensity = 1.0;

  // !!!!!!!!!!!!
  // !!! Note that this linear search can be changed to binary search,
  // !!! thus, the complexity from O(N) to O(log N).
  // !!!

#define USE_BINARY_SEARCH_PHI
  
#ifdef USE_BINARY_SEARCH_PHI

  if (phi > ALPHA[HORIZONTAL_COUNT]) {
    cerr << "LocateDigonStretch: FATAL ERROR FOR phi = " << phi << endl;
    exit(3);
  }
  else {
    // use binary search to locate correct band
    int li = 0, ri = HORIZONTAL_COUNT-1;
    for (;;) {
      int b = (li + ri)/2;
      if (phi < ALPHA[b]) { // search in the left half
	ri = b - 1;
      }
      else {
	if (phi > ALPHA[b+1]) { // search in the right half
	  li = b + 1;
	}
	else {// (rr >= ff) && (rr <= gg)
	  left = b;
	  right = left+1;
	  break; // found correct band
	}
      }
    } // for
  }
#else
  // linear search
  for (left = 0; left < HORIZONTAL_COUNT; left++)
    if ( (ALPHA[left] <= phi) &&
	 (phi <= ALPHA[(left+1) ]) ) {
      right = (left + 1);
      break; // found correct index for phi angle
    }
#endif

  assert( (phi > ALPHA[left]) && (phi < ALPHA[left+1]) );

  if (right == -1) {
    cerr << "SampleOne: FATAL ERROR FOR phi = " << phi << endl;
    exit(3);
  }

  // for purposes of visualization
  last_sector = left;
  
  // the interpolation coeficient in the horizontal direction (phi)   
  // between left and right phi-angles.
  // It holds for alpha:
  //   'left' contribution is multiplied by 'alpha', 'right' by '1.0-alpha'
  // interpolated value is then: V = VLeft * alpha + VRight * (1-alpha)
  double GG_LEFT = GG[GetIndex(VERTICAL_COUNT - 1, left)];
  double GG_RIGHT = GG[GetIndex(VERTICAL_COUNT - 1, (left+1) % HORIZONTAL_COUNT)];

  // what is the linear interpolation coefficient in the input domain.
  double A = ((double)phi - ALPHA[left])/(ALPHA[right] - ALPHA[left]);

  if (fabs(GG_LEFT - GG_RIGHT) < 1e-7) {
    // simple case, linear interpolation of alpha
    alpha = A;
    effectiveDensity = sqr(MULT) / GG_LEFT;
  }
  else {
    // common case, different densities

    double a, b, c, det, GammaA, alpha2;
    
    GammaA = phi - ALPHA[left];
    assert(GammaA >= 0.0);

    double PHI_SEGMENT =
      HORIZONTAL_ANGLES[left+1] - HORIZONTAL_ANGLES[left];
    
    // solve the quadratic equation
    a = PHI_SEGMENT / 2.0 * (GG_RIGHT - GG_LEFT);
    b = PHI_SEGMENT * GG_LEFT;
    c = -GammaA * sqr(MULT);

    det = b * b - 4.0 * a * c;
    if (det < 0.0) {
      cerr << "det = " << det << endl;
      cerr << "a = " << a << " b = " << b << " c = " << c << endl;
      cerr << "Determinant cannot be principially negative - implem. bug" << endl;
      exit(3);
    }
    det = sqrt(det);
    // compute the first root of the quadratic equation
    alpha = (-b + det)/(2.0 * a);
    
    // cout << "alpha = " << alpha << endl;
#if 0
    if ((alpha <= 0.0) || (alpha > 1.0)) {
      // compute the second root
      cout << "AAA" << endl;
      alpha = (-b - det)/(2.0 * a);
      if (GG_RIGHT > GG_LEFT)
	alpha = (-b + det)/(2.0 * a);
      else
	alpha = (-b - det)/(2.0 * a);
    }
#endif
#if 1
    alpha2 = (-b - det)/(2.0 * a);
    if ((alpha2 >= 0.0) && (alpha2 <= 1.0)) {
      cout << "alpha2 = " << alpha2 << " alpha = " << alpha << endl;
    }
#endif
    // compute the effective density
    effectiveDensity = sqr(MULT)/( (1.0 - alpha) * GG_LEFT + alpha * GG_RIGHT);

#if 0    
    cout << " phi= " << phi << " a= " <<a<< " b= " << b << " c= " << c
	 << " det= " << det << " alpha= " << alpha
	 << " ed= " << effectiveDensity
	 << " alpha2= " << alpha2 << endl;
#endif    
  } // -------------

#ifdef _DEBUG
  if ((alpha < 0.0) || (alpha > 1.0)) {
    cout << "Alpha = " << alpha << " left = " << left << endl;
    cout << " radius = " << radius << " phi = " << phi << endl;
    cout << "FATAL ERROR" << endl;
    exit(3);
  }

  const float epsilon = 1e-15;
  if (effectiveDensity < 1.0 - epsilon) {
    cout << "effectiveDensity = " << effectiveDensity << " left = " << left << endl;
    cout << " radius = " << radius << " phi = " << phi << endl;
    cout << "FATAL ERROR" << endl;
    exit(3);
  }  
#endif // _DEBUG

  // possibly, compute correct value of right index, that is further used
  // to access GG array.
  right = right % HORIZONTAL_COUNT;

  // compute new value of phi
  phi = HORIZONTAL_ANGLES[left];
  phi += alpha * (HORIZONTAL_ANGLES[left+1] - HORIZONTAL_ANGLES[left]);

  // modify also radius by effective density for given 'A' and thus 'alpha'
  radius /= sqrt(effectiveDensity);
  
  return;
}

// -------------------------------------------------------------------------------
// Performs that density stretching to fit the description of goniometric
// light source:
// modifies 'radius' and 'phi' according to the goniometric light source diagram
// Returns false, when the sample is valid.
// Returns true, when the sample cannot be taken as the input of
// the goniometric light source diagram -- ERROR for current implementation.
//
// Assumes radius to be in the range <0.0-1.0>  ... radius in 2D disk
//         phi to be in the range <0.0-2*PI> ... phi angle
bool
CHemisphereGonio::SampleOne(double &radius, double &phi, int index)
{
  index = index; // satisfy the compiler  

  // we have to first multiply radius by a precomputed constant
  radius *= MULT;
  
  // the indeces of the digon
  int left, right;
  double alpha;
  double effectiveDensity = 1.0;

  // locate simple digon and compute effective density, avoiding
  // rejection sampling completely.
  LocateDigonStretch(phi, radius, left, right, alpha, effectiveDensity);

  // the bands to be used
  int b;
  double rr = sqr(radius);
  double ff = 0.0, gg;
  // the indeces
  int indexMinRight, indexMaxRight;
  int indexMinLeft, indexMaxLeft;

#define USE_BINARY_SEARCH_THETA
  
#ifdef USE_BINARY_SEARCH_THETA
  // More sophisticated and faster approach, especially for higher
  // number of angles subdividing the theta angle range.
  //
  // first test if the sample is not out of range at all
  indexMaxLeft  = GetIndex(VERTICAL_COUNT - 1, left);
  indexMaxRight = GetIndex(VERTICAL_COUNT - 1, right);
  gg = GG[indexMaxLeft] * (1.0 - alpha) + GG[indexMaxRight] * alpha;

  if (rr > gg) {
    // handling numerical inacuracy
    const double epsilon = 1e-15;
    if (rr > gg + epsilon) {
      cout << "ERROR - NOT FOUND AA " << " rr = " << rr << " gg(interpolated) = " << gg;
      cout << " gg_left = " << GG[indexMaxLeft]
	   << " gg_right = " << GG[indexMaxRight] << endl;
      return true; // it has no solution, mapped outside possible range
    }
    // the last band of cells (closest to equator, theta = M_PI/2)
    b = VERTICAL_COUNT-2;
  }
  else {
    // use binary search to locate correct band
    int li = 0, ri = VERTICAL_COUNT-2;
    for (;;) {
      b = (li + ri)/2;
      indexMinLeft = GetIndex(b, left);
      indexMinRight = GetIndex(b, right);
      ff = GG[indexMinLeft] * (1.0 - alpha) + GG[indexMinRight] * alpha;
      if (rr < ff) { // search in the left half
	ri = b - 1;
      }
      else {
	indexMaxLeft  = GetIndex(b+1, left);
	indexMaxRight = GetIndex(b+1, right);
	gg = GG[indexMaxLeft] * (1.0 - alpha) + GG[indexMaxRight] * alpha;
	if (rr > gg) { // search in the right half
	  li = b + 1;
	}
	else // (rr >= ff) && (rr <= gg)
	  break; // found correct band
      }
    } // for
  }
#else
  // simpler approach .. linear search
  for (b = 1; b < VERTICAL_COUNT; b++) {
    int lindex = GetIndex(b, left);
    int rindex = GetIndex(b, right);

    gg = GG[lindex] * (1.0 - alpha) + GG[rindex] * alpha;
    if ( (rr >= ff) && (rr <= gg)) {
      b--;
      break; // found correct band
    }
    // go to the next band
    ff = gg;
  } // for b
  
  if (b == VERTICAL_COUNT) { // not found
    // handling numerical inacuracy
    const double epsilon = 1e-15;
    if (rr > gg + epsilon) {
      cout << "ERROR - NOT FOUND A " << "f = " << sqrt(ff) << "   gg = " << gg
	   << " radius = " << radius << endl;
      return true; // invalid sample
    }
  }
#endif // binary search

  // for purposes of visualization
  last_band = b;

  // handling numerical inaccuracy
  const float epsilon = 1e-15;
  if (rr < ff + epsilon) {
    radius = SIN_VERTICAL_ANGLES[b];
    return false;
  }
  if (rr > gg - epsilon) {
    radius = SIN_VERTICAL_ANGLES[b+1];
    return false;
  }
  
  // cout << "rr - ff = " << rr - ff << " gg - rr = " << gg - rr << endl;
  
  // now the right patch is found, perform the stretching itself
  
  // now interpolate, since it is always possible
  // d and e is required density/intensity
  double d = INTENSITY[GetIndex(b, left)] * (1.0 - alpha) +
             INTENSITY[GetIndex(b, right)] * alpha;
  double e = INTENSITY[GetIndex(b+1, left)] * (1.0 - alpha) +
             INTENSITY[GetIndex(b+1, right)] * alpha;

  // constant density/intensity in the patch for given phi ?
  if (fabs(d-e) < 1e-12) {
    // constant remapping ... easy case
    if (fabs(d) > 1e-12)
      radius = sqrt(SIN_VERTICAL_ANGLES_SQR[b] + (rr - ff)/d);
    else
      radius = SIN_VERTICAL_ANGLES[b]; // quite unlikely

#ifdef _DEBUG
    double a = SIN_VERTICAL_ANGLES[b];
    double c = SIN_VERTICAL_ANGLES[b+1];
    if (radius < a) {
      cout << "error A radius = " << radius << " a = " << a << " c = " << c << " band = " << b
	   << " lmin = " << INTENSITY[GetIndex(b, left)]
	   << " rmin = " << INTENSITY[GetIndex(b, right)]
	   << " lmax = " << INTENSITY[GetIndex(b+1, left)]
	   << " rmax = " << INTENSITY[GetIndex(b+1, right)]
	   << " d = " << d << " e = " << e << " alpha = " << alpha << endl;
      return true; // error
    }

    if (radius > c) {
      cout << "error B radius = " << radius << " a = " << a << " c = " << c << " band = " << b
	   << " lmin = " << INTENSITY[GetIndex(b, left)]
	   << " rmin = " << INTENSITY[GetIndex(b, right)]
	   << " lmax = " << INTENSITY[GetIndex(b+1, left)]
	   << " rmax = " << INTENSITY[GetIndex(b+1, right)]
	   << " d = " << d << " e = " << e << " alpha = " << alpha << endl;
      return true; // error ... sample is not usable
    }
#endif // _DEBUG
  }
  else { // different densities/intensities
    // linear remapping, it corresponds to solving a cubic equation
    double a = SIN_VERTICAL_ANGLES[b];
    double c = SIN_VERTICAL_ANGLES[b+1];
    double aa = SIN_VERTICAL_ANGLES_SQR[b]; // or sqr(a);
    double tt = (c - a)/(d - e);

    // Prepare input parameters for the cubic equation
    double B = (-3.0 * (e * tt + c))/2.0;
    double C = 3.0/2.0 * ((e * aa + (rr - ff)) * tt + c * aa) - aa * a;    
    bool ed = (e < d) ? true : false;
    // Solve the cubic equation
    radius = SolveCubicRoot(ed, B, C);

#ifdef _DEBUG
    if (radius < a) {
      cout << "error C radius = " << radius << " a = " << a << " c = " << c
	   << " band = " << b
	   << " lmin = " << INTENSITY[GetIndex(b, left)]
	   << " rmin = " << INTENSITY[GetIndex(b, right)]
	   << " lmax = " << INTENSITY[GetIndex(b+1, left)]
	   << " rmax = " << INTENSITY[GetIndex(b+1, right)]
	   << " d = " << d << " e = " << e << " alpha = " << alpha << endl;
      return true; // error
    }

    // handling numerical inaccuracy
    const float epsilon = 1e-15;
    if (radius > c + epsilon) {
      // this can rarely happen (inaccuracy)
      cout << " index = " << index << endl;
      cout << "error D radius = " << radius << " a = " << a << " c = " << c
	   << " band = " << b
	   << " lmin = " << INTENSITY[GetIndex(b, left)]
	   << " rmin = " << INTENSITY[GetIndex(b, right)]
	   << " lmax = " << INTENSITY[GetIndex(b+1, left)]
	   << " rmax = " << INTENSITY[GetIndex(b+1, right)]
	   << " d = " << d << " e = " << e << " alpha = " << alpha << endl;
      return true; // error ... sample is not usable
    }
#endif // _DEBUG
  } // different densities/intensities

#ifdef _DEBUG
  if (radius > 1.0) {
    cout << "error E" << " radius = " << radius << " phi = " << phi
	 << " effDensity = " << effectiveDensity << endl;
    // this sample cannot be used, some bug or numerical inaccuracy
    // It should never happen!
    return true; // error
  }
#endif  

  return false; // valid sample
}

// -------------------------------------------------------------------------------
// Searches a digon with remapping of 'phi' variable for inverse setting (going
// from the hemisphere direction to the 2D rectangle)
void
CHemisphereGonio::LocateDigonStretchInverse(
  double &phi, double & /*radius*/,
  int &left, int &right,
  double &A, double &effectiveDensity)
{
  // locate the slice for which fits phi angle
  right = -1;
  // how 's' is remapped
  effectiveDensity = 1.0;

  // Get the phi into the correct range
  while (phi > 2.0 * M_PI) {
    phi -= 2.0 * M_PI;
  }
  while (phi < 0.0) {
    phi += 2.0 * M_PI;
  }
  
#define USE_BINARY_SEARCH_PHI
  
#ifdef USE_BINARY_SEARCH_PHI

  // use binary search to locate correct band
  int li = 0, ri = HORIZONTAL_COUNT-1;
  for (;;) {
    int b = (li + ri)/2;
    if (phi < HORIZONTAL_ANGLES[b]) { // search in the left half
      ri = b - 1;
    }
    else {
      if (phi > HORIZONTAL_ANGLES[b+1]) { // search in the right half
	li = b + 1;
      }
      else {// (rr >= ff) && (rr <= gg)
	left = b;
	right = left+1;
	break; // found correct band
      }
    }
  } // for
#else
  // linear search
  for (left = 0; left < HORIZONTAL_COUNT; left++)
    if ( (HORIZONTAL_ANGLES[left] <= phi) &&
	 (phi <= HORIZONTAL_ANGLES[(left+1) ]) ) {
      right = (left + 1);
      break; // found correct index for phi angle
    }
#endif

  assert(phi >= HORIZONTAL_ANGLES[left]);
  assert(phi <= HORIZONTAL_ANGLES[left+1]);

  if (right == -1) {
    cerr << "SampleOne: FATAL ERROR FOR phi = " << phi << endl;
    exit(3);
  }

  // for purposes of visualization
  last_sector = left;
  
  // the interpolation coeficient in the horizontal direction (phi)   
  // between left and right phi-angles.
  // It holds for alpha:
  //   'left' contribution is multiplied by 'alpha', 'right' by '1.0-alpha'
  // interpolated value is then: V = VLeft * alpha + VRight * (1-alpha)
  double GG_LEFT = GG[GetIndex(VERTICAL_COUNT - 1, left)];
  double GG_RIGHT = GG[GetIndex(VERTICAL_COUNT - 1, (left+1) % HORIZONTAL_COUNT)];

  // what is the linear interpolation coefficient in the input domain.
  double alpha = ((double)phi - HORIZONTAL_ANGLES[left]) /
                  (HORIZONTAL_ANGLES[right] - HORIZONTAL_ANGLES[left]);
  
  if (fabs(GG_LEFT - GG_RIGHT) < 1e-7) {
    // simple case, linear interpolation of alpha
    A = alpha;
    effectiveDensity = sqr(MULT) / GG_LEFT;
  }
  else {
    // common case, different densities
    // compute the effective density
    effectiveDensity = sqr(MULT)/( (1.0 - alpha) * GG_LEFT + alpha * GG_RIGHT);
    double PHI_SEGMENT =
      HORIZONTAL_ANGLES[left+1] - HORIZONTAL_ANGLES[left];
    assert(PHI_SEGMENT >= 0.0);
    
    double a, b;

    // solve the quadratic equation
    a = PHI_SEGMENT / 2.0 * (GG_RIGHT - GG_LEFT);
    b = PHI_SEGMENT * GG_LEFT;

    // compute A from the alpha
    double gammaA = 1.0/sqr(MULT) * (sqr(alpha)*a + alpha*b);

    A = gammaA / (ALPHA[right] - ALPHA[left]);
    assert(A >= 0.0);
    assert(A <= 1.0);
  } // -------------

  // possibly, compute correct value of right index, that is further used
  // to access GG array.
  right = right % HORIZONTAL_COUNT;

  // compute new value of phi
  phi = ALPHA[left];
  phi += A * (ALPHA[left+1] - ALPHA[left]);
  phi /= _maxPhiRescale;

  // for the output, we actually need alpha
  A = alpha;
  
  // now 'radius' and 'phi' should be correct - with respect to the forward
  // mapping just the opposite.
  
  return;
}

// The inverse function to the one above:
//
// Assumes input:
//         radius to be in the range <0.0-1.0>  ... radius in 2D disk
//         phi to be in the range <0.0-2*PI> ... phi angle
bool
CHemisphereGonio::SampleOneInverse(double &radius, double &phi, int index)
{
  index = index; // satisfy the compiler  

  // the indeces of the digon
  int left, right;
  double alpha;
  double effectiveDensity = 1.0;

  // locate simple digon and compute effective density, avoiding
  // rejection sampling completely.
  LocateDigonStretchInverse(phi, radius, left, right, alpha, effectiveDensity);
  
  // now phi is correctly remapped, we have to remap the radius
  // The radius input is now in the range <0.0-GG_LEFT> or <0.0-GG_RIGHT>
  // or interpolated by alpha.
  
  // the band to be searched for from the angles
  int b;
  double sinMin = 0., sinMax = 0.;
  
#ifdef USE_BINARY_SEARCH_THETA
  // More sophisticated and faster approach, especially for higher
  // number of angles subdividing the theta angle range.

  // the sin of maximum angle
  sinMax = SIN_VERTICAL_ANGLES[VERTICAL_COUNT - 1];
  
  if (radius > sinMax) {
    // handling numerical inacuracy
    const double epsilon = 1e-15;
    if (radius > sinMax + epsilon) {
      cout << "ERROR - NOT FOUND AAINVERSE rardius = " << radius << endl;
      return true; // it has no solution, mapped outside possible range
    }
    // the last band of cells (closest to equator, theta = M_PI/2)
    b = VERTICAL_COUNT-2;
  }
  else {
    // use binary search to locate correct band
    int li = 0, ri = VERTICAL_COUNT-2;
    for (;;) {
      b = (li + ri)/2;
      sinMin = SIN_VERTICAL_ANGLES[b];
      sinMax = SIN_VERTICAL_ANGLES[b+1];
      if (radius < sinMin) { // search in the left half
	ri = b - 1;
      }
      else {
	if (radius > sinMax) { // search in the right half
	  li = b + 1;
	}
	else {// (rr >= ff) && (rr <= gg)
	  break; // found correct band
	}
      }
    } // for
  }
  
#else
  // simpler approach .. linear search
  sinMin = IN_VERTICAL_ANGLES[0];
  for (b = 1; b < VERTICAL_COUNT; b++) {
    double sinMax = SIN_VERTICAL_ANGLES[b];
    if ( (rr >= sinMin) && (rr <= sinMax)) {
      b--;
      break; // found correct band
    }
    // go to the next band
    sinMin = sinMax;
  } // for b
  
  if (b == VERTICAL_COUNT) { // not found
    // handling numerical inacuracy
    const double epsilon = 1e-15;
    if (radius > sinMax + epsilon) {
      cout << "ERROR - NOT FOUND A - radius = " << radius << endl;
      return true; // invalid sample
    }
  }
#endif // binary search

  // for purposes of visualization
  last_band = b;
  
  // the values required to be computed
  int indexMaxLeft  = GetIndex(b+1, left);
  int indexMinLeft  = GetIndex(b, left);
  int indexMaxRight = GetIndex(b+1, right);
  int indexMinRight = GetIndex(b, right);
  double ff = GG[indexMinLeft] * (1.0 - alpha) + GG[indexMinRight] * alpha;
  double gg = GG[indexMaxLeft] * (1.0 - alpha) + GG[indexMaxRight] * alpha;
  
  // handling numerical inaccuracy
  const float epsilon = 1e-15;
  if (radius < sinMin - epsilon) {
    radius = sinMin;
    cerr << "AAAA" << endl;
    return false;
  }
  if (radius > sinMax + epsilon) {
    radius = sinMax;
    cerr << "BBBB" << endl;
    return false;
  }
  
  // cout << "rr - ff = " << rr - ff << " gg - rr = " << gg - rr << endl;
  
  // now the right patch is found, perform the stretching itself
  
  // now interpolate, since it is always possible
  // d and e is required density/intensity
  double d = INTENSITY[GetIndex(b, left)] * (1.0 - alpha) +
             INTENSITY[GetIndex(b, right)] * alpha;
  double e = INTENSITY[GetIndex(b+1, left)] * (1.0 - alpha) +
             INTENSITY[GetIndex(b+1, right)] * alpha;
  double rr = 0.0;
  
  // constant density/intensity in the patch for given phi ?
  if (fabs(d-e) < 1e-12) {
    // constant remapping ... easy case
    if (fabs(d) > 1e-12)
      rr = ((radius - SIN_VERTICAL_ANGLES[b]) * d + ff);
    else
      rr = SIN_VERTICAL_ANGLES_SQR[b]; // quite unlikely

#ifdef _DEBUG
    double a = SIN_VERTICAL_ANGLES[b];
    double c = SIN_VERTICAL_ANGLES[b+1];
    if (rr < a) {
      cout << "error A radius = " << radius << " a = " << a << " c = "
	   << c << " band = " << b
	   << " lmin = " << INTENSITY[GetIndex(b, left)]
	   << " rmin = " << INTENSITY[GetIndex(b, right)]
	   << " lmax = " << INTENSITY[GetIndex(b+1, left)]
	   << " rmax = " << INTENSITY[GetIndex(b+1, right)]
	   << " d = " << d << " e = " << e << " alpha = " << alpha << endl;
      return true; // error
    }

    if (rr > c) {
      cout << "error B radius = " << radius << " a = " << a << " c = "
	   << c << " band = " << b
	   << " lmin = " << INTENSITY[GetIndex(b, left)]
	   << " rmin = " << INTENSITY[GetIndex(b, right)]
	   << " lmax = " << INTENSITY[GetIndex(b+1, left)]
	   << " rmax = " << INTENSITY[GetIndex(b+1, right)]
	   << " d = " << d << " e = " << e << " alpha = " << alpha << endl;
      return true; // error ... sample is not usable
    }
#endif // _DEBUG

    radius = rr;
  }
  else { // different densities/intensities
    // linear remapping, it corresponds to solving a cubic equation
    double a = SIN_VERTICAL_ANGLES[b];
    double c = SIN_VERTICAL_ANGLES[b+1];
    double aa = SIN_VERTICAL_ANGLES_SQR[b]; // or sqr(a);
    double tt = (c - a)/(d - e);
    // Prepare input parameters to compute the equation
    double B = (-3.0 * (e * tt + c))/2.0;

    rr = ff - ( (radius*radius*radius + B * radius * radius) * 2.0/3.0 +
		aa * (e * tt + c) - aa*a * 2.0/3.0) / tt;
    assert(rr >= 0.0);
    // compute the new value of radius
    radius = sqrt(rr);

#ifdef _DEBUG
    if (rr < ff) {
      cout << "error C radius = " << radius
	   << " phi = " << phi
	   << " a = " << a << " c = " << c
	   << " band = " << b
	   << " lmin = " << INTENSITY[GetIndex(b, left)]
	   << " rmin = " << INTENSITY[GetIndex(b, right)]
	   << " lmax = " << INTENSITY[GetIndex(b+1, left)]
	   << " rmax = " << INTENSITY[GetIndex(b+1, right)]
	   << " d = " << d << " e = " << e << " alpha = " << alpha << endl;
      return true; // error
    }

    // handling numerical inaccuracy
    const float epsilon = 1e-15;
    if (rr > gg + epsilon) {
      // this can rarely happen (inaccuracy)
      cout << " index = " << index << endl;
      cout << "error D radius = " << radius
	   << " phi = " << phi
	   << " a = " << a << " c = " << c
	   << " band = " << b
	   << " lmin = " << INTENSITY[GetIndex(b, left)]
	   << " rmin = " << INTENSITY[GetIndex(b, right)]
	   << " lmax = " << INTENSITY[GetIndex(b+1, left)]
	   << " rmax = " << INTENSITY[GetIndex(b+1, right)]
	   << " d = " << d << " e = " << e << " alpha = " << alpha << endl;
      return true; // error ... sample is not usable
    }
#endif // _DEBUG
  } // different densities/intensities

  // modify also radius by effective density for given 'A' and thus 'alpha'
  radius *= sqrt(effectiveDensity);
  
  // we now have radius in the range <0.0-MULT>
  // and phi in the range <0.0-2*PI>
  radius /= MULT;

#ifdef _DEBUG
  if (radius > MULT) {
    cout << "error E" << " radius = " << radius << " phi = " << phi
	 << " effDensity = " << effectiveDensity << endl;
    // this sample cannot be used, some bug or numerical inaccuracy
    // It should never happen!
    return true; // error
  }
#endif  
  
  return false; // valid sample
}

// Given the vector (vecZ >= 0.0), computes the interpolated intensity with
// respect to the definition of goniometric diagram as the linear interpolation
// between the four corners of a patch on the goniometric web.
float
CHemisphereGonio::GetIntensity(double vecX, double vecY, double vecZ)
{
  double magnitude = sqrt(vecX * vecX + vecY * vecY + vecZ * vecZ);
  assert(magnitude > 0.0);
  magnitude = 1.0 / magnitude;
  vecX *= magnitude;
  vecY *= magnitude;
  vecZ *= magnitude;
  assert(vecZ<=1.0);
  assert(vecZ>=0.0);

  // the value in the horizontal direction
  double phi = 0.0;
  
  // now rad is in the range <0.0-1.0>
  if (vecZ > 1e-20) {
    double cosPHI = vecX / vecZ;
    double sinPHI = vecY / vecZ;
    // compute phi by arcusstangent over the two arguments
    phi = atan2(sinPHI, cosPHI);
    if (phi < 0.0)
      phi = 2.0 * M_PI + phi;
  }

  // so 'phi' is computed, range is <0.0-2*PI>

  // now compute 'theta' angle
  double theta = acos(vecZ);
  // now theta is in range <0.0-PI/2>
  
  // compute now the indeces - the patch position where the
  // direction belongs to.
  int left, band;

  // For high number of 'phi' angles in subdivision, binary search would
  // be much faster!
  for (left = 0; left < HORIZONTAL_COUNT; left++) {
    if (HORIZONTAL_ANGLES[left+1] > phi)
      break;
  }
  assert(left < HORIZONTAL_COUNT);
  // the interpolation coeficient in horizontal direction (phi)
  float alpha = (HORIZONTAL_ANGLES[left+1] - phi) /
    (HORIZONTAL_ANGLES[left+1] - HORIZONTAL_ANGLES[left]);
  int right = (left+1) % HORIZONTAL_COUNT;

  // For high number of 'theta' angles in subdivision, binary search would
  // be much faster!
  for (band = 0; band < VERTICAL_COUNT-1; band++) {
    if (VERTICAL_ANGLES[band+1] > theta)
      break;
  }
  // check if for the hemisphere not specified till the equator - if we
  // are not out of range
  if (band >= VERTICAL_COUNT-1)
    return 0.0; // yes, we are, zero intensity
  
  // the interpolation coeficient in vertical direction (theta)
  float beta = (VERTICAL_ANGLES[band+1] - theta) /
    (VERTICAL_ANGLES[band+1] - VERTICAL_ANGLES[band]);

  // now compute the intensity
  double intensity =
    alpha * (beta * INTENSITY[GetIndex(band, left)] +
	     (1.0 - beta) * INTENSITY[GetIndex(band+1, left)])
    + (1.0 - alpha) * (beta * INTENSITY[GetIndex(band, right)] +
		       (1.0 - beta) * INTENSITY[GetIndex(band+1, right)]);

  return intensity; // return the interpolated value
}


CHemisphereGonio *CHemisphereGonio::Dup() const
{
  CHemisphereGonio *ret = new CHemisphereGonio;

  if(VERTICAL_COUNT)
    ret->VERTICAL_COUNT = VERTICAL_COUNT;
  if(HORIZONTAL_COUNT)
    ret->HORIZONTAL_COUNT = HORIZONTAL_COUNT;

  if(VERTICAL_ANGLES)
  {
    ret->VERTICAL_ANGLES = new double[VERTICAL_COUNT];
    ret->VERTICAL_ANGLES = VERTICAL_ANGLES;
  }

  if(HORIZONTAL_ANGLES)
  {
    ret->HORIZONTAL_ANGLES = new double[HORIZONTAL_COUNT+1];
    ret->HORIZONTAL_ANGLES = HORIZONTAL_ANGLES;
  }

  if(INTENSITY)
  {
    ret->INTENSITY = new double[VERTICAL_COUNT*HORIZONTAL_COUNT];
    ret->INTENSITY = INTENSITY;
  }

  if(SIN_VERTICAL_ANGLES)
  {
    ret->SIN_VERTICAL_ANGLES = new double[VERTICAL_COUNT];
    ret->SIN_VERTICAL_ANGLES = SIN_VERTICAL_ANGLES;
  }

  if(SIN_VERTICAL_ANGLES_SQR)
  {
    ret->SIN_VERTICAL_ANGLES_SQR = new double[VERTICAL_COUNT];
    ret->SIN_VERTICAL_ANGLES_SQR = SIN_VERTICAL_ANGLES_SQR;
  }

  if(MULT) ret->MULT = MULT;

  if(MULT_index) ret->MULT_index = MULT_index;

  if(GG)
  {
    ret->GG = new double[VERTICAL_COUNT*HORIZONTAL_COUNT];
    ret->GG = GG;
  }

  if(ALPHA)
  {
    ret->ALPHA = new double[HORIZONTAL_COUNT+1];
    ret->ALPHA = ALPHA;
  }

  if(_maxPhiRescale) ret->_maxPhiRescale = _maxPhiRescale;

  if(last_band) ret->last_band = last_band;
  if(last_sector) ret->last_sector = last_sector;

  if(maximumIntensity) ret->maximumIntensity = maximumIntensity;

  if(totalPower) ret->totalPower = totalPower;

  return ret;
}

__END_BRDFLIB_SOURCE
