/* ************************************************************************* **
** AUTHOR: Alexis Angelidis						     **
** STARTED: Thu Jul 17 20:29:44 NZST 2003				     **
** ************************************************************************* */
#include <iostream>
#include <math.h>

/* ************************************************************************* **
** DESCRIPTION: Basic vertex class.					     **
** ************************************************************************* */
class Vertex3d
{
  double _m_x, _m_y, _m_z;
public:
  Vertex3d(const double x, const double y, const double z):
    _m_x(x), _m_y(y), _m_z(z) {}
  Vertex3d	operator+(const Vertex3d &v) const
  { return Vertex3d(_m_x + v._m_x, _m_y + v._m_y, _m_z + v._m_z); }
  Vertex3d	operator-(const Vertex3d &v) const
  { return Vertex3d(_m_x - v._m_x, _m_y - v._m_y, _m_z - v._m_z); }
  Vertex3d	operator/(const double d) const
  { return Vertex3d(_m_x / d, _m_y / d, _m_z / d); }
  double	operator*(const Vertex3d &v) const
  { return _m_x * v._m_x + _m_y * v._m_y + _m_z * v._m_z; }
  double dist(const Vertex3d &v) const { return sqrt(dist2(v)); }
  double dist2(const Vertex3d &v) const { return (*this - v).length2(); }
  double length2() const { return _m_x * _m_x + _m_y * _m_y + _m_z * _m_z; }
  friend Vertex3d	operator*(const double d, const Vertex3d &v)
  { return Vertex3d(d * v._m_x, d * v._m_y, d * v._m_z); }
  friend std::ostream &operator<<(std::ostream &os, const Vertex3d &v)
  { return os << v._m_x << " " << v._m_y << " " << v._m_z;}
};

/* ************************************************************************* **
** DESCRIPTION: Computes the field by convolution of r^2/d^2 along the       **
**   segment (v0, v1), where r = r0 at v0 and r = r1 at v1.		     **
* ************************************************************************* */
double	field(const Vertex3d	&p,
	      const Vertex3d	&v0,
	      const Vertex3d	&v1,
	      double		r0,
	      double		r1)
{
  const Vertex3d h(v0 + (((p - v0) * (v1 - v0)) / v1.dist2(v0)) * (v1 - v0));
  const double	 d2(h.dist2(p)), d(sqrt(d2)), invd(1.0 / d);
  const double	 d0(h.dist(v0)), d1(h.dist(v1));
  double	 a0, a1, rswap;;

  if (d0 < d1)
    {
      a1 = d1;
      a0 = ((v1 - h) * (v0 - h) < 0.0)? -d0: d0;
    }
  else
    {
      a1 = d0;
      a0 = ((v1 - h) * (v0 - h) < 0.0)? -d1: d1;
      rswap = r1;
      r1 = r0;
      r0 = rswap;
    }
  return (invd * (a1 * r0 - a0 * r1) * (atan(a0 * invd) - atan(a1 * invd))
	  + 0.5 * (r0 - r1) * log((a1 * a1 + d2)/(a0 * a0 + d2))) / (a0 - a1);
}

/* ************************************************************************* */
int	main(int	ac,
	     char	**av)
{
  std::cout << field(Vertex3d(0.5, 1.0, 1.0),
		     Vertex3d(0.0, 0.0, 0.0),
		     Vertex3d(1.0, 0.0, 0.0),
		     0.10, 0.10) << std::endl;
  return 0;
}