#ifndef XYZH
#define XYZH
#include <iostream>
#include <cstdlib>
#include <cmath>
using namespace std;

//4 * atan2(1, 1);
const double pi = 3.141592653589793238462643383279502884197169399375104;
inline double myrand() {return static_cast<double>(rand()) / RAND_MAX;}

inline double weighted_rand() {const double d = myrand(); return d * d;}

struct rtp;

struct xyz {
	double x, y, z;

	xyz(double initial_x = 0, double initial_y = 0, double initial_z = 0)
		: x(initial_x), y(initial_y), z(initial_z) {}

	xyz& operator+=(const xyz& other) {
		x += other.x;
		y += other.y;
		z += other.z;
		return *this;
	}

	xyz& operator-=(const xyz& other) {
		x -= other.x;
		y -= other.y;
		z -= other.z;
		return *this;
	}

	xyz& operator*=(double d) {
		x *= d;
		y *= d;
		z *= d;
		return *this;
	}

	xyz& operator/=(double d) {
		if (d == 0) {
			cerr << "Can't divide " << *this << " by 0.\n";
			exit(EXIT_FAILURE);
		}
		x /= d;
		y /= d;
		z /= d;
		return *this;
	}

	template<double xyz::*A, double xyz::*B>
	xyz& rotate(double theta) {
		double& a = this->*A;
		double& b = this->*B;

		const double r = sqrt(a * a + b * b);
		theta += (a == 0 && b == 0 ? 0 : atan2(b, a));

		a = r * cos(theta);
		b = r * sin(theta);
		return *this;
	}

	//Rotate around x axis, follow right hand rule
	xyz& xrot(double theta) {return rotate<&xyz::y, &xyz::z>(theta);}

	//Rotate around y axis, follow right hand rule
	xyz& yrot(double theta) {return rotate<&xyz::z, &xyz::x>(theta);}

	//Rotate around z axis, follow right hand rule
	xyz& zrot(double theta) {return rotate<&xyz::x, &xyz::y>(theta);}

	//Rotate around line from direction to viewpoint.
	xyz& rot(const xyz& direction, const xyz& viewpoint, double theta);

	friend ostream& operator<<(ostream& ost, const xyz& point) {
		return ost << "("
			<< point.x << ", "
			<< point.y << ", "
			<< point.z << ")";
	}

	operator rtp() const;
};

inline const xyz operator+(xyz xyz1, const xyz& xyz2) {return xyz1 += xyz2;}
inline const xyz operator-(xyz xyz1, const xyz& xyz2) {return xyz1 -= xyz2;}

inline const xyz operator*(xyz xyz, double d) {return xyz *= d;}
inline const xyz operator*(double d, xyz xyz) {return xyz *= d;}
inline const xyz operator/(xyz xyz, double d) {return xyz /= d;}

//Rotate around the line connecting direction and viewpoint.

inline xyz& xyz::rot(const xyz& direction, const xyz& viewpoint, double theta)
{
	const xyz d = viewpoint - direction;

	const double xtheta = d.y == 0 && d.z == 0 ? 0 : atan2(d.y, d.z);
	const double ytheta = d.x == 0 && d.z == 0 ? 0 : atan2(d.x, d.z);

	xrot( xtheta);
	yrot( ytheta);
	zrot(  theta);
	yrot(-ytheta);
	xrot(-xtheta);
	
	return *this;
}


struct rtp {
	const double r;		//radius
	const double theta;	//longitude
	const double phi;	//colatitude: angle down from north pole

	rtp(double initial_r = 0, double initial_theta = 0,
		double initial_phi = pi / 2)
		: r(initial_r), theta(initial_theta), phi(initial_phi) {}

	operator xyz() const {
		const double s = sin(phi);
		return r * xyz(s * cos(theta), s * sin(theta), cos(phi));
	}
};

inline xyz::operator rtp() const
{
	return rtp(
		sqrt(x * x + y * y + z * z),
		x == 0 && y == 0 ? 0 : atan2(y, x),
		z == 0 ? pi / 2 : atan2(sqrt(x * x + y * y), z)
	);
}
#endif