/****************************************************************************
*                   media.c
*
*  This module contains all functions for participating media.
*
*  from Persistence of Vision(tm) Ray Tracer
*  Copyright 1996,1999 Persistence of Vision Team
*---------------------------------------------------------------------------
*  NOTICE: This source code file is provided so that users may experiment
*  with enhancements to POV-Ray and to port the software to platforms other
*  than those supported by the POV-Ray Team.  There are strict rules under
*  which you are permitted to use this file.  The rules are in the file
*  named POVLEGAL.DOC which should be distributed with this file.
*  If POVLEGAL.DOC is not available or for more info please contact the POV-Ray
*  Team Coordinator by email to team-coord@povray.org or visit us on the web at
*  http://www.povray.org. The latest version of POV-Ray may be found at this site.
*
* This program is based on the popular DKB raytracer version 2.12.
* DKBTrace was originally written by David K. Buck.
* DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
*
*****************************************************************************/

#include "frame.h"
#include "vector.h"
#include "povproto.h"
#include "chi2.h"
#include "colour.h"
#include "povray.h"
#include "texture.h"
#include "pigment.h"
#include "objects.h"
#include "lighting.h"
#include "matrices.h"
#include "media.h"
#include "pattern.h"
#include "povray.h"
#include "point.h"
#include "texture.h"
#include "ray.h"



/*****************************************************************************
* Local preprocessor defines
******************************************************************************/

/*****************************************************************************
* Local typedefs
******************************************************************************/

typedef struct Light_List_Struct LIGHT_LIST;
typedef struct Media_Interval_Struct MEDIA_INTERVAL;
typedef struct Lit_Interval_Struct LIT_INTERVAL;

struct Lit_Interval_Struct
{
  int lit;
  DBL s0, s1, ds;
};

struct Media_Interval_Struct
{
  int lit;
  int samples;
  DBL s0, s1, ds;
  COLOUR od;
  COLOUR te;
  COLOUR te2;
};

struct Light_List_Struct
{
  int active;
  DBL s0, s1;
  LIGHT_SOURCE *Light;
};



/*****************************************************************************
* Local variables
******************************************************************************/



/*****************************************************************************
* Static functions
******************************************************************************/

static void sample_media (LIGHT_LIST *, RAY *, IMEDIA **, MEDIA_INTERVAL *, int );
static void get_light_list (LIGHT_LIST *, RAY *, INTERSECTION *);
static void get_lit_interval (int *, LIT_INTERVAL *, int, LIGHT_LIST *, INTERSECTION *);
static int set_up_sampling_intervals (MEDIA_INTERVAL *,
  int, LIT_INTERVAL *, IMEDIA *);

static int intersect_spotlight (RAY *Ray, LIGHT_SOURCE *Light, DBL *d1, DBL *d2);
static int intersect_cylinderlight (RAY *Ray, LIGHT_SOURCE *Light, DBL *d1, DBL *d2);

static int CDECL compdoubles (CONST void *in_a, CONST void *in_b);

static void evaluate_density_pattern (IMEDIA *, VECTOR , COLOUR);


/*****************************************************************************
*
* FUNCTION
*
*   Simulate_Media
*
* INPUT
*
*   Ray       - Current ray, start point P0
*   Inter     - Current intersection, end point P1
*   Colour    - Color emitted at P1 towards P0
*   light_ray - TRUE if we are looking at a light source ray
*
* OUTPUT
*
*   Colour    - Color arriving at the end point
*
* RETURNS
*
* AUTHOR
*
*   Dieter Bayer
*
* DESCRIPTION
*
*   Simulate participating media using volume sampling.
*
*   The effects of participating media on the light emitted at P1
*   towards P0 are determined using Monte Carlo integration.
*
*   The effects include: emission, absoprtion and scattering.
*
*   Currently one global medium with constant coefficients is implemented.
*
*   Ideas for the atmospheric scattering were taken from:
*
*     - M. Inakage, "An Illumination Model for Atmospheric Environments", ..
*
*     - Nishita, T., Miyawaki, Y. and Nakamae, E., "A Shading Model for
*       Atmospheric Scattering Considering Luminous Intensity Distribution
*       of Light Sources", Computer Graphics, 21, 4, 1987, 303-310
*
* CHANGES
*
*   Nov 1994 : Creation.
*
*   Jan 1995 : Added support of cylindrical light sources. [DB]
*
*   Jun 1995 : Added code for alpha channel support. [DB]
*
******************************************************************************/

