a_point.cxx

Go to the documentation of this file.

/*
         Copyright 2010-17 Pierre SMARS (smars@yuntech.edu.tw)
         This program is free software: you can redistribute it and/or modify
         it under the terms of the GNU General Public License as published by
         the Free Software Foundation, either version 2 of the License, or
         (at your option) any later version.
 
         This program 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 General Public License for more details.
 
         You should have received a copy of the GNU General Public License
         along with this program.  If not, see <http://www.gnu.org/licenses/>.
         */
#include <math.h>
#include <sstream>
 
#include "a_point.h"
 
//---------------------------------------------------------------------------
const std::string a_point::help()
{
        std::ostringstream o;
        o << "********" << std::endl;
        o << "a_point:" << std::endl;
        o << "********" << std::endl;
        o << "This is a 3D 'point' or 'vector' class" << std::endl;
        o << std::endl;
        o << "Create a point:" << std::endl;
        o << " a_point a" << std::endl;
        o << " a_point b 1. 2.5 0." << std::endl;
        o << std::endl;
        o << "Commands:" << std::endl;
        o << "+, -: add or subtract points." << std::endl;
        o << "*: scalar product (or multiply by a scalar)." << std::endl;
        o << "cross: cross product." << std::endl;
        o << "/: divide by a scalar." << std::endl;
        o << "x, y, z, x _x, y _y, z _z: get or set coordinates." << std::endl;
        o << "set _x _y _z: set coordinates." << std::endl;
        o << "set_cylindrical _r _phi _h: set coordinates (using cylindrical coordinate input)." << std::endl;
        o << "set_spherical _r _theta _phi: set coordinates (using spherical coordinate input)." << std::endl;
        o << "clear: set to (0,0,0)." << std::endl;
        o << "translate _x _y _z: translate." << std::endl;
        o << "rotate _xaxis _yaxis _zaxis: rotate." << std::endl;
        o << "rotate _axis _angle: rotate around an axis of a given angle _angle." << std::endl;
        o << "sumsq: sum of the square of the coordinates." << std::endl;
        o << "norm: returns norm." << std::endl;
        o << "normalise: normalise point." << std::endl;
        o << "norm1: taxi cab norm." << std::endl;
        o << "norm2: euclidian norm." << std::endl;
        o << "normI: maximum norm." << std::endl;
        o << "dist _pt: returns distance to _pt." << std::endl;
        o << "max: set maximum abs value of component to 1." << std::endl;
        o << "_pt1 angle _pt2: returns angle between 2 points." << std::endl;
        o << "average _pt1 _pt2 _f: returns intermediate point; average coordinates; proportion is given by _f." << std::endl;
        o << "average_rot _pt1 _pt2 _pos: returns intermediate point; average norm and average orientation; proportion is given by _f." << std::endl;
        o << std::endl;
        o << "Examples:" << std::endl;
        o << " puts \"a: [a print]\"" << std::endl;
        o << " puts \"a-b: [ [a - b ] print]\"" << std::endl;   
        o << " puts \"a*b: [ a * b]\"" << std::endl;    
        o << " puts \"||a||: [a norm]\"" << std::endl;  
        o << " a normalise" << std::endl;       
        o << " puts \"distance ab: [a dist b]\"" << std::endl;  
        o << "..." << std::endl;        
        o << a_geom_base::help();
        return o.str();
}
//---------------------------------------------------------------------------
a_point& a_point::operator-()
{
        x_ = -x_;
        y_ = -y_;
        z_ = -z_;
        return *this;
}
//---------------------------------------------------------------------------
a_point a_point::operator+(const a_point& p)
{
        a_point p2 = *this;
        return p2 += p;
}
//---------------------------------------------------------------------------
a_point a_point::operator-(const a_point& p)
{
        a_point p2 = *this;
        return p2 -= p;
}
//---------------------------------------------------------------------------
double a_point::operator*(const a_point& p)
{
        return x_*p.x()+y_*p.y()+z_*p.z();
}
//---------------------------------------------------------------------------
a_point a_point::cross(const a_point& p) const
{
        a_point r;
        r.x(y_*p.z()-z_*p.y());
        r.y(z_*p.x()-x_*p.z());
        r.z(x_*p.y()-y_*p.x());
        return r;
}
//---------------------------------------------------------------------------
/*
         a_point a_point::operator*(double v)
         {
         a_point p2 = *this;
         return p2*v;
         }
         */
