/*
 *  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 "shapes/cone.h"

namespace RAYPP {

const float8 CONE::Epsilon = 1e-7;

inline VECTOR CONE::Get_Normal (const VECTOR &rawnorm) const
  {
  VECTOR tmp = Trans.TransNormal (rawnorm);
  tmp.Normalize();
  if (Inverted) tmp.Flip();
  return tmp;
  }

CONE::CONE()
  : Trans(Translation_Transform(VECTOR(0,1,0)))
  {}

//virtual
void CONE::Init ()
  {
  initialized = true;
  }

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

//virtual
AXISBOX CONE::BBox () const
  {
  ci();
  TRANSFORM t2 = Trans;
  VECTOR lcnt = t2.TransPoint (VECTOR(0,-1,0));
  VECTOR xdir = t2.TransDirection (VECTOR(1,0,0));
  VECTOR zdir = t2.TransDirection (VECTOR(0,0,1));

  float8 phi_x = atan2 (zdir.x, xdir.x);
  float8 phi_y = atan2 (zdir.y, xdir.y);
  float8 phi_z = atan2 (zdir.z, xdir.z);
  VECTOR offset (abs (xdir.x*cos(phi_x) + zdir.x*sin(phi_x)),
                 abs (xdir.y*cos(phi_y) + zdir.y*sin(phi_y)),
                 abs (xdir.z*cos(phi_z) + zdir.z*sin(phi_z)));

  AXISBOX Box (lcnt-offset, lcnt+offset);
  Box.Include (t2.TransPoint(VECTOR(0,0,0)));

  return Box;
  }

//virtual
bool CONE::Test (const GEOM_RAY &Ray, float8 &dist, bool &realhit) const
  {
  ci();
  const VECTOR s = Trans.InvTransPoint (Ray.start),
               d = Trans.InvTransDirection (Ray.dir);
  float8 maxd = Ray.maxdist;
  bool found = false;

  if (abs(d.y) > Epsilon)
    {
    dist = (-1-s.y)/d.y;
    if ((dist > Ray.mindist) && (dist < maxd))
      {
      const float8 tmp1 = s.x + dist*d.x;
      const float8 tmp2 = s.z + dist*d.z;
      if ((tmp1*tmp1 + tmp2*tmp2) < 1)
        {
        maxd = dist;
        found = true;
        }
      }
    }

  const float8 dxyz_quad = d.x*d.x + d.z*d.z - d.y*d.y;
  if (abs(dxyz_quad) > Epsilon)
    {
    const float8 sxyz_quad = s.x*s.x + s.z*s.z - s.y*s.y,
                 dxyz_sxyz = s.x*d.x + s.z*d.z - s.y*d.y;

    float8 discrim = dxyz_sxyz*dxyz_sxyz - dxyz_quad*sxyz_quad;
    if (discrim <= Epsilon) return false;

    discrim = sqrt (discrim);
    dist = (-dxyz_sxyz - discrim)/dxyz_quad;
    if ((dist > Ray.mindist) && (dist < maxd))
      {
      VECTOR tmpvec = s + dist*d;
      if ((tmpvec.y > -1) && (tmpvec.y <= 0))
        {
        maxd = dist;
        found = true;
        }
      }
    dist = (-dxyz_sxyz + discrim)/dxyz_quad;
    if ((dist > Ray.mindist) && (dist < maxd))
      {
      VECTOR tmpvec = s + dist*d;
      if ((tmpvec.y > -1) && (tmpvec.y <= 0))
        {
        maxd = dist;
        found = true;
        }
      }
    }

  if (!found) return false;
  realhit = true;
  dist = maxd;
  return true;
  }

//virtual
bool CONE::Inside (const VECTOR &Loc) const
  {
  ci();
  const VECTOR tmp = Trans.InvTransPoint (Loc);
  if ((tmp.y<(-1)) || (tmp.y>0)) return Inverted;
  return (((tmp.x*tmp.x + tmp.z*tmp.z - tmp.y*tmp.y) <= 0) ^ Inverted);
  }

//virtual
bool CONE::Intersect (const GEOM_RAY &Ray,
  float8 &dist, VECTOR &Normal) const
  {
  ci();

  const VECTOR s = Trans.InvTransPoint (Ray.start),
               d = Trans.InvTransDirection (Ray.dir);
  float8 maxd = Ray.maxdist;
  bool found = false;

  if (abs(d.y) > Epsilon)
    {
    dist = (-1-s.y)/d.y;
    if ((dist > Ray.mindist) && (dist < maxd))
      {
      const float8 tmp1 = s.x + dist*d.x;
      const float8 tmp2 = s.z + dist*d.z;
      if ((tmp1*tmp1 + tmp2*tmp2) < 1)
        {
        maxd = dist;
        found = true;
        Normal = VECTOR (0,-1,0);
        }
      }
    }

  const float8 dxyz_quad = d.x*d.x + d.z*d.z - d.y*d.y;
  if (abs(dxyz_quad) > Epsilon)
    {
    const float8 sxyz_quad = s.x*s.x + s.z*s.z -s.y*s.y,
                 dxyz_sxyz = s.x*d.x + s.z*d.z -s.y*d.y;

    float8 discrim = dxyz_sxyz*dxyz_sxyz - dxyz_quad*sxyz_quad;
    if (discrim <= Epsilon) return false;

    discrim = sqrt (discrim);
    dist = (-dxyz_sxyz - discrim)/dxyz_quad;
    if ((dist > Ray.mindist) && (dist < maxd))
      {
      VECTOR tmpvec = s + dist*d;
      if ((tmpvec.y > -1) && (tmpvec.y <= 0))
        {
        maxd = dist;
        found = true;
        Normal = VECTOR (tmpvec.x,-tmpvec.y,tmpvec.z);
        }
      }
    dist = (-dxyz_sxyz + discrim)/dxyz_quad;
    if ((dist > Ray.mindist) && (dist < maxd))
      {
      VECTOR tmpvec = s + dist*d;
      if ((tmpvec.y > -1) && (tmpvec.y <= 0))
        {
        maxd = dist;
        found = true;
        Normal = VECTOR (tmpvec.x,-tmpvec.y,tmpvec.z);
        }
      }
    }

  if (!found) return false;

  dist = maxd;
  Normal = Get_Normal (Normal);
  return true;
  }

//virtual
void CONE::All_Intersections (const GEOM_RAY &Ray,
  vector<INTER> &Inter) const
  {
  ci();

  const VECTOR s = Trans.InvTransPoint (Ray.start),
               d = Trans.InvTransDirection (Ray.dir);

  float8 dist;

  if (abs(d.y) > Epsilon)
    {
    dist = (-1-s.y)/d.y;
    if ((dist > Ray.mindist) && (dist < Ray.maxdist))
      {
      const float8 tmp1 = s.x + dist*d.x;
      const float8 tmp2 = s.z + dist*d.z;
      if ((tmp1*tmp1 + tmp2*tmp2) < 1)
        Inter.push_back (INTER (dist, Get_Normal (VECTOR (0,-1,0))));
      }
    }

  const float8 dxyz_quad = d.x*d.x + d.z*d.z - d.y*d.y;
  if (abs(dxyz_quad) > Epsilon)
    {
    const float8 sxyz_quad = s.x*s.x + s.z*s.z - s.y*s.y,
                 dxyz_sxyz = s.x*d.x + s.z*d.z - s.y*d.y;

    float8 discrim = dxyz_sxyz*dxyz_sxyz - dxyz_quad*sxyz_quad;
    if (discrim <= Epsilon) return;
    discrim = sqrt (discrim);
    dist = (-dxyz_sxyz - discrim)/dxyz_quad;
    if ((dist > Ray.mindist) && (dist < Ray.maxdist))
      {
      VECTOR tmpvec = s + dist*d;
      if ((tmpvec.y > -1) && (tmpvec.y <= 0))
        Inter.push_back
          (INTER (dist, Get_Normal (VECTOR (tmpvec.x,-tmpvec.y,tmpvec.z))));
      }
    dist = (-dxyz_sxyz + discrim)/dxyz_quad;
    if ((dist > Ray.mindist) && (dist < Ray.maxdist))
      {
      VECTOR tmpvec = s + dist*d;
      if ((tmpvec.y > -1) && (tmpvec.y <= 0))
        Inter.push_back
          (INTER (dist, Get_Normal (VECTOR (tmpvec.x,-tmpvec.y,tmpvec.z))));
      }
    }
  }

} // namespace RAYPP