codebase / stimlib

#ifndef RTS_VECTOR_H
#define RTS_VECTOR_H

#include <iostream>
//#include "rts/point.h"

namespace rts
{



template <class T, int N>
struct rtsVector
{
	T v[N];

	CUDA_CALLABLE rtsVector()
	{
		//memset(v, 0, sizeof(T) * N);
		for(int i=0; i<N; i++)
			v[i] = 0;
	}

	//efficiency constructor, makes construction easier for 1D-4D vectors
	CUDA_CALLABLE rtsVector(T x, T y = (T)0.0, T z = (T)0.0, T w = (T)0.0)
	{
		if(N >= 1)
			v[0] = x;
		if(N >= 2)
			v[1] = y;
		if(N >= 3)
			v[2] = z;
		if(N >= 4)
			v[3] = w;
	}

	CUDA_CALLABLE rtsVector(const T(&data)[N])
	{
		memcpy(v, data, sizeof(T) * N);
	}

	CUDA_CALLABLE T len()
	{
        //compute and return the vector length
        T sum_sq = (T)0;
        for(int i=0; i<N; i++)
        {
            sum_sq += v[i] * v[i];
        }
        return std::sqrt(sum_sq);

	}

	CUDA_CALLABLE rtsVector<T, N> cart2sph()
	{
		//convert the vector from cartesian to spherical coordinates
		//x, y, z -> r, theta, phi (where theta = 0 to 2*pi)

		rtsVector<T, N> sph;
		sph[0] = std::sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
		sph[1] = std::atan2(v[1], v[0]);
		sph[2] = std::acos(v[2] / sph[0]);

		return sph;
	}

	CUDA_CALLABLE rtsVector<T, N> sph2cart()
	{
		//convert the vector from cartesian to spherical coordinates
		//r, theta, phi -> x, y, z (where theta = 0 to 2*pi)

		rtsVector<T, N> cart;
		cart[0] = v[0] * std::cos(v[1]) * std::sin(v[2]);
		cart[1] = v[0] * std::sin(v[1]) * std::sin(v[2]);
		cart[2] = v[0] * std::cos(v[2]);

		return cart;
	}

	CUDA_CALLABLE rtsVector<T, N> norm()
	{
        //compute and return the vector norm
        rtsVector<T, N> result;

        //compute the vector length
        T l = len();

        //normalize
        for(int i=0; i<N; i++)
        {
            result.v[i] = v[i] / l;
        }

        return result;
	}

	CUDA_CALLABLE rtsVector<T, 3> cross(rtsVector<T, 3> rhs)
	{
		rtsVector<T, 3> result;

		//compute the cross product (only valid for 3D vectors)
		result[0] = v[1] * rhs[2] - v[2] * rhs[1];
		result[1] = v[2] * rhs[0] - v[0] * rhs[2];
		result[2] = v[0] * rhs[1] - v[1] * rhs[0];

		return result;
	}

    CUDA_CALLABLE T dot(rtsVector<T, N> rhs)
    {
        T result = (T)0;

        for(int i=0; i<N; i++)
            result += v[i] * rhs.v[i];

        return result;

    }

	//arithmetic
	CUDA_CALLABLE rtsVector<T, N> operator+(rtsVector<T, N> rhs)
	{
        rtsVector<T, N> result;

        for(int i=0; i<N; i++)
            result.v[i] = v[i] + rhs.v[i];

        return result;
	}
	CUDA_CALLABLE rtsVector<T, N> operator-(rtsVector<T, N> rhs)
	{
        rtsVector<T, N> result;

        for(int i=0; i<N; i++)
            result.v[i] = v[i] - rhs.v[i];

        return result;
	}
	CUDA_CALLABLE rtsVector<T, N> operator*(T rhs)
	{
        rtsVector<T, N> result;

        for(int i=0; i<N; i++)
            result.v[i] = v[i] * rhs;

        return result;
	}
	CUDA_CALLABLE rtsVector<T, N> operator/(T rhs)
	{
        rtsVector<T, N> result;

        for(int i=0; i<N; i++)
            result.v[i] = v[i] / rhs;

        return result;
	}

	std::string toStr()
	{
		std::stringstream ss;

		ss<<"[";
		for(int i=0; i<N; i++)
		{
			ss<<v[i];
			if(i != N-1)
				ss<<", ";
		}
		ss<<"]";

		return ss.str();
	}

	//bracket operator
	CUDA_CALLABLE T& operator[](int i)
	{
        return v[i];
    }

};


}	//end namespace rts

template <typename T, int N>
std::ostream& operator<<(std::ostream& os, rts::rtsVector<T, N> v)
{
    os<<v.toStr();
    return os;
}

//arithmetic operators
template <typename T, int N>
CUDA_CALLABLE rts::rtsVector<T, N> operator-(rts::rtsVector<T, N> v)
{
    rts::rtsVector<T, N> r;

    //negate the vector
    for(int i=0; i<N; i++)
        r.v[i] = -v.v[i];

    return r;
}

template <typename T, int N>
CUDA_CALLABLE rts::rtsVector<T, N> operator*(T lhs, rts::rtsVector<T, N> rhs)
{
    rts::rtsVector<T, N> r;

    return rhs * lhs;
}

#endif