#ifndef INCLUDE_ONCE_VECTOR_FIELD_FOUR_CENTERS_2D
#define INCLUDE_ONCE_VECTOR_FIELD_FOUR_CENTERS_2D

#include <cmath>
#include <algorithm>

typedef std::pair<double, double> Vec2d;

class fourcenters2d
{
public:
	fourcenters2d(double Scale = 1, double Al_t = 1) : scale(Scale), al_t(Al_t) {}

	// Samples the vector field of the four centers flow
	Vec2d sample(double x, double y, double t)
	{
		return Vec2d(
			exp(-y * y - x * x)*(al_t*y*exp(y*y + x * x) - 6 * scale*cos(al_t*t)*sin(al_t*t)*y*y*y + (12 * scale*(cos(al_t*t)*cos(al_t*t)) - 6 * scale)*x*y*y + (6 * scale*cos(al_t*t)*sin(al_t*t)*x*x + 6 * scale*cos(al_t*t)*sin(al_t*t))*y + (3 * scale - 6 * scale*(cos(al_t*t)*cos(al_t*t)))*x),
			-exp(-y * y - x * x)*(al_t*x*exp(y*y + x * x) - 6 * scale*cos(al_t*t)*sin(al_t*t)*x*y*y + ((12 * scale*(cos(al_t*t)*cos(al_t*t)) - 6 * scale)*x*x - 6 * scale*(cos(al_t*t)*cos(al_t*t)) + 3 * scale)*y + 6 * scale*cos(al_t*t)*sin(al_t*t)*x*x*x - 6 * scale*cos(al_t*t)*sin(al_t*t)*x));
	}

	// Samples the x-partial of the vector field of the four centers flow
	Vec2d sample_dx(double x, double y, double t)
	{
		return Vec2d(
			3 * scale*(4 * cos(al_t*t)*sin(al_t*t)*x*y*y*y - 8 * (cos(al_t*t)*cos(al_t*t)) * x*x * y*y + 4 * x*x * y*y + 4 * (cos(al_t*t)*cos(al_t*t)) * y*y - 2 * y*y - 4 * cos(al_t*t)*sin(al_t*t)*x*x*x * y + 4 * (cos(al_t*t)*cos(al_t*t)) * x*x - 2 * x*x - 2 * (cos(al_t*t)*cos(al_t*t)) + 1)*exp(-y * y - x * x),
			-exp(-y * y - x * x)*(al_t*exp(y*y + x * x) + 12 * scale*cos(al_t*t)*sin(al_t*t)*x*x * y*y - 6 * scale*cos(al_t*t)*sin(al_t*t)*y*y - 24 * scale*(cos(al_t*t)*cos(al_t*t)) * x*x*x * y + 12 * scale*x*x*x * y + 36 * scale*(cos(al_t*t)*cos(al_t*t)) * x*y - 18 * scale*x*y - 12 * scale*cos(al_t*t)*sin(al_t*t)*x*x*x*x + 30 * scale*cos(al_t*t)*sin(al_t*t)*x*x - 6 * scale*cos(al_t*t)*sin(al_t*t)));
	}

	// Samples the y-partial of the vector field of the four centers flow
	Vec2d sample_dy(double x, double y, double t)
	{
		return Vec2d(
			exp(-y * y - x * x)*(al_t*exp(y*y + x * x) + 12 * scale*cos(al_t*t)*sin(al_t*t)*y*y*y*y - 24 * scale*(cos(al_t*t)*cos(al_t*t)) * x*y*y*y + 12 * scale*x*y*y*y - 12 * scale*cos(al_t*t)*sin(al_t*t)*x*x * y*y - 30 * scale*cos(al_t*t)*sin(al_t*t)*y*y + 36 * scale*(cos(al_t*t)*cos(al_t*t)) * x*y - 18 * scale*x*y + 6 * scale*cos(al_t*t)*sin(al_t*t)*x*x + 6 * scale*cos(al_t*t)*sin(al_t*t)),
			-3 * scale*(4 * cos(al_t*t)*sin(al_t*t)*x*y*y*y - 8 * (cos(al_t*t)*cos(al_t*t)) * x*x * y*y + 4 * x*x * y*y + 4 * (cos(al_t*t)*cos(al_t*t)) * y*y - 2 * y*y - 4 * cos(al_t*t)*sin(al_t*t)*x*x*x * y + 4 * (cos(al_t*t)*cos(al_t*t)) * x*x - 2 * x*x - 2 * (cos(al_t*t)*cos(al_t*t)) + 1)*exp(-y * y - x * x));
	}

	// Samples the t-partial of the vector field of the four centers flow
	Vec2d sample_dt(double x, double y, double t)
	{
		return Vec2d(
			6 * al_t*scale*((sin(al_t*t)*sin(al_t*t)) * y*y*y - (cos(al_t*t)*cos(al_t*t)) * y*y*y - 4 * cos(al_t*t)*sin(al_t*t)*x*y*y - (sin(al_t*t)*sin(al_t*t)) * x*x * y + (cos(al_t*t)*cos(al_t*t)) * x*x * y - (sin(al_t*t)*sin(al_t*t)) * y + (cos(al_t*t)*cos(al_t*t)) * y + 2 * cos(al_t*t)*sin(al_t*t)*x)*exp(-y * y - x * x),
			-6 * al_t*scale*((sin(al_t*t)*sin(al_t*t)) * x*y*y - (cos(al_t*t)*cos(al_t*t)) * x*y*y - 4 * cos(al_t*t)*sin(al_t*t)*x*x * y + 2 * cos(al_t*t)*sin(al_t*t)*y - (sin(al_t*t)*sin(al_t*t)) * x*x*x + (cos(al_t*t)*cos(al_t*t)) * x*x*x + (sin(al_t*t)*sin(al_t*t)) * x - (cos(al_t*t)*cos(al_t*t)) * x)*exp(-y * y - x * x));
	}

private:
	double scale;
	double al_t;
};

#endif