// -------------------------------------------------
// FILE NAME   bigdlight.cpp
// PURPOSE     Interface with goniodiagram point light
// CREATED     '00.07.01 Rkd (Kirill Dmitriev)
//
// This file is a property of Rkd (Kirill Dmitriev),
// do not change or distribute it without permission
// of Rkd.
// -------------------------------------------------

#ifdef WIN32
#include <windows.h>
#endif
#include <GL/gl.h>

//#include <biglobal/biglobal.h>
#include <base/time.hpp>
#include <math/sort.hpp>
#include <math/solve.hpp>

#include "bigdlight.h"
#include <bifixture/bifixture.h>
#include <bifmcrtray/bifmcrtray.h>
#include <bisampler/bisampler.h>
#include <bihit/bihit.h>

#include <bidrafter/bicontext.h>
#include <bigeforce/bigeforce.h>

biPointer _stdcall biGDLight::Init(void)
  { return new biGDLight; }
biType *biGDLight::st_type = NULL;

biGDLight::biGDLight() : biPointLight(biGDLight::getStaticType())
  {
  }
biGDLight::biGDLight(biType *type0) : biPointLight(type0)
  {
  }

biGDLight::~biGDLight()
  {
  }

extern void GenerateSample(double vx, double vy, double *dir);
extern void InitIntensities(int phi_num, double *phi, int theta_num, double *theta, double *diagram);

biType *biGDLight::getStaticType()
  {
  if (st_type == NULL)
    st_type = new biType(biLight::getStaticType(), BI_GD_LIGHT, Init);
  return st_type;
  }

#define IS_INDIRECT

// Switch between Vlastimil's approach (NUMERIC undefined)
// and Kirill's approach (NUMERIC DEFINED)
//#define NUMERIC

void biGDLight::FMCRTShootRay(FixtureRef &fixture, QRNG &rng, biFMCRTRay &ray) const
  {
  double dtheta = rng.QRnd();
  double dphi = rng.QRnd();

#ifdef NUMERIC
  int iphi, itheta, ii;
  //  choose iphi and itheta of the cell
  iphi = TSort<ThetaLine>::BinarySearchD(grid.Data(), phi_num, dphi) - 1;
  dphi = (dphi - grid[iphi].probab) / (grid[iphi + 1].probab - grid[iphi].probab);
  CELL *cells = grid[iphi].cells;
  itheta = TSort<CELL>::BinarySearchD(cells, theta_num, dtheta) - 1;
  dtheta = (dtheta - cells[itheta].probab) / (cells[itheta + 1].probab - cells[itheta].probab);

  //  aquire the intensity values in the corners
  double a1 = diagram[ii = iphi * theta_num + itheta];
  double b1 = diagram[ii + 1];
  double a2 = diagram[ii += theta_num];
  double b2 = diagram[ii + 1];

#ifdef FAST_SHOOT
  double costh;
  a1 = ((b1 * (1.0 - dphi) + b2 * dphi) * dtheta + (a1 * (1.0 - dphi) + a2 * dphi) * (1.0 - dtheta));
  dphi = phi[iphi] + (phi[iphi + 1] - phi[iphi]) * dphi;
  costh = cos_theta[itheta] + (cos_theta[itheta + 1] - cos_theta[itheta]) * dtheta;
  ray.SetWeight(a1 / cells[itheta].area_func);
#else
  double data[9];

  //  numerically solve non-linear equation to find theta, equation is:
  //  ((a1 + a2) * (-Sin(DT1) + Sin(dt0) + DT * Cos(dt0)) + (B1 + B2) * ( Sin(DT1) - Sin(dt0) - DT * Cos(DT1))) / DT - 2 * cell_area * dtheta / dp = 0
  //  (a1 + a2) * Cos(dt0) - (B1 + B2) * Cos(DT1) - (Sin(DT1) - Sin(dt0)) * (a1 + a2 - (B1 + B2)) / DT - 2 * cell_area * dtheta / dp = 0
  if (dtheta < MathF::TOLERANCE)
    {
    dtheta = theta[itheta] + dtheta * (theta[itheta + 1] - theta[itheta]);
    }
  else
    {
    data[0] = theta[itheta];
    data[1] = theta[itheta + 1];
    data[2] = data[1] - data[0];
    data[3] = a1 + a2;
    data[4] = b1 + b2;
    data[5] = Sin(data[0]);
    data[6] = Cos(data[0]);
    data[7] = data[3] * data[6] - cells[itheta].area_func * dtheta;
    dtheta = Solve::FindRoot((double (*)(void *actor, double x))area_func, data, data[0], -cells[itheta].area_func * dtheta, data[1], cells[itheta].area_func * (1 - dtheta));
    }
  dphi = phi[iphi] + dphi * (phi[iphi + 1] - phi[iphi]);
  ray.SetWeight(1);
#endif
  //  setup ray parameters
  a1 = Sqrt(1 - Sqr(costh));
  Vect3d dir(a1 * Cos(dphi), a1 * Sin(dphi), -costh);
#else
  Vect3d dir;
  //  initialize vlastimil's gd sampler
  GenerateSample(dphi, dtheta, dir.x);
  ray.SetWeight(1);
#endif
  ray.SetOrig(fixture->ToWorldGlobal()[3].Point());
  ray.SetDir(dir * fixture->ToWorldGlobalWoScale());
  ray.SetColor(color);
  ray.eta = 1;

#ifdef IS_INDIRECT
  ray.flags = (UINT)biBRTRay::FL_INDIRECT;
#else
  ray.flags = 0;
#endif
  }

