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

namespace RAYPP {

inline VECTOR TUBE::Get_Normal (const VECTOR &rawnorm) const
  {
  VECTOR tmp = Trans.TransNormal (rawnorm);
  tmp.Normalize();
  return tmp;
  }

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

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

//virtual
AXISBOX TUBE::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 TUBE::Test (const GEOM_RAY &Ray, float8 &dist, bool &realhit) const
  {
  ci();
  realhit = true;

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

  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.maxdist) return false;
    if (dist > Ray.mindist)
      {
      VECTOR tmpvec = s + dist*d;
      if ((tmpvec.y > 0) && (tmpvec.y < 1))
        return true;
      }
    dist = (-dxz_sxz + discrim)/(dxz_quad);
    if ((dist > Ray.mindist) && (dist < Ray.maxdist))
      {
      VECTOR tmpvec = s + dist*d;
      if ((tmpvec.y > 0) && (tmpvec.y < 1))
        return true;
      }
    }

  return false;
  }

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

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

  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.maxdist) return false;
    if (dist > Ray.mindist)
      {
      VECTOR tmpvec = s + dist*d;
      if ((tmpvec.y > 0) && (tmpvec.y < 1))
        {
        Normal = Get_Normal (VECTOR (tmpvec.x,0,tmpvec.z));
        return true;
        }
      }
    dist = (-dxz_sxz + discrim)/(dxz_quad);
    if ((dist > Ray.mindist) && (dist < Ray.maxdist))
      {
      VECTOR tmpvec = s + dist*d;
      if ((tmpvec.y > 0) && (tmpvec.y < 1))
        {
        Normal = Get_Normal (VECTOR (tmpvec.x,0,tmpvec.z));
        return true;
        }
      }
    }

  return false;
  }

//virtual
void TUBE::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;

  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 tmpvec = s + dist*d;
      if ((tmpvec.y > 0) && (tmpvec.y < 1))
        Inter.push_back
          (INTER (dist, Get_Normal (VECTOR (tmpvec.x,0,tmpvec.z))));
      }
    dist = (-dxz_sxz + discrim)/(dxz_quad);
    if ((dist > Ray.mindist) && (dist < Ray.maxdist))
      {
      VECTOR tmpvec = s + dist*d;
      if ((tmpvec.y > 0) && (tmpvec.y < 1))
        Inter.push_back
          (INTER (dist, Get_Normal (VECTOR (tmpvec.x,0,tmpvec.z))));
      }
    }
  }

} // namespace RAYPP