void Simulate_Media(IMEDIA **Media_List, RAY *Ray, INTERSECTION *Inter, COLOUR Colour, int light_ray)
{
  int i, j, intervals, use_extinction;
  int lit_interval_entries;
  DBL n;
  COLOUR Od, Te, Va;
  LIGHT_LIST *Light_List = NULL;
  LIT_INTERVAL *Lit_Interval;
  IMEDIA *IMedia, **Tmp, *Local;
  MEDIA_INTERVAL *Media_Interval, *curr;

  /* Why are we here? */

  if ((Media_List == NULL) || (Media_List[0] == NULL))
  {
    return;
  }

  /* Find media with the largest number of intervals. */

  intervals = 0;

  use_extinction = FALSE;

  IMedia = Media_List[0];

  for (Tmp = Media_List; (*Tmp) != NULL; Tmp++)
  {
    for (Local = *Tmp; Local != NULL; Local = Local->Next_Media)
    {
      if (Local->Intervals > intervals)
      {
        intervals = Local->Intervals;

        IMedia = Local;
      }

      use_extinction |= Local->use_extinction;
    }
  }

  /* If this is a light ray and no extinction is used we can return. */

  if ((light_ray) && (!use_extinction))
  {
    return;
  }

  /*
   * Prepare the Monte Carlo integration along the ray from P0 to P1.
   */

  if (light_ray || (Frame.Number_Of_Light_Sources==0))
  {
    Lit_Interval = (LIT_INTERVAL *)POV_MALLOC(sizeof(LIT_INTERVAL), "lit interval");

    lit_interval_entries = 1;

    Lit_Interval[0].lit = FALSE;

    Lit_Interval[0].s0 = 0.0;
    Lit_Interval[0].s1 =
    Lit_Interval[0].ds = Inter->Depth;
  }
  else
  {
    /* Get light list. */

    Light_List = (LIGHT_LIST *)POV_MALLOC(Frame.Number_Of_Light_Sources*sizeof(LIGHT_LIST), "light list");

    get_light_list(Light_List, Ray, Inter);

    /* Get lit intervals. */

    Lit_Interval = (LIT_INTERVAL *)POV_MALLOC((2*Frame.Number_Of_Light_Sources+1)*sizeof(LIT_INTERVAL), "lit interval");

    get_lit_interval(&lit_interval_entries, Lit_Interval, Frame.Number_Of_Light_Sources, Light_List, Inter);
  }

  /* Set up sampling intervals. */

  Media_Interval = (MEDIA_INTERVAL *)POV_MALLOC(IMedia->Intervals*sizeof(MEDIA_INTERVAL), "media intervals");

  intervals = set_up_sampling_intervals(Media_Interval, lit_interval_entries, Lit_Interval, IMedia);

  /* Sample all intervals. */

  for (i = 0; i < intervals; i++)
  {
    /* Sample current interval. */

    Increase_Counter(stats[Media_Intervals]);

    for (j = 0; j < IMedia->Min_Samples; j++)
    {
      sample_media(Light_List, Ray, Media_List, &Media_Interval[i], light_ray);
    }
  }

  /* Cast additional samples if necessary. */

  if ((!light_ray) && (IMedia->Max_Samples > IMedia->Min_Samples))
  {
    curr = &Media_Interval[0];

    for (i = 0; i < intervals; i++)
    {
      if (curr->samples < IMedia->Max_Samples)
      {
        /* Get variance of samples. */

        n = (DBL)curr->samples;

        Va[0] = (curr->te2[0] / n - Sqr(curr->te[0] / n)) / n;
        Va[1] = (curr->te2[1] / n - Sqr(curr->te[1] / n)) / n;
        Va[2] = (curr->te2[2] / n - Sqr(curr->te[2] / n)) / n;

        /* Take additional samples until variance is small enough. */

        while ((Va[0] >= IMedia->Sample_Threshold[curr->samples-1]) ||
               (Va[1] >= IMedia->Sample_Threshold[curr->samples-1]) ||
               (Va[2] >= IMedia->Sample_Threshold[curr->samples-1]))
        {
          /* Sample current interval again. */

          sample_media(Light_List, Ray, Media_List, curr, light_ray);

          /* Have we reached maximum number of samples. */

          if (curr->samples > IMedia->Max_Samples)
          {
            break;
          }

          /* Get variance of samples. */

          n = (DBL)curr->samples;

          Va[0] = (curr->te2[0] / n - Sqr(curr->te[0] / n)) / n;
          Va[1] = (curr->te2[1] / n - Sqr(curr->te[1] / n)) / n;
          Va[2] = (curr->te2[2] / n - Sqr(curr->te[2] / n)) / n;
        }
      }

      curr++;
    }
  }

  /* Sum the influences of all intervals. */

  Make_Colour(Od, 0.0, 0.0, 0.0);
  Make_Colour(Te, 0.0, 0.0, 0.0);

  curr = &Media_Interval[0];

  for (i = 0; i < intervals; i++)
  {
    /* Add total emission. */

    Te[0] += curr->te[0] / (DBL)curr->samples * exp(-Od[0]);
    Te[1] += curr->te[1] / (DBL)curr->samples * exp(-Od[1]);
    Te[2] += curr->te[2] / (DBL)curr->samples * exp(-Od[2]);

    /* Add optical depth of current interval. */

    Od[0] += curr->od[0] / (DBL)curr->samples;
    Od[1] += curr->od[1] / (DBL)curr->samples;
    Od[2] += curr->od[2] / (DBL)curr->samples;

    curr++;
  }

  /* Add contribution estimated for the participating media. */

  Colour[0] = Colour[0] * exp(-Od[0]) + Te[0];
  Colour[1] = Colour[1] * exp(-Od[1]) + Te[1];
  Colour[2] = Colour[2] * exp(-Od[2]) + Te[2];

  if (!(light_ray || (Frame.Number_Of_Light_Sources==0)))
  {
    POV_FREE(Light_List);
  }

  POV_FREE(Lit_Interval);

  POV_FREE(Media_Interval);
}