void biGDLight::ShadeSampling(biSampler &sampler, biFixture *fixture, biRay &ray, biHit &hit, RGBColor &lum) const
  {
#ifndef IS_INDIRECT
  Vect3d dir = fixture->ToWorldGlobal()[3] - hit.Point().Vector();
  double phi0, theta0, xx, yy, a1, a2, b1, b2, len = DotProd(dir, dir);
  int itheta, iphi, icell;

  theta0 = DotProd(dir, fixture->ToWorldGlobalWoScale()[2]);
  if (MathF::AboutEqual(theta0 * theta0, len))
    {
    if (theta0 > 0)
      {
      if (theta[0] >= 0 || (xx = diagram[0]) == 0) return;
      }
    else if (theta[theta_num - 1] < PI || (xx = diagram[theta_num - 1]) == 0)
      {
      return;
      }
    len = Sqrt(len);
    sampler.ShadeDirectional(ray, hit, dir /= len, len, ccolor * (xx / (len * len)), lum);
    return;
    }
  if ((theta0 = ACos(theta0 / (len = Sqrt(len)))) <= theta[0] || theta0 > theta[theta_num - 1])
    return;

  xx = DotProd(dir, fixture->ToWorldGlobalWoScale()[0]);
  yy = DotProd(dir, fixture->ToWorldGlobalWoScale()[1]);
  phi0 = Sqrt(xx * xx + yy * yy);
  if (yy < 0)
    phi0 = PI - ACos(xx / phi0);
  else
    phi0 = PI + ACos(xx / phi0);
  if (phi0 <= phi[0] || phi0 > phi[phi_num - 1])
    return;

  itheta = TSort<double>::BinarySearch(theta.Data(), theta_num, theta0) - 1;
  iphi = TSort<double>::BinarySearch(phi.Data(), phi_num, phi0) - 1;
  dir /= len;

  //  cell intensity is constant
  if (grid[iphi].cells[itheta].max_intensity == 0)
    return;
  if (grid[iphi].cells[itheta].max_intensity  < 0)
    {
    sampler.ShadeDirectional(ray, hit, dir, len, ccolor * (-grid[iphi].cells[itheta].max_intensity / (len * len)), lum);
    }
  //  cell intensity is variable
  else
    {
    xx = (phi0 - phi[iphi]) / (phi[iphi + 1] - phi[iphi]);
    yy = (theta0 - theta[itheta]) / (theta[itheta + 1] - theta[itheta]);
    a1 = diagram[icell = iphi * theta_num + itheta];
    b1 = diagram[icell + 1];
    a2 = diagram[icell += theta_num];
    b2 = diagram[icell + 1];
    xx = ((b1 * (1.0 - xx) + b2 * xx) * yy + (a1 * (1.0 - xx) + a2 * xx) * (1.0 - yy));
    sampler.ShadeDirectional(ray, hit, dir, len, ccolor * (xx / (len * len)), lum);
    }
#endif
  }