//---------------------------------------------------------------------------
a_point& a_point::operator+=(const a_point& p)
{
        x_ += p.x();
        y_ += p.y();
        z_ += p.z();
        return *this;
}
//---------------------------------------------------------------------------
a_point& a_point::operator-=(const a_point& p)
{
        x_ -= p.x();
        y_ -= p.y();
        z_ -= p.z();
        return *this;
}
//---------------------------------------------------------------------------
double a_point::operator*=(const a_point& p)
{
        return x_*p.x()+y_*p.y()+z_*p.z();
}
//---------------------------------------------------------------------------
bool a_point::operator==(const a_point& p) const
{
        return ((x_ == p.x())&&(y_ == p.y())&&(z_ == p.z()));
}
//---------------------------------------------------------------------------
a_point & a_point::operator=(const a_point& p)
{
        x_ = p.x();
        y_ = p.y();
        z_ = p.z();
        return *this;
}
//-operator*=----------------------------------------------------------------
a_point & a_point::operator*=(double v)
{
        x_ *= v;
        y_ *= v;
        z_ *= v;
        return *this;
}
//-operator/=----------------------------------------------------------------
a_point & a_point::operator/=(double v)
{
        x_ /= v;
        y_ /= v;
        z_ /= v;
        return *this;
}
//-sumsq---------------------------------------------------------------------
double a_point::sumsq() const
{
        return x_*x_+y_*y_+z_*z_;
}
//-norm----------------------------------------------------------------------
double a_point::norm() const
{
        return sqrt(sumsq());
}
double a_point::norm1() const
{
        return fabs(x_)+fabs(y_)+fabs(z_);
}
double a_point::norm2() const
{
        return sqrt(sumsq());
}
double a_point::normI() const
{
        double max(fabs(x_));
        if (fabs(y_) > max)
                max = fabs(y_);
        if (fabs(z_) > max)
                max = fabs(z_);
        return max;
}
//-normalise-----------------------------------------------------------------
a_point& a_point::normalise()
{
        double n = norm();
        if (n != 0) operator/=(n);
        return *this;
}
//---------------------------------------------------------------------------
a_mat_sq a_point::inertia() const
{
        a_mat_sq a(3);
        for (short i=0; i<3; i++)
        { 
                double xi = (*this)[i];
                for (short j=0; j<=i; j++)
                { 
                        double xj = (*this)[j];
                        double v = xi*xj;
                        a(i,j) = v;
                        a(j,i) = v;
                }
        }
        return a;
}
//---------------------------------------------------------------------------
a_point& a_point::max()
{
        double max = fabs(x_);
        if (fabs(y_)>max)
                max = fabs(y_);
        if (fabs(z_)>max)
                max = fabs(z_);
        if (max != 0) operator/=(max);
        return *this;
}
//-rotate--------------------------------------------------------------------
a_point& a_point::rotate(const a_point& x_axis, const a_point& y_axis, const a_point& z_axis)
{
        double xn = (*this)*x_axis;
        double yn = (*this)*y_axis;
        double zn = (*this)*z_axis;
        x_ = xn; y_ = yn; z_ = zn;
        return *this;
}
//-rotate--------------------------------------------------------------------
a_point& a_point::rotate(const a_point& axis, const double angle)
{
        const a_point r = *this;
        a_point v = (r*axis)*axis;
        a_point x = r-v;
        a_point y = axis.cross(r);
        a_point rn = v+cos(angle)*x+sin(angle)*y;
        x_ = rn.x();
        y_ = rn.y();
        z_ = rn.z();
        return *this;
}
//***************************************************************************
//-operator+-----------------------------------------------------------------
a_point operator+(const a_point& a, const a_point& b)
{
        a_point r = a;
        return r += b;
}
//-operator------------------------------------------------------------------
a_point operator-(const a_point& a, const a_point& b)
{
        a_point r = a;
        return r -= b;
}
//-operator*-----------------------------------------------------------------
double operator*(const a_point& a, const a_point& b)
{
        a_point r = a;
        return r *= b;
}
//-operator*-----------------------------------------------------------------
a_point operator*(const double v, const a_point& a)
{
        a_point r = a;
        return r *= v;
}
//-operator*-----------------------------------------------------------------
a_point operator*(const a_point& a, double v)
{
        a_point r = a;
        return r *= v;
}
//-operator/-----------------------------------------------------------------
a_point operator/(const a_point& a, double v)
{
        a_point r = a;
        return r /= v;
}
//-cross---------------------------------------------------------------------
a_point cross(const a_point& a, const a_point& b)
{
        a_point r;
        r.x(a.y()*b.z()-a.z()*b.y());
        r.y(a.z()*b.x()-a.x()*b.z());
        r.z(a.x()*b.y()-a.y()*b.x());
        return r;
}
//---------------------------------------------------------------------------
a_point circle_centre(const a_point& p1, const a_point& p2, const a_point& p3)
{
        a_point dp1 = p2-p1;//segments between points
        a_point dp2 = p3-p2;
        a_point n = cross(dp1,dp2).normalise();//normal to the plane
        a_point n1 = cross(dp1,n).normalise();//normal to the segments
        a_point n2 = cross(dp2,n).normalise();
        a_point pm1 = (p1+p2)/2.;//mid points of the segments
        a_point pm2 = (p2+p3)/2.;
        a_point dpp = pm2-pm1;//side of triangle
        double d = dpp.norm();
        double a1 = acos(n1*n2);
        double a2 = acos(dpp*n2/d);
        double r;
        if (a1 != 0)
                r = d*sin(a2)/sin(a1);//sine law
        else
                r = 1.e10;//something big
        return pm1-n1*r;
}
//***************************************************************************
//---------------------------------------------------------------------------
void a_point::read(std::istream &in)
{
        in >> x_;
        in >> y_;
        in >> z_;
}
//---------------------------------------------------------------------------
void a_point::write(std::ostream &o) const
{
        o << this->x() << " ";
        o << this->y() << " ";
        o << this->z() << " ";
}
//---------------------------------------------------------------------------
double angle(a_point a, a_point b)
{
        double pi = 3.141592653589793238;
        a.normalise();
        b.normalise();
        double c = a*b;
        double s = cross(a,b).norm();
        if (c<0)
        {
                if (s<0)
                        return acos(-c)+pi;
                else
                        return pi-acos(-c);
        }
        else
        {
                if (s<0)
                        return -acos(c);
                else
                        return acos(c);
        }
}
//---------------------------------------------------------------------------
a_point average_rot(const a_point& a, const a_point& b, double f)
{
        double l1 = a.norm();
        double l2 = b.norm();
        double l = l1+f*(l2-l1);
        a_point axis = cross(a,b).normalise();
        double phi = angle(a,b)*f;
        a_point inter = a;
        a_point pn = l*inter.rotate(axis,phi).normalise();
        return pn;
}