/*****************************************************************************
*
* FUNCTION
*
*   sample_media
*
* INPUT
*
*   dist  - distance of current sample
*   Ray   - pointer to ray
*   IMedia - pointer to media to use
*
* OUTPUT
*
*   Col          - color of current sample
*
* RETURNS
*
* AUTHOR
*
*   Dieter Bayer
*
* DESCRIPTION
*
*   Calculate the color of the current media sample.
*
* CHANGES
*
*   Nov 1994 : Creation.
*
******************************************************************************/

static void sample_media(LIGHT_LIST *Light_List, RAY *Ray, IMEDIA **Media_List, MEDIA_INTERVAL *Interval, int light_ray)
{
  int i, n, use_scattering;
  DBL alpha, d0, d1, len, k, g, g2;
  VECTOR P, H;
  COLOUR C0, Light_Colour, Te;
  COLOUR Em, Ex, Sc;
  RAY Light_Ray;
  IMEDIA **Tmp, *Local;

  Increase_Counter(stats[Media_Samples]);

  /* Set up sampling location. */

  d0 = Interval->ds * FRAND();

  d1 = Interval->s0 + d0;

  VEvaluateRay(H, Ray->Initial, d1, Ray->Direction);

  /* Get coefficients in current sample location. */

  Make_Colour(Em, 0.0, 0.0, 0.0);
  Make_Colour(Ex, 0.0, 0.0, 0.0);
  Make_Colour(Sc, 0.0, 0.0, 0.0);

  use_scattering = FALSE;

  for (Tmp = Media_List; (*Tmp) != NULL; Tmp++)
  {
    for (Local = *Tmp; Local != NULL; Local = Local->Next_Media)
    {
      Assign_Vector(P, H);

      evaluate_density_pattern(Local, P, C0);

      Ex[0] += C0[0] * Local->Extinction[0];
      Ex[1] += C0[1] * Local->Extinction[1];
      Ex[2] += C0[2] * Local->Extinction[2];

      if (!light_ray)
      {
        Em[0] += C0[0] * Local->Emission[0];
        Em[1] += C0[1] * Local->Emission[1];
        Em[2] += C0[2] * Local->Emission[2];

        Sc[0] += C0[0] * Local->Scattering[0];
        Sc[1] += C0[1] * Local->Scattering[1];
        Sc[2] += C0[2] * Local->Scattering[2];
      }

      use_scattering |= Local->use_scattering;
    }
  }

  /* Get estimate for the total optical depth of the current interval. */

  Interval->od[0] += Ex[0] * Interval->ds;
  Interval->od[1] += Ex[1] * Interval->ds;
  Interval->od[2] += Ex[2] * Interval->ds;

  /* Get estimate for the total emission of the current interval. */

  Te[0] = Em[0];
  Te[1] = Em[1];
  Te[2] = Em[2];

  if ((!light_ray) && (use_scattering) && (Interval->lit))
  {
    /* Process all light sources. */

    for (i = 0; i < Frame.Number_Of_Light_Sources; i++)
    {
      /* Use light only if active and within it's boundaries. */

      if ((Light_List[i].active) && (d1 >= Light_List[i].s0) && (d1 <= Light_List[i].s1))
      {
        if (!(Test_Shadow(Light_List[i].Light, &len, &Light_Ray, Ray, P, Light_Colour)))
        {
          /* Get attenuation due to scattering. */

          k = 0.0;

          for (n = 0, Tmp = Media_List; (*Tmp) != NULL; n++, Tmp++)
          {
            for (Local = *Tmp; Local != NULL; Local = Local->Next_Media)
            {
              switch (Local->Type)
              {
                case RAYLEIGH_SCATTERING:

                  VDot(alpha, Light_Ray.Direction, Ray->Direction);

                  k += 0.799372013 * (1.0 + Sqr(alpha));

                  break;

                case MIE_HAZY_SCATTERING:

                  VDot(alpha, Light_Ray.Direction, Ray->Direction);

                  k += 0.576655375 * (1.0 + 9.0 * pow(0.5 * (1.0 + alpha), 8.0));

                  break;

                case MIE_MURKY_SCATTERING:

                  VDot(alpha, Light_Ray.Direction, Ray->Direction);

                  k += 0.495714547 * (1.0 + 50.0 * pow(0.5 * (1.0 + alpha), 32.0));

                  break;

                case HENYEY_GREENSTEIN_SCATTERING:

                  VDot(alpha, Light_Ray.Direction, Ray->Direction);

                  g = Local->Eccentricity;

                  g2 = Sqr(g);

                  k += (1.0 - g2) / pow(1.0 + g2 - 2.0 * g * alpha, 1.5);

                  break;

                case ISOTROPIC_SCATTERING:
                default:

                  k += 1.0;

                break;
              }
            }
          }

          k /= (DBL)n;

          Te[0] += k * Sc[0] * Light_Colour[0];
          Te[1] += k * Sc[1] * Light_Colour[1];
          Te[2] += k * Sc[2] * Light_Colour[2];
        }
      }
    }
  }

  Te[0] *= Interval->ds * exp(-Ex[0] * d0);
  Te[1] *= Interval->ds * exp(-Ex[1] * d0);
  Te[2] *= Interval->ds * exp(-Ex[2] * d0);

  /* Add emission. */

  Interval->te[0] += Te[0];
  Interval->te[1] += Te[1];
  Interval->te[2] += Te[2];

  Interval->te2[0] += Sqr(Te[0]);
  Interval->te2[1] += Sqr(Te[1]);
  Interval->te2[2] += Sqr(Te[2]);

  Interval->samples++;
}



