Vector.cc

Go to the documentation of this file.

//Copyright (c) 2013-2020, The MercuryDPM Developers Team. All rights reserved.
//For the list of developers, see <http://www.MercuryDPM.org/Team>.
//
//Redistribution and use in source and binary forms, with or without
//modification, are permitted provided that the following conditions are met:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name MercuryDPM nor the
//    names of its contributors may be used to endorse or promote products
//    derived from this software without specific prior written permission.
//
//THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
//ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
//WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
//DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
//DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
//(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
//LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
//ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
//(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
//SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

#include "Vector.h"
#include "SmallVector.h"

Vec3D::Vec3D(const SmallVector<3>& vector)
{
    X = vector[0];
    Y = vector[1];
    Z = vector[2];
}

void Vec3D::setZero()
{
    X = 0.0;
    Y = 0.0;
    Z = 0.0;
}

void Vec3D::setNaN()
{
    X = constants::NaN;
    Y = constants::NaN;
    Z = constants::NaN;
}

bool Vec3D::isNaN() const
{
    return std::isnan(X) || std::isnan(Y) || std::isnan(Z);
}

Mdouble Vec3D::dot(const Vec3D& a, const Vec3D& b)
{
    return a.X * b.X + a.Y * b.Y + a.Z * b.Z;
}

Vec3D Vec3D::max(const Vec3D& a, const Vec3D& b)
{
    return Vec3D(std::max(a.X, b.X), std::max(a.Y, b.Y), std::max(a.Z, b.Z));
}

Vec3D Vec3D::min(const Vec3D& a, const Vec3D& b)
{
    return Vec3D(std::min(a.X, b.X), std::min(a.Y, b.Y), std::min(a.Z, b.Z));
}

Vec3D Vec3D::square(const Vec3D& a)
{
    return Vec3D(a.X * a.X, a.Y * a.Y, a.Z * a.Z);
}

void Vec3D::normalise()
{
    Mdouble length2 = this->getLengthSquared();
    if (length2 == 0)
    {
        logger(ERROR, "Normalizing a vector of length 0");
    }
    *this /= std::sqrt(length2);
}

void Vec3D::setLength(Mdouble length)
{
    this->normalise();
    *this *= length;
}

Vec3D Vec3D::sqrt(const Vec3D& a)
{
    return Vec3D(std::sqrt(a.X), std::sqrt(a.Y), std::sqrt(a.Z));
}

Vec3D Vec3D::cross(const Vec3D& a, const Vec3D& b)
{
    return Vec3D(a.Y * b.Z - a.Z * b.Y, a.Z * b.X - a.X * b.Z, a.X * b.Y - a.Y * b.X);
}

Mdouble Vec3D::getDistance(const Vec3D& a, const Vec3D& b)
{
    return std::sqrt(getDistanceSquared(a, b));
}

Mdouble Vec3D::getLengthSquared() const
{
    return (X * X + Y * Y + Z * Z);
}

Mdouble Vec3D::getComponent(const int index) const
{
    switch (index)
    {
        case 0:
            return X;
        case 1:
            return Y;
        case 2:
            return Z;
        default:
            logger(ERROR, "[Vector::getComponent] Index = %, which is too high for a 3D vector (should be 0-2).",
                   index);
            return 0;
    }
}

void Vec3D::setComponent(const int index, const double val)
{
    switch (index)
    {
        case 0:
            X = val;
            break;
        case 1:
            Y = val;
            break;
        case 2:
            Z = val;
            break;
        default:
            logger(ERROR, "[Vector::setComponent] Index = %, which is too high for a 3D vector (should be 0-2).",
                   index);
    }
}


Mdouble Vec3D::getRadialCoordinateSquared() const
{
    return X * X + Y * Y;
}

Mdouble Vec3D::getRadialCoordinate() const
{
    return std::sqrt(X * X + Y * Y);
}

Vec3D Vec3D::getCylindricalCoordinates() const
{
    return Vec3D(std::sqrt(X * X + Y * Y), std::atan2(Y, X), Z);
}

Vec3D Vec3D::getCylindricalTensorField(const Vec3D& p) const
{
    //define sin(A)=y/r, cos(A)=x/r
    Mdouble r = std::sqrt(p.X * p.X + p.Y * p.Y);
    Mdouble s = p.Y / r;
    Mdouble c = p.X / r;
    if (r == 0)
    {
        s = 0;
        c = 1;
    }
    return Vec3D(X * c + Y * s, -X * s + Y * c, Z);
}

Vec3D Vec3D::getFromCylindricalCoordinates() const
{
    return Vec3D(X * std::cos(Y), X * std::sin(Y), Z);
}

bool Vec3D::isEqualTo(const Vec3D& other, const double tol) const
{
    if ((Vec3D::getLengthSquared(*this - other)) <= tol * tol)
    {
        return true;
    }
    else
    {
        return false;
    }
}

//  void ConvertToCylindricalCoordinates()
//  {
//      double R = sqrt(X*X+Y*Y); Y = atan2(Y,X); X = R; return;
//  }
// 
//  void ConvertFromCylindricalCoordinates()
//  {
//      double Xnew = X * cos(Y); Y = X * sin(Y); X = Xnew; return;
//  }

Mdouble Vec3D::getLength() const
{
    return std::sqrt(getLengthSquared());
}

Mdouble Vec3D::getLength(const Vec3D& a)
{
    return a.getLength();
}

Vec3D Vec3D::getUnitVector(const Vec3D& a)
{
    Mdouble Length2 = a.getLengthSquared();
    if (Length2 != 0.0)
        return a / std::sqrt(Length2);
    else
        return Vec3D(0, 0, 0);
}

std::ostream& operator<<(std::ostream& os, const Vec3D& a)
{
    os << a.X << ' ' << a.Y << ' ' << a.Z;
    return os;
}

std::istream& operator>>(std::istream& is, Vec3D& a)
{
    //TW: clearing the stream avoids the nasty problem that the failbit is set to true if numbers below DBL_MIN=1e-308 are read.
    is >> a.X; is.clear();
    is >> a.Y; is.clear();
    is >> a.Z; //is.clear();
    return is;
}