/* xxdf.c : some general functions regarding edf,brdf,btdf,bsdf */

#include "error.h"
#include "xxdf.h"
#include "vector.h"

/* Compute an approximate Geometric IOR from a complex IOR (cfr. Gr.Gems II, p289) */

float ComplexToGeometricRefractionIndex(REFRACTIONINDEX nc)
{
  float sqrtF, f1,f2;

  f1 = (nc.nr - 1.0); 
  f1 = f1*f1 + nc.ni * nc.ni;

  f2 = (nc.nr + 1.0); 
  f2 = f2*f2 + nc.ni * nc.ni;

  sqrtF = sqrt(f1 / f2);

  return (1.0 + sqrtF) / (1.0 - sqrtF);
}

/* Calculate the ideal reflected ray direction (independent of the brdf) */
VECTOR IdealReflectedDirection(VECTOR *in, VECTOR *normal)
{
  double tmp;
  VECTOR result;

  tmp = 2 * VECTORDOTPRODUCT(*normal, *in);
  VECTORSCALE(tmp, *normal, result);
  VECTORSUBTRACT(*in, result, result);
  VECTORNORMALISE(result);

  return result;
}

DVECTOR DIdealReflectedDirection(DVECTOR *in, DVECTOR *normal)
{
  double tmp;
  DVECTOR result;

  tmp = 2.0 * VECTORDOTPRODUCT(*normal, *in);
  VECTORSCALE(tmp, *normal, result);
  VECTORSUBTRACT(*in, result, result);
  VECTORNORMALISE(result);

  return result;
}



/* Calculate the perfect refracted ray direction.
 * Sets totalInternalReflection to TRUE or FALSE accordingly.
 * Cfr. An Introduction to RayTracing (Glassner)
 */

VECTOR IdealRefractedDirection(VECTOR *in, VECTOR *normal, 
				   REFRACTIONINDEX inIndex, 
				   REFRACTIONINDEX outIndex, 
				   int *totalInternalReflection)
{
  double Ci,Ct, Ct2;
  double refractionIndex;
  /*  float nr, ni; */
  double normalScale;
  VECTOR result;

  /* !! Only real part of n for now */

  refractionIndex = inIndex.nr / outIndex.nr;

  Ci = - VECTORDOTPRODUCT(*in, *normal);
  Ct2 = 1 + refractionIndex * refractionIndex * ( Ci * Ci - 1);

  if(Ct2 < 0)
  {
    *totalInternalReflection = TRUE;
    return(IdealReflectedDirection(in, normal));
  }
  *totalInternalReflection = FALSE;

  Ct = sqrt(Ct2);

  normalScale = refractionIndex * Ci - Ct;

  VECTORSCALE(refractionIndex, *in, result);
  VECTORADDSCALED(result, normalScale, *normal, result);

  /* Normalise */
  VECTORNORMALISE(result);

  return result;
}


DVECTOR DIdealRefractedDirection(DVECTOR *in, DVECTOR *normal, 
				 REFRACTIONINDEX inIndex, 
				 REFRACTIONINDEX outIndex, 
				 int *totalInternalReflection)
{
  double Ci,Ct, Ct2;
  double refractionIndex;
  /* float nr, ni; */
  double normalScale;
  DVECTOR result;

  /* !! Only real part of n for now */
  refractionIndex = inIndex.nr / outIndex.nr;

  Ci = - VECTORDOTPRODUCT(*in, *normal);
  Ct2 = 1 + refractionIndex * refractionIndex * ( Ci * Ci - 1);

  if(Ct2 < 0)
  {
    *totalInternalReflection = TRUE;
    return(DIdealReflectedDirection(in, normal));
  }
  *totalInternalReflection = FALSE;

  Ct = sqrt(Ct2);

  normalScale = refractionIndex * Ci - Ct;

  VECTORSCALE(refractionIndex, *in, result);
  VECTORADDSCALED(result, normalScale, *normal, result);

  /* Normalise */
  VECTORNORMALISE(result);

  return result;
}