/*****************************************************************************
*
* FUNCTION
*
*   evaluate_density_pattern
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
*   Dieter Bayer
*
* DESCRIPTION
*
* CHANGES
*
*   Dec 1996 : Creation.
*
******************************************************************************/

static void evaluate_density_pattern(IMEDIA *IMedia, VECTOR P, COLOUR C)
{
  COLOUR Local_Color;
  PIGMENT *Temp=IMedia->Density;
  
  Make_Colour (C, 1.0, 1.0, 1.0);

  while (Temp != NULL)
  {
    Make_Colour (Local_Color, 0.0, 0.0, 0.0);
    
    Compute_Pigment (Local_Color, Temp, P);

    C[0] *= Local_Color[0];
    C[1] *= Local_Color[1];
    C[2] *= Local_Color[2];

    Temp=(PIGMENT *)Temp->Next;
  }
}



/*****************************************************************************
*
* FUNCTION
*
*   get_light_list
*
* INPUT
*
*   Light_List - array containing light source information
*   dist       - distance of current sample
*   Ray        - pointer to ray
*   IMedia      - pointer to media to use
*
* OUTPUT
*
*   Col        - color of current sample
*
* RETURNS
*
* AUTHOR
*
*   Dieter Bayer
*
* DESCRIPTION
*
* CHANGES
*
*   Nov 1994 : Creation.
*
******************************************************************************/

static void get_light_list(LIGHT_LIST *Light_List, RAY *Ray, INTERSECTION *Inter)
{
  int i, insert;
  DBL t1, t2;
  LIGHT_SOURCE *Light;

  /* Get depths for all light sources and disconnected sampling intervals. */

  t1 = t2 = 0.0;

  for (i = 0, Light = Frame.Light_Sources; Light != NULL; Light = Light->Next_Light_Source, i++)
  {
    /* Init interval. */

    Light_List[i].active = FALSE;
    Light_List[i].s0     = 0.0;
    Light_List[i].s1     = Max_Distance;
    Light_List[i].Light  = NULL;

    insert = FALSE;

    Light_List[i].Light = Light;

    if (!Light->Media_Interaction)
    {
      continue;
    }

    switch (Light->Light_Type)
    {
      case CYLINDER_SOURCE:

        if (intersect_cylinderlight(Ray, Light, &t1, &t2))
        {
          if ((t1 < Inter->Depth) && (t2 > Small_Tolerance))
          {
            insert = TRUE;
          }
        }

        break;

      case POINT_SOURCE:

        t1 = 0.0;
        t2 = Inter->Depth;

        insert = TRUE;

        break;

      case SPOT_SOURCE:

        if (intersect_spotlight(Ray, Light, &t1, &t2))
        {
          if ((t1 < Inter->Depth) && (t2 > Small_Tolerance))
          {
            insert = TRUE;
          }
        }

        break;
    }

    /* Insert distances into sampling interval list. */

    if (insert)
    {
      /* Insert light source intersections into light list. */

      t1 = max(t1, 0.0);
      t2 = min(t2, Inter->Depth);

      Light_List[i].active = TRUE;
      Light_List[i].s0 = t1;
      Light_List[i].s1 = t2;
    }
  }
}



/*****************************************************************************
*
* FUNCTION
*
*   get_lit_interval
*
* INPUT
*
*   Light_List - array containing light source information
*   dist       - distance of current sample
*   Ray        - pointer to ray
*   IMedia      - pointer to media to use
*
* OUTPUT
*
*   Col        - color of current sample
*
* RETURNS
*
* AUTHOR
*
*   Dieter Bayer
*
* DESCRIPTION
*
* CHANGES
*
*   Nov 1994 : Creation.
*
******************************************************************************/