void biGDLight::GLApply(biFixture *fixture, int index) const
  {
  float data[4] = {0, 0, 0, 1};

  if(biContext::crrnt_cxt->ge_force3)
    {
    RGBColor specular(ccolor * max_intensity * biContext::crrnt_cxt->adp_mult);
#ifdef IS_INDIRECT
    RGBColor diffuse(0);
#else
    RGBColor diffuse = specular / PI;
#endif
    Vect3f::Cast(data) = ConvF(fixture->ToWorldGlobal()[3]);

    //  FIXME: hack ...
    glLightfv(GL_LIGHT0 + index, GL_POSITION, data);
    glGetLightfv(GL_LIGHT0 + index, GL_POSITION, data);
    Vect3f position = Vect3f::Cast(data);

    biContext::crrnt_cxt->ge_force3->setLightParameter(index, position, diffuse, specular, false);
    }
  else
    {
    Vect3f::Cast(data) = ccolor * (max_intensity * biContext::crrnt_cxt->adp_mult);
    glLightfv(GL_LIGHT0 + index, GL_SPECULAR, data);
#ifdef IS_INDIRECT
    Vect3f::Cast(data).Set(0);
#else
    Vect3f::Cast(data) /= PI; //  for  diffuse component we have to devide luminance by PI
#endif
    glLightfv(GL_LIGHT0 + index, GL_DIFFUSE, data);

    Vect3f::Cast(data) = ConvF(fixture->ToWorldGlobal()[3]);
    glLightfv(GL_LIGHT0 + index, GL_POSITION, data);

    //  define goniodiagram as spot light
    Vect3f::Cast(data) = ConvF(-fixture->ToWorldGlobal()[2]);
    glLightfv(GL_LIGHT0 + index, GL_SPOT_DIRECTION, data);
    glLightf(GL_LIGHT0 + index, GL_SPOT_CUTOFF, 90);
    glLightf(GL_LIGHT0 + index, GL_SPOT_EXPONENT, Clipped((PI2 * max_intensity) / flux - 1, 0, 128));

    glLightf(GL_LIGHT0 + index, GL_CONSTANT_ATTENUATION, 0);
    glLightf(GL_LIGHT0 + index, GL_LINEAR_ATTENUATION, 0);
    glLightf(GL_LIGHT0 + index, GL_QUADRATIC_ATTENUATION, 1);
    }
  }

void biGDLight::SetFlux(const double flux0)
  {
  int iphi, itheta;

  //  goniodiagram must be rescaled
  for (iphi = phi_num * theta_num - 1; iphi >= 0; --iphi)
    {
    diagram[iphi] *= flux0 / flux;
    }
  //  cells area_func and max_intensity values must be rescaled
  for (itheta = theta_num - 2; itheta >= 0; --itheta)
    {
    for (iphi = phi_num - 2; iphi >= 0; --iphi)
      {
      grid[iphi].cells[itheta].area_func *= flux0 / flux;
      grid[iphi].cells[itheta].max_intensity *= flux0 / flux;
      }
    }
  //  global max_intensity value must be rescaled
  max_intensity *= flux0 / flux;

  biPointLight::SetFlux(flux0);
  }

const Str IESNA_91_STR("IESNA91");
const Str IESNA_95_STR("IESNA:LM-63-1995");

#define TILT          "TILT="
#define TILT_NONE     "TILT=NONE"
#define TILT_INCLUDE  "TILT=INCLUDE"

