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

namespace RAYPP {

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

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

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

//virtual
AXISBOX CYLINDER::BBox () const
  {
  ci();
  TRANSFORM t2 = Trans;
  VECTOR lcnt = t2.TransPoint (VECTOR(0,0,0));
  VECTOR ucnt = 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)));

  VECTOR Min = lcnt;
  Min.Minimize (ucnt);
  VECTOR Max = lcnt;
  Max.Maximize (ucnt);
  return AXISBOX (Min-offset, Max+offset);
  }

//virtual
bool CYLINDER::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) > Small_dist)
    {
    dist = -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;
        }
      }
    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 dxz_quad = d.x*d.x + d.z*d.z;
  if (dxz_quad > Small_dist*Small_dist)
    {
    const float8 sxz_quad = s.x*s.x + s.z*s.z,
                 dxz_sxz  = s.x*d.x + s.z*d.z;

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

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

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

//virtual
bool CYLINDER::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) > Small_dist)
    {
    dist = -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);
        }
      }
    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 dxz_quad = d.x*d.x + d.z*d.z;
  if (dxz_quad > Small_dist*Small_dist)
    {
    const float8 sxz_quad = s.x*s.x + s.z*s.z,
                 dxz_sxz  = s.x*d.x + s.z*d.z;

    float8 discrim = dxz_sxz*dxz_sxz - dxz_quad*(sxz_quad - 1);
    if (discrim <= Small_dist*Small_dist) return false;

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

  if (!found) return false;

  dist = maxd;
  Normal = Get_Normal (Normal);

  return true;
  }

//virtual
void CYLINDER::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) > Small_dist)
    {
    dist = -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))));
      }
    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 dxz_quad = d.x*d.x + d.z*d.z;
  if (dxz_quad > Small_dist*Small_dist)
    {
    const float8 sxz_quad = s.x*s.x + s.z*s.z,
                 dxz_sxz  = s.x*d.x + s.z*d.z;

    float8 discrim = dxz_sxz*dxz_sxz - dxz_quad*(sxz_quad - 1);
    if (discrim <= Small_dist*Small_dist) return;
    discrim = sqrt (discrim);
    dist = (-dxz_sxz - discrim)/(dxz_quad);
    if ((dist > Ray.mindist) && (dist < Ray.maxdist))
      {
      VECTOR tmp = s + dist*d;
      if ((tmp.y > 0) && (tmp.y < 1))
        Inter.push_back (INTER (dist, Get_Normal (VECTOR (tmp.x,0,tmp.z))));
      }
    dist = (-dxz_sxz + discrim)/(dxz_quad);
    if ((dist > Ray.mindist) && (dist < Ray.maxdist))
      {
      VECTOR tmp = s + dist*d;
      if ((tmp.y > 0) && (tmp.y < 1))
        Inter.push_back (INTER (dist, Get_Normal (VECTOR (tmp.x,0,tmp.z))));
      }
    }
  }

} // namespace RAYPP