static void get_lit_interval(int *number, LIT_INTERVAL *Lit_Interval, int entries, LIGHT_LIST *Light_List, INTERSECTION *Inter)
{
  int a, i, n;
  DBL *s0, *s1;
  LIT_INTERVAL *curr, *prev;

  s0 = (DBL *)POV_MALLOC(Frame.Number_Of_Light_Sources*sizeof(DBL), "temp data");
  s1 = (DBL *)POV_MALLOC(Frame.Number_Of_Light_Sources*sizeof(DBL), "temp data");

  for (i = a = 0; i < entries; i++)
  {
    if (Light_List[i].active)
    {
      s0[a] = Light_List[i].s0;
      s1[a] = Light_List[i].s1;

      a++;
    }
  }

  n = 0;

  curr = &Lit_Interval[0];

  if (a)
  {
    QSORT((void *)(&s0[0]), (unsigned long)a, sizeof(DBL), compdoubles);
    QSORT((void *)(&s1[0]), (unsigned long)a, sizeof(DBL), compdoubles);

    if (s0[0] > 0.0)
    {
      curr->lit = FALSE;

      curr->s0 = 0.0;
      curr->s1 = s0[0];

      curr++;

      n++;
    }

    curr->lit = TRUE;

    curr->s0 = s0[0];
    curr->s1 = s1[0];

    prev = curr;

    curr++;

    n++;

    for (i = 1; i < a; i++)
    {
      if (s0[i] > prev->s1)
      {
        curr->lit = FALSE;

        curr->s0 = prev->s1;
        curr->s1 = s0[i];

        prev++;
        curr++;

        n++;

        curr->lit = TRUE;

        curr->s0 = s0[i];
        curr->s1 = s1[i];

        prev++;
        curr++;

        n++;
      }
      else
      {
        if (s1[i] > prev->s1)
        {
          prev->s1 = s1[i];
        }
      }
    }

    if (prev->s1 < Inter->Depth)
    {
      curr->lit = FALSE;

      curr->s0 = prev->s1;
      curr->s1 = Inter->Depth;

      curr++;

      n++;
    }
  }
  else
  {
    curr->lit = FALSE;

    curr->s0 = 0.0;
    curr->s1 = Inter->Depth;

    curr++;

    n++;
  }

  curr = &Lit_Interval[0];

  for (i = 0; i < n; i++)
  {
    curr->ds = curr->s1 - curr->s0;

    curr++;
  }

  POV_FREE(s0);
  POV_FREE(s1);

  *number = n;
}



/*****************************************************************************
*
* FUNCTION
*
*   set_up_sampling_intervals
*
* INPUT
*
*   interval - array containing media intervals
*   number   - number of lit intervals
*   list     - array of lit intervals
*   media    - media to use
*
* OUTPUT
*
*   interval - array containing media intervals
*
* RETURNS
*
*   int      - number of media intervals created
*
* AUTHOR
*
*   Dieter Bayer
*
* DESCRIPTION
*
*   Distribute samples along an interval according to the maximum
*   number of samples and the ratio of samples in lit and unlit
*   areas as given by the participating media.
*
* CHANGES
*
*   Nov 1994 : Creation.
*
******************************************************************************/

static int set_up_sampling_intervals(MEDIA_INTERVAL *interval, int number, LIT_INTERVAL *list, IMEDIA *media)
{
  int i, j, n, r, remaining, intervals;
  DBL delta, sum, weight;
  MEDIA_INTERVAL *curr;
  LIT_INTERVAL *entry;

  /* Set up sampling intervals. */

  intervals = media->Intervals;

  /* Use one interval if no lit intervals and constant media. */

  if ((number == 0) && (media->is_constant))
  {
    intervals = 1;

    delta = list[0].ds;

    curr = interval;

    curr->lit = TRUE;

    curr->samples = 0;

    curr->s0 = list[0].s0;
    curr->s1 = list[0].s0 + delta;
    curr->ds = delta;

    Make_Colour(curr->od,  0.0, 0.0, 0.0);
    Make_Colour(curr->te,  0.0, 0.0, 0.0);
    Make_Colour(curr->te2, 0.0, 0.0, 0.0);

    return(intervals);
  }

  /* Choose intervals. */

  if (number == 1)
  {
    /* Use uniform intervals. */

    delta = list[0].ds / (DBL)intervals;

    curr = interval;

    for (i = 0; i < intervals; i++)
    {
      curr->lit = TRUE;

      curr->samples = 0;

      curr->s0 = list[0].s0 + delta * (DBL)i;
      curr->s1 = list[0].s0 + delta * (DBL)(i + 1);
      curr->ds = delta;

      Make_Colour(curr->od,  0.0, 0.0, 0.0);
      Make_Colour(curr->te,  0.0, 0.0, 0.0);
      Make_Colour(curr->te2, 0.0, 0.0, 0.0);

      curr++;
    }
  }
  else
  {
    /* Choose intervals according to the specified ratio. */

    if (number > intervals)
    {
      Error("Too few sampling intervals.\n");
    }

    sum = 0.0;

    entry = list;

    for (i = 0; i < number; i++)
    {
/*
      sum += ((entry->lit) ? (0.9) : (0.1)) * entry->ds;
*/
      sum += ((entry->lit) ? (media->Ratio) : (1.0 - media->Ratio));

      entry++;
    }

    remaining = intervals;

    curr = interval;

    entry = list;

    for (i = 0; i < number; i++)
    {
/*
      weight = ((entry->lit) ? (0.9) : (0.1)) * entry->ds;
*/
      weight = ((entry->lit) ? (media->Ratio) : (1.0 - media->Ratio));

      n = (int)(weight / sum * (DBL)intervals) + 1;

      r = remaining - number + i + 1;

      if (n > r)
      {
        n = r;
      }

      delta = entry->ds / (DBL)n;

      for (j = 0; j < n; j++)
      {
        curr->lit = entry->lit;

        curr->samples = 0;

        curr->s0 = entry->s0 + delta * (DBL)j;
        curr->s1 = entry->s0 + delta * (DBL)(j + 1);
        curr->ds = delta;

        Make_Colour(curr->od,  0.0, 0.0, 0.0);
        Make_Colour(curr->te,  0.0, 0.0, 0.0);
        Make_Colour(curr->te2, 0.0, 0.0, 0.0);

        curr++;
      }

      remaining -= n;

      entry++;
    }
  }

  return(intervals);
}



