/*
 *  Ray++ - Object-oriented ray tracing library
 *  Copyright (C) 1998-2001 Martin Reinecke and others.
 *  See the AUTHORS file for more information.
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Library General Public
 *  License as published by the Free Software Foundation; either
 *  version 2 of the License, or (at your option) any later version.
 *
 *  This library is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  Library General Public License for more details.
 *
 *  You should have received a copy of the GNU Library General Public
 *  License along with this library; if not, write to the Free
 *  Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 *
 *  See the README file for more information.
 */

#include "volumes/gas2.h"

namespace RAYPP {

extern HANDLE<RENDERER> Renderer;

GAS2::GAS2 ()
  : phase(RAYLEIGH),
    min_step(float4(0.01)), max_step(float4(0.2)),
    abs_tol(float4(1.0/256)), rel_tol(float4(1.0/256))
  {}

GAS2::GAS2 (phase_func ph, const HANDLE<OPACITY_FUNCTOR> &Opac)
  : phase(ph), Opacfunc (Opac),
    min_step(float4(0.01)), max_step(float4(0.2)),
    abs_tol(float4(1.0/256)), rel_tol(float4(1.0/256))
  {}

COLOUR GAS2::Integrate_Opacity_Function(const VECTOR &Lower,
  const VECTOR &Higher, const COLOUR &Ingoing) const
  {
  VECTOR diff = Higher - Lower;
  float8 step = diff.Length();
  if (step < Small_float4) return COLOUR(0,0,0);

  diff.Normalize();

  const float8 Pgrow   = -0.33;
  const float8 Pshrink = -0.5;
  const float8 Fcor    = 1.0/3.0;
  const float8 Safety  = 0.9;
  const float8 Errcon  = exp((1.0/Pgrow)*log(4.0/Safety));

  VECTOR x    = Higher;
  float8 xpos = 0.0;
  float4 h    = min_step;
  float4 hnext;

  COLOUR y (Ingoing);
  COLOUR dynew = -y*Opacfunc->Get_Opacity(x);

  while (true)
    {
    COLOUR dyl = dynew;
    COLOUR dy1, dy2, delta;
    if ((xpos+h) > step) h = step-xpos;

    while (true)
      {
      float4 hh = 0.5*h;
      VECTOR xr = x - h*diff;
      COLOUR dyh = -(y-dyl*h)*Opacfunc->Get_Opacity(xr);
      COLOUR dym = -(y-dyl*hh)*Opacfunc->Get_Opacity(x-hh*diff);
      COLOUR dyh2 = -(y-(dyl+dym)*hh)*Opacfunc->Get_Opacity(xr);
      dy1 = hh*(dyl+dyh);
      dy2 = hh*(0.5*(dyl+dyh2) + dym);
      delta = dy2-dy1;
      float8 err = delta.AbsoluteSum()/max(abs_tol, rel_tol*y.AbsoluteSum());
      if (err <= 1.0)
        {
        hnext = ((err>Errcon) ?
          (Safety*h*exp(Pgrow*log(err))) : (4.0*h));
        hnext = max (min_step, min(max_step, hnext));
        dynew = dyh;
        break;
        }
      if (h < (min_step + Small_float4))
        {
        hnext = min_step;
        dynew = dyh;
        break;
        }
      hnext = Safety*h*exp(Pshrink*log(err));
      hnext = max (min_step, min(max_step, hnext));
      h = hnext;
      }

    y += dy2 + Fcor*delta;
    x -= h*diff;
    xpos += h;
    h = hnext;
    if (abs(xpos-step) < Small_float8)
      return y;
    }
  }

//virtual
void GAS2::Init ()
  {
  if (initialized) return;
  if (!Opacfunc)
    error ("GAS2: no opacity function given!");
  initialized = true;
  }

//virtual
void GAS2::Transform (const TRANSFORM &trans)
  {
  cni();
  Trans.Add_Transform (trans);
  }

//virtual 
COLOUR GAS2::Calc_new_Importance (const RAY &Ray) const
  {
  ci();

  VECTOR Near = Trans.InvTransPoint (Ray.evalmin());
  VECTOR Far  = Trans.InvTransPoint (Ray.evalmax());

  return Integrate_Opacity_Function(Near,Far,Ray.Importance);
  }

//private
COLOUR GAS2::Phase_Function(const float4 cosangle) const
  {
  switch (phase)
    {
    case LAMBERT:   //Lambert Surface Phase Function
      {      
      const float8 a = acos(cosangle);
      return COLOUR(1,1,1) * ( 8/(3*Pi) * (sin(a) + (Pi - a) * cosangle));
      }

    case ANISOTROPIC: //Anisotropic
      {
// integration over all angles must be <= 1
      const float4 x = float4(0.8);
      return COLOUR(1,1,1) * (1 + x * cosangle);
      }

    case RAYLEIGH: //Rayleigh Scattering
      return COLOUR(1,1,1)*( 0.75*(1+cosangle*cosangle));

    default:
      return COLOUR(1,1,1);
    }
  }

COLOUR GAS2::Get_Contrib (const COLOUR &ActInt, const VECTOR &Interest,
  const VECTOR &dir) const
  {
  COLOUR term2(0,0,0);
  COLOUR term3 = Opacfunc->Get_Opacity(Interest);
  INCIDENT_ARRAY Arr;

  Renderer->Calc_Illumination (Interest, COLOUR(1,1,1), Arr);
  for(uint4 i=0; i<Arr.size(); ++i)
    {
    float8 cosangle = Dot(Arr[i].Direction,dir);
    term2 += Arr[i].Intensity * Phase_Function(cosangle);
    }
  term2 += Arr.Ambient;

  return  term3*(-ActInt + term2);
  }

//virtual
COLOUR GAS2::Calc_Modified_Colour (const RAY &Ray, const COLOUR &Ingoing) const
  {
  ci();

  VECTOR Lower  = Trans.InvTransPoint (Ray.evalmin());
  VECTOR Higher = Trans.InvTransPoint (Ray.evalmax());
  VECTOR dir    = Ray.dir;
  VECTOR diff = Higher - Lower;
  float8 step = diff.Length();
  if (step < Small_float4) return COLOUR(0,0,0);

  diff.Normalize();

  const float8 Pgrow   = -0.33;
  const float8 Pshrink = -0.5;
  const float8 Fcor    = 1.0/3.0;
  const float8 Safety  = 0.9;
  const float8 Errcon  = exp((1.0/Pgrow)*log(4.0/Safety));

  VECTOR x    = Higher;
  float8 xpos = 0.0;
  float4 h    = min_step;
  float4 hnext;

  COLOUR y (Ingoing);
  COLOUR dynew=Get_Contrib(y,x,dir);

  while (true)
    {
    COLOUR dyl = dynew;
    COLOUR dy1, dy2, delta;
    if ((xpos+h) > step) h = step-xpos;

    while (true)
      {
      float4 hh = 0.5*h;
      VECTOR xr = x - h*diff;
      COLOUR dyh  = Get_Contrib (y+dyl*h, xr, dir);
      COLOUR dym  = Get_Contrib (y+dyl*hh, x-hh*diff, dir);
      COLOUR dyh2 = Get_Contrib (y+(dyl+dym)*hh, xr, dir);
      dy1 = hh*(dyl+dyh);
      dy2 = hh*(0.5*(dyl+dyh2) + dym);
      delta = dy2-dy1;
      float8 err = delta.AbsoluteSum()/max(abs_tol, rel_tol*y.AbsoluteSum());
      if (err <= 1.0)
        {
        hnext = ((err>Errcon) ?
          (Safety*h*exp(Pgrow*log(err))) : (4.0*h));
        hnext = max (min_step, min(max_step, hnext));
        dynew = dyh2;
        break;
        }
      if (h < (min_step+Small_float4))
        {
        hnext = min_step;
        dynew = dyh2;
        break;
        }
      hnext = Safety*h*exp(Pshrink*log(err));
      hnext = max (min_step, min(max_step, hnext));

      h = hnext;
      }

    y += dy2 + Fcor*delta;
    x -= h*diff;
    xpos += h;
    h = hnext;
    if (abs(xpos-step) < Small_float8)
      return y;
    }
  }

//virtual 
COLOUR GAS2::Get_Attenuated_Light (const RAY &Ray, const COLOUR &Ingoing) const
  {
  ci ();

  VECTOR Near = Trans.InvTransPoint (Ray.evalmin());
  VECTOR Far  = Trans.InvTransPoint (Ray.evalmax());

  return Integrate_Opacity_Function(Near,Far, Ingoing);
  }

void GAS2::Set_Phase (phase_func pf)
  {
  cni();
  phase = pf;
  }

void GAS2::Set_Opacity  (const HANDLE<OPACITY_FUNCTOR> &Opac)
  {
  cni();
  Opacfunc = Opac;
  }

void GAS2::Set_Constants (float4 min_st, float4 max_st,
                         float4 abstol, float4 reltol)
  {
  cni();
  min_step = min_st;
  max_step = max_st;
  abs_tol  = abstol;
  rel_tol  = reltol;
  }

} // namespace RAYPP