#define IES_UNDEFINED             0
#define IES_AXISYMMETRICAL        1   /* 0-0    */
#define IES_QUADRANT_SYMMETRICAL  2   /* 0-90   */
#define IES_PLANESYMMETRICAL_1    3   /* 0-180  */
#define IES_PLANESYMMETRICAL_2    4   /* 90-270 */
#define IES_ASYMMETRICAL          5   /* 0-360  */
#define IES_NONSTANDARD           6   /* 0-x, 180< x < 360 */

BOOL biGDLight::ReadGD(const PathStr &file_name)
  {
  Timer debug_timer;
  biGlobal::EnterFunction(NT_DEBUG, "biGDLight::ReadGD()");

  Str line;
  InFile file(file_name);
  int lamps_num, unit;
  double lumens, width, length, height, ball_factor, future_used, watts;

  //  open the IESNA data file
  if (file.Open() != SUCCESS)
    {
    biGlobal::InsideFunction(NT_FATAL, "error reading %s", file_name);
    return FALSE;
    }

  //  read IESNA standard identifier line
  if (file.SkipText(IESNA_91_STR) == SUCCESS)
    file.ReadLine(line);  //  read the rest of line
  if (file.SkipText(IESNA_95_STR) == SUCCESS)
    file.ReadLine(line);  //  read the rest of line

  //  search for TILT= line
  for ( ; ; )
    {
    if (file.ReadLine(line) != SUCCESS)
      {
      biGlobal::InsideFunction(NT_FATAL, "error reading %s", file_name);
      return FALSE;
      }
    if (line.SubStrCmp(TILT) == 0)
      break;
    }

  //  TILT= line was read
  //  skip info specific for TILT=INCLUDE
  if (line.SubStrCmp(TILT_INCLUDE) == 0)
    {
    file.Read(unit);    //  lamp-to-luminaire Geometry
    for (unit = 0, file.Read(unit); unit > 0; --unit)
      file.Read(watts); //  TILT angles
    }

  //  dimension data reading
  if ((file.Read(lamps_num)     != SUCCESS) ||
      (file.Read(lumens)        != SUCCESS) ||
      (file.Read(multiplier)    != SUCCESS) ||
      (file.Read(ies_theta_num) != SUCCESS) ||
      (file.Read(ies_phi_num)   != SUCCESS) ||
      (file.Read((int &)gonio_type) != SUCCESS) ||
      (file.Read(unit)          != SUCCESS) ||
      (file.Read(width)         != SUCCESS) ||
      (file.Read(length)        != SUCCESS) ||
      (file.Read(height)        != SUCCESS) ||
      (file.Read(ball_factor)   != SUCCESS) ||
      (file.Read(future_used)   != SUCCESS) ||
      (file.Read(watts)         != SUCCESS))
    {
    biGlobal::InsideFunction(NT_FATAL, "error reading %s", file_name);
    return FALSE;
    }
  multiplier *= ball_factor * future_used;

  if (!ies_theta.Allocate(ies_theta_num) || !ies_phi.Allocate(ies_phi_num) || !ies_diagram.Allocate(ies_theta_num * ies_phi_num))
    {
    biGlobal::InsideFunction(NT_FATAL, "error reading %s", file_name);
    return FALSE;
    }

  if (!LoadArray(file, ies_theta.Data(), ies_theta_num) ||
      !LoadArray(file, ies_phi.Data(), ies_phi_num) ||
      !LoadArray(file, ies_diagram.Data(), ies_theta_num * ies_phi_num))
    {
    biGlobal::InsideFunction(NT_FATAL, "error reading %s", file_name);
    return FALSE;
    }

  if (!IesnaToGonio())
    {
    biGlobal::InsideFunction(NT_FATAL, "error reading %s", file_name);
    return FALSE;
    }

  biGlobal::LeaveFunction(NT_DEBUG, "biGDLight::ReadGD() done in %d ms", debug_timer.Elapsed());
  return TRUE;
  }