/*****************************************************************************
*
* FUNCTION
*
*   intersect_spotlight
*
* INPUT
*
*   Ray    - current ray
*   Light  - current light source
*
* OUTPUT
*
*   d1, d2 - intersection depths
*
* RETURNS
*
*   int - TRUE, if hit
*
* AUTHOR
*
*   Dieter Bayer
*
* DESCRIPTION
*
*   Intersect a ray with the light cone of a spotlight.
*
* CHANGES
*
*   Nov 1994 : Creation.
*
******************************************************************************/

static int intersect_spotlight(RAY *Ray, LIGHT_SOURCE *Light, DBL *d1, DBL  *d2)
{
  int viewpoint_is_in_cone;
  DBL a, b, c, d, m, l, l1, l2, t, t1, t2, k1, k2, k3, k4;
  VECTOR V1;

  /* Get cone's slope. Note that cos(falloff) is stored in Falloff! */

  a = acos(Light->Falloff);

  /* This only works for a < 180 degrees! */

  m = tan(a);

  m = 1.0 + Sqr(m);

  VSub(V1, Ray->Initial, Light->Center);

  VDot(k1, Ray->Direction, Light->Direction);

  VDot(k2, V1, Light->Direction);

  VLength(l, V1);

  if (l > EPSILON)
  {
    viewpoint_is_in_cone = (k2 / l >= Light->Falloff);
  }
  else
  {
    viewpoint_is_in_cone = FALSE;
  }

  if ((k1 <= 0.0) && (k2 < 0.0))
  {
    return (FALSE);
  }

  VDot(k3, V1, Ray->Direction);

  VDot(k4, V1, V1);

  a = 1.0 - Sqr(k1) * m;

  b = k3 - k1 * k2 * m;

  c = k4 - Sqr(k2) * m;

  if (a != 0.0)
  {
    d = Sqr(b) - a * c;

    if (d > EPSILON)
    {
      d = sqrt(d);

      t1 = (-b + d) / a;
      t2 = (-b - d) / a;

      if (t1 > t2)
      {
        t = t1; t1 = t2; t2 = t;
      }

      l1 = k2 + t1 * k1;
      l2 = k2 + t2 * k1;

      if ((l1 <= 0.0) && (l2 <= 0.0))
      {
        return (FALSE);
      }

      if ((l1 <= 0.0) || (l2 <= 0.0))
      {
        if (l1 <= 0.0)
        {
          if (viewpoint_is_in_cone)
          {
            t1 = 0.0;
            t2 = (t2 > 0.0) ? (t2) : (Max_Distance);
          }
          else
          {
            t1 = t2;
            t2 = Max_Distance;
          }
        }
        else
        {
          if (viewpoint_is_in_cone)
          {
            t2 = t1;
            t1 = 0.0;
          }
          else
          {
            t2 = Max_Distance;
          }
        }
      }

      *d1 = t1;
      *d2 = t2;

      return (TRUE);
    }
    else
    {
      if (d > -EPSILON)
      {
        if (viewpoint_is_in_cone)
        {
          *d1 = 0.0;
          *d2 = -b / a;
        }
        else
        {
          *d1 = -b / a;
          *d2 = Max_Distance;
        }

        return(TRUE);
      }
    }
  }
  else
  {
    if (viewpoint_is_in_cone)
    {
      *d1 = 0.0;
      *d2 = -c/b;

      return(TRUE);
    }
  }

  return (FALSE);
}



/*****************************************************************************
*
* FUNCTION
*
*   intersect_cylinderlight
*
* INPUT
*
*   Ray    - current ray
*   Light  - current light source
*
* OUTPUT
*
*   d1, d2 - intersection depths
*
* RETURNS
*
*   int - TRUE, if hit
*
* AUTHOR
*
*   Dieter Bayer
*
* DESCRIPTION
*
*   Intersect a ray with the light cylinder of a cylinderlight.
*
* CHANGES
*
*   Jan 1995 : Creation.
*
******************************************************************************/

static int intersect_cylinderlight(RAY *Ray, LIGHT_SOURCE *Light, DBL *d1, DBL  *d2)
{
  DBL a, b, c, d, l1, l2, t, t1, t2, k1, k2, k3, k4;
  VECTOR V1;

  VSub(V1, Ray->Initial, Light->Center);

  VDot(k1, Ray->Direction, Light->Direction);

  VDot(k2, V1, Light->Direction);

  if ((k1 <= 0.0) && (k2 < 0.0))
  {
    return (FALSE);
  }

  a = 1.0 - Sqr(k1);

  if (a != 0.0)
  {
    VDot(k3, V1, Ray->Direction);

    VDot(k4, V1, V1);

    b = k3 - k1 * k2;

    c = k4 - Sqr(k2) - Sqr(Light->Falloff);

    d = Sqr(b) - a * c;

    if (d > EPSILON)
    {
      d = sqrt(d);

      t1 = (-b + d) / a;
      t2 = (-b - d) / a;

      if (t1 > t2)
      {
        t = t1; t1 = t2; t2 = t;
      }

      l1 = k2 + t1 * k1;
      l2 = k2 + t2 * k1;

      if ((l1 <= 0.0) && (l2 <= 0.0))
      {
        return (FALSE);
      }

      if ((l1 <= 0.0) || (l2 <= 0.0))
      {
        if (l1 <= 0.0)
        {
          t1 = 0.0;
        }
        else
        {
          t2 = (Max_Distance - k2) / k1;
        }
      }

      *d1 = t1;
      *d2 = t2;

      return (TRUE);
    }
  }

  return (FALSE);
}



/*****************************************************************************
*
* FUNCTION
*
*   Post_Media
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
*   Dieter Bayer
*
* DESCRIPTION
*
*   Initialize media.
*
* CHANGES
*
*   Dec 1996 : Creation.
*
******************************************************************************/

void Post_Media(IMEDIA *IMedia)
{
  int i;
  DBL t;

  if (IMedia == NULL)
  {
    return;
  }
  
  /* Get extinction coefficient. */

  CLinComb2(IMedia->Extinction, 1.0, IMedia->Absorption, IMedia->sc_ext, IMedia->Scattering);

  /* Determine used effects. */

  IMedia->use_absorption = (IMedia->Absorption[0] != 0.0) ||
                          (IMedia->Absorption[1] != 0.0) ||
                          (IMedia->Absorption[2] != 0.0);

  IMedia->use_emission = (IMedia->Emission[0] != 0.0) ||
                        (IMedia->Emission[1] != 0.0) ||
                        (IMedia->Emission[2] != 0.0);

  IMedia->use_scattering = (IMedia->Scattering[0] != 0.0) ||
                          (IMedia->Scattering[1] != 0.0) ||
                          (IMedia->Scattering[2] != 0.0);

  IMedia->use_extinction = IMedia->use_absorption || IMedia->use_scattering;

  IMedia->is_constant = (IMedia->Density == NULL);

  /* Init sample threshold array. */

  if (IMedia->Sample_Threshold != NULL)
  {
    POV_FREE(IMedia->Sample_Threshold);
  }

  /* Create list of thresholds for confidence test. */

  IMedia->Sample_Threshold = (DBL *)POV_MALLOC(IMedia->Max_Samples*sizeof(DBL), "sample threshold list");

  if (IMedia->Max_Samples > 1)
  {
    t = chdtri((DBL)(IMedia->Max_Samples-1), IMedia->Confidence);

    if (t > 0.0)
    {
      t = IMedia->Variance / t;
    }
    else
    {
      t = IMedia->Variance * EPSILON;
    }

    for (i = 0; i < IMedia->Max_Samples; i++)
    {
      IMedia->Sample_Threshold[i] = t * chdtri((DBL)(i+1), IMedia->Confidence);

/*
      fprintf(stderr, "threshold for n = %3d: %f\n", i+1, IMedia->Sample_Threshold[i]);
*/
    }
  }
  else
  {
    IMedia->Sample_Threshold[0] = 0.0;
  }

  if (IMedia->Density != NULL)
  {
    Post_Pigment(IMedia->Density);
  }
  
  Post_Media(IMedia->Next_Media);  
}



/*****************************************************************************
*
* FUNCTION
*
*   Create_Media
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
*   IMEDIA * - created media
*
* AUTHOR
*
*   Dieter Bayer
*
* DESCRIPTION
*
*   Create a media.
*
* CHANGES
*
*   Dec 1994 : Creation.
*
******************************************************************************/