void biGDLight::EditDone()
  {
  phi_num = phi.Length();
  theta_num = theta.Length();

  int itheta, iphi, ii;
  double fluxp, fluxt, a1, a2, b1, b2, dp, dt, dt0, dt1, aa, bb;

  grid.Allocate(phi_num);
  for (iphi = phi_num - 1; iphi >= 0; --iphi)
    {
    delete []grid[iphi].cells;
    grid[iphi].cells = new CELL[theta_num];
    }

  //             ^ theta
  //             |
  //          b1 ---------------- b2
  //             |               |
  //             |               |
  //             |               |
  //             ----------------    ->phi
  //          a1                  a2

  for (fluxp = 0, max_intensity = 0, iphi = 0; iphi < phi_num - 1; ++iphi)
    {
    for (fluxt = 0, itheta = 0; itheta < theta_num - 1; ++itheta)
      {
      a1 = diagram[ii = iphi * theta_num + itheta];
      b1 = diagram[ii + 1];
      a2 = diagram[ii += theta_num];
      b2 = diagram[ii + 1];

      if ((dt = Max(Max(a1, a2), Max(b1, b2))) > max_intensity)
        max_intensity = dt;
      if ((dt == a1) && (dt == a2) && (dt == b1) && (dt == b2))
        grid[iphi].cells[itheta].max_intensity = -dt;
      else
        grid[iphi].cells[itheta].max_intensity =  dt;

      dt = (dt1 = theta[itheta + 1]) - (dt0 = theta[itheta]);
      dp = phi[iphi + 1] - phi[iphi];

      aa = 0.5 * dp / dt * (-Sin(dt1) + Sin(dt0) + dt * Cos(dt0));
      bb = 0.5 * dp / dt * ( Sin(dt1) - Sin(dt0) - dt * Cos(dt1));

#ifdef FAST_SHOOT
      grid[iphi].cells[itheta].area_func = (dt = (a1 + a2) * aa + (b1 + b2) * bb) / (2 * aa + 2 * bb);
#else
      grid[iphi].cells[itheta].area_func = (dt = (a1 + a2) * aa + (b1 + b2) * bb) * 2 / dp;
#endif
      grid[iphi].cells[itheta].probab = fluxt;
      fluxt += dt;
      }
    //  normalize probabilities array
    grid[iphi].cells[0].probab = -MathF::TOLERANCE;
    for (itheta = 1; itheta < theta_num - 1; ++itheta)
      {
      grid[iphi].cells[itheta].probab /= fluxt;
      }
    grid[iphi].cells[itheta].probab = 1;

    grid[iphi].probab = fluxp;
    fluxp += fluxt;
    }
  //  normalize probabilities array
  grid[0].probab = -MathF::TOLERANCE;
  for (iphi = 1; iphi < phi_num - 1; ++iphi)
    {
    grid[iphi].probab /= fluxp;
    }
  grid[iphi].probab = 1;

  //  fill sinth array
  cos_theta.Allocate(theta_num);
  for (itheta = theta_num - 1; itheta >= 0; --itheta)
    {
    cos_theta[itheta] = theta[itheta] > 0 ? Cos(theta[itheta]) : 1;
    }

  //  now we will not need any checks after BinarySearch
  phi[0] -= MathF::TOLERANCE;
  theta[0] -= MathF::TOLERANCE;

  //  update light flux according to goniodiagram
  biPointLight::SetFlux(fluxp);

#ifndef NUMERIC
  //  initialize vlasimils gd sampler
  InitIntensities(phi_num, phi.Data(), theta_num, theta.Data(), diagram.Data());
#endif
  }

INLINE int biGDLight::IES_IND(int itheta, int iphi)
  { return iphi * ies_theta_num + itheta; }
INLINE int biGDLight::IND(int itheta, int iphi)
  { return iphi * theta_num + itheta; }

BOOL biGDLight::LoadArray(InFile &file, double *array, int num)
  {
  for (int i = 0; i < num; i++)
    {
    //  read the next float number
    if (file.Read(array[i]) != SUCCESS)
      {
      return FALSE;
      }
    //  skip commas
    file.SkipText(",");
    }
  return TRUE;
  }