IMEDIA *Create_Media()
{
  IMEDIA *New;

  New = (IMEDIA *)POV_MALLOC(sizeof(IMEDIA), "media");

  New->Type = ISOTROPIC_SCATTERING;

  New->Intervals    = 10;
  New->Min_Samples  = 1;
  New->Max_Samples  = 1;
  New->Eccentricity = 0.0;

  Make_Colour(New->Absorption, 0.0, 0.0, 0.0);
  Make_Colour(New->Emission,   0.0, 0.0, 0.0);
  Make_Colour(New->Extinction, 0.0, 0.0, 0.0);
  Make_Colour(New->Scattering, 0.0, 0.0, 0.0);

  New->use_absorption = FALSE;
  New->use_emission   = FALSE;
  New->use_extinction = FALSE;
  New->use_scattering = FALSE;

  New->is_constant = FALSE;

  New->sc_ext     = 1.0;
  New->Ratio      = 0.9;
  New->Confidence = 0.9;
  New->Variance   = 1.0 / 128.0;

  New->Sample_Threshold = NULL;

  New->Density = NULL;

  New->Next_Media = NULL;

  return (New);
}



/*****************************************************************************
*
* FUNCTION
*
*   Copy_Media
*
* INPUT
*
*   Old - media to copy
*
* OUTPUT
*
* RETURNS
*
*   IMEDIA * - new media
*
* AUTHOR
*
*   Dieter Bayer
*
* DESCRIPTION
*
*   Copy an media.
*
* CHANGES
*
*   Dec 1994 : Creation.
*
******************************************************************************/

IMEDIA *Copy_Media(IMEDIA *Old)
{
  int i;
  IMEDIA *New, *First, *Previous, *Local_Media;

  Previous = First = NULL;

  if (Old != NULL)
  {
    for (Local_Media = Old; Local_Media != NULL; Local_Media = Local_Media->Next_Media)
    {
      New = Create_Media();

      *New = *Local_Media;

      if (Local_Media->Sample_Threshold != NULL)
      {
        if (New->Intervals > 0)
        {
          New->Sample_Threshold = (DBL *)POV_MALLOC(New->Intervals*sizeof(DBL), "sample threshold list");

          for (i = 0; i < New->Intervals; i++)
          {
            New->Sample_Threshold[i] =  Local_Media->Sample_Threshold[i];
          }
        }
      }

      New->Density = Copy_Pigment(Local_Media->Density);

      if (First == NULL)
      {
        First = New;
      }

      if (Previous != NULL)
      {
        Previous->Next_Media = New;
      }

      Previous = New;
    }
  }

  return(First);
}



/*****************************************************************************
*
* FUNCTION
*
*   Destroy_Media
*
* INPUT
*
*   IMedia - media to destroy
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
*   Dieter Bayer
*
* DESCRIPTION
*
*   Destroy an media.
*
* CHANGES
*
*   Dec 1994 : Creation.
*
******************************************************************************/

void Destroy_Media(IMEDIA *IMedia)
{
  IMEDIA *Local_Media, *Temp;

  if (IMedia != NULL)
  {
    Local_Media = IMedia;

    while (Local_Media != NULL)
    {
      if (Local_Media->Sample_Threshold != NULL)
      {
        POV_FREE(Local_Media->Sample_Threshold);
      }

      /* Note Destroy_Pigment also handles Density->Next */
      Destroy_Pigment(Local_Media->Density);

      Temp = Local_Media->Next_Media;

      POV_FREE(Local_Media);

      Local_Media = Temp;
    }
  }
}



/*****************************************************************************
*
* FUNCTION
*
*   compdoubles
*
* INPUT
*
*   in_a, in_b - Elements to compare
*
* OUTPUT
*
* RETURNS
*
*   int - result of comparison
*
* AUTHOR
*
*   Dieter Bayer
*
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   Dec 1996 : Creation.
*
******************************************************************************/

static int CDECL compdoubles(CONST void *in_a, CONST void *in_b)
{
  DBL *a, *b;

  a = (DBL *)in_a;
  b = (DBL *)in_b;

  if (*a < *b)
  {
    return (-1);
  }
  else
  {
    if (*a == *b)
    {
      return (0);
    }
    else
    {
      return (1);
    }
  }
}



/*****************************************************************************
*
* FUNCTION
*
*   Transform_Media
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
*   Dieter Bayer
*
* DESCRIPTION
*
*   Transform media.
*
* CHANGES
*
*   Dec 1996 : Creation.
*
******************************************************************************/

void Transform_Media(IMEDIA *IMedia, TRANSFORM *Trans)
{
  IMEDIA *Temp;

  if (IMedia != NULL)
  {
    for (Temp = IMedia; Temp != NULL; Temp = Temp->Next_Media)
    {
      Transform_Density(Temp->Density, Trans);
    }
  }
}

void Transform_Density(PIGMENT *Density, TRANSFORM *Trans)
{
  TPATTERN *Temp = (TPATTERN *)Density;
  
  while (Temp != NULL)
  {
    Transform_Tpattern(Temp, Trans);
    Temp = Temp->Next;
  }
}