BOOL biGDLight::IesnaToGonio()
  {
  if (ies_theta_num == 0 || ies_phi_num == 0)
    return FALSE;

  double dt;
  int iphi, itheta, ii, hor_symmetry = IES_UNDEFINED;

  //  check that theta angles are ok
  if (ies_theta[0] == 0.0)
    {
    if ((ies_theta[ies_theta_num - 1] != 90.0) && (ies_theta[ies_theta_num - 1] != 180.0))
      return FALSE;
    }
  else if (ies_theta[0] == 90.0)
    {
    if (ies_theta[ies_theta_num - 1] != 180.0)
      return FALSE;
    }
  else
    {
    return FALSE;
    }
  theta_num = ies_theta_num;

  if (ies_phi[0] == 0.0)
    {
    if (ies_phi[ies_phi_num - 1] == 0.0)
      {
      hor_symmetry = IES_AXISYMMETRICAL;
      phi_num = 2;
      }
    else if (ies_phi[ies_phi_num - 1] == 90.0)
      {
      hor_symmetry = IES_QUADRANT_SYMMETRICAL;
      phi_num = 4 * ies_phi_num - 3;
      }
    else if (ies_phi[ies_phi_num - 1] == 180.0)
      {
      hor_symmetry = IES_PLANESYMMETRICAL_1;
      phi_num = 2 * ies_phi_num - 1;
      }
    else if (ies_phi[ies_phi_num - 1] == 360.0)
      {
      hor_symmetry = IES_ASYMMETRICAL;
      phi_num = ies_phi_num;
      }
    else if ((ies_phi[ies_phi_num - 1] > 180.0) && (ies_phi[ies_phi_num - 1] < 360.0))
      {
      hor_symmetry = IES_NONSTANDARD;
      phi_num = ies_phi_num + 1;
      }
    }
  else if (ies_phi[0] == 90.0 && ies_phi[ies_phi_num - 1] == 270.0)
    {
    hor_symmetry = IES_PLANESYMMETRICAL_2;
    phi_num = 2 * ies_phi_num - 1;
    }
    
  if (hor_symmetry == IES_UNDEFINED)
    return FALSE;

  if (!theta.Allocate(theta_num) || !phi.Allocate(phi_num) || !diagram.Allocate(theta_num * phi_num))
    return FALSE;

  //  theta values setting
  for (itheta = 0; itheta < theta_num; itheta++)
    theta[itheta] = Rad(ies_theta[itheta]);

  //  phi values setting
  switch (hor_symmetry)
    {
    case IES_AXISYMMETRICAL:
      phi[0] = 0;
      phi[1] = 2.0 * PI;
      break;
    case IES_QUADRANT_SYMMETRICAL:     //  0.0 >= ies_phi <= 90.0
      for (iphi = 0; iphi < ies_phi_num; ++iphi)
        {
        Assert(2 * ies_phi_num - 2 - iphi < phi_num);
        Assert(2 * ies_phi_num - 2 + iphi < phi_num);
        Assert(phi_num - iphi - 1 < phi_num);

        phi[iphi] = ies_phi[iphi];
        phi[2 * ies_phi_num - 2 - iphi] = 180.0 - ies_phi[iphi];
        phi[2 * ies_phi_num - 2 + iphi] = 180.0 + ies_phi[iphi];
        phi[phi_num - 1 - iphi] = 360.0 - ies_phi[iphi];
        }
      break;
    case IES_PLANESYMMETRICAL_1:       //  0.0 >= ies_phi <= 180.0
      for (iphi = 0; iphi < ies_phi_num; iphi++)
        {
        phi[iphi] = ies_phi[iphi];
        phi[phi_num - 1 - iphi] = 360.0 - ies_phi[iphi];
        }
      break;
    case IES_ASYMMETRICAL:
      memcpy(phi.Data(), ies_phi.Data(), phi_num * sizeof(double));
      break;
    case IES_PLANESYMMETRICAL_2:
      ii = ies_phi_num / 2;
      for (iphi = 0; iphi < ies_phi_num; iphi++)
        {
        phi[ii + iphi] = ies_phi[iphi];
        if (iphi <= ii)
          phi[ii - iphi] = 90.0 - ies_phi[iphi];
        if (iphi >= ii)
          phi[phi_num - (iphi - ii)] = 360.0 - ies_phi[iphi];
        }
      break;
    case IES_NONSTANDARD:
      memcpy(phi.Data(), ies_phi.Data(), ies_phi_num * sizeof(double));
      phi[ies_phi_num] = 360.0;
      break;
    }
  for (iphi = 0; iphi < phi_num; ++iphi)
    phi[iphi] = Rad(phi[iphi]);

  //  diagram values setting
  if (hor_symmetry == IES_ASYMMETRICAL)
    memcpy(diagram.Data(), ies_diagram.Data(), phi_num * theta_num * sizeof(double));
  else if (hor_symmetry == IES_NONSTANDARD)
    {
    memcpy(diagram.Data(), ies_diagram.Data(), ies_phi_num * ies_theta_num * sizeof(double));
    memcpy(&diagram[ies_phi_num * ies_theta_num], ies_diagram.Data(), ies_theta_num * sizeof(double));
    }
  else
    {
    for (itheta = 0; itheta < theta_num; itheta++)
      {
      switch (hor_symmetry)
        {
        case IES_AXISYMMETRICAL:
          diagram[IND(itheta, 0)] = ies_diagram[IES_IND(itheta, 0)];
          diagram[IND(itheta, 1)] = ies_diagram[IES_IND(itheta, 0)];
          phi[1] = 2.0 * PI;
          break;
        case IES_QUADRANT_SYMMETRICAL:       //  0.0 >= ies_phi <= 90.0
          for (iphi = 0; iphi < ies_phi_num; iphi++)
            {
            Assert(IND(itheta, iphi) < phi_num * theta_num);
            Assert(IND(itheta, 2 * ies_phi_num - 2 - iphi) < phi_num * theta_num);
            Assert(IND(itheta, 2 * ies_phi_num - 2 + iphi) < phi_num * theta_num);
            Assert(IND(itheta, phi_num - 1 - iphi) < phi_num * theta_num);

            diagram[IND(itheta, iphi)] = (dt = ies_diagram[IES_IND(itheta, iphi)]);
            diagram[IND(itheta, 2 * ies_phi_num - 2 - iphi)] = dt;
            diagram[IND(itheta, 2 * ies_phi_num - 2 + iphi)] = dt;
            diagram[IND(itheta, phi_num - 1 - iphi)] = dt;
            }
          break;
        case IES_PLANESYMMETRICAL_1:       //  0.0 >= ies_phi <= 180.0
          for (iphi = 0; iphi < ies_phi_num; iphi++)
            {
            diagram[IND(itheta, iphi)] = (dt = ies_diagram[IES_IND(itheta, iphi)]);
            diagram[IND(itheta, phi_num - 1 - iphi)] = dt;
            }
          break;
        case IES_PLANESYMMETRICAL_2:
          ii = ies_phi_num / 2;
          for (iphi = 0; iphi < ies_phi_num; iphi++)
            {
            diagram[IND(itheta, ii + iphi)] = (dt = ies_diagram[IES_IND(itheta, iphi)]);
            if (iphi <= ii)
              diagram[IND(itheta, ii - iphi)] = dt;
            if (iphi >= ii)
              diagram[IND(itheta, phi_num - (iphi - ii))] = dt;
            }
          break;
        default:
          Assert(FALSE);
        }
      }
    }

  //  initialize cells array
  EditDone();

  return TRUE;
  }

double biGDLight::area_func(double *data, double x)
  {
  double b_sum = (data[3] * (data[1] - x) + data[4] * (x - data[0])) / data[2];
  return data[7] - b_sum * Cos(x) - (Sin(x) - data[5]) * (data[3] - b_sum) / (x - data[0]);
  }
