#ifndef INCLUDE_ONCE_VECTOR_FIELD_DOUBLE_GYRE_2D
#define INCLUDE_ONCE_VECTOR_FIELD_DOUBLE_GYRE_2D

#include <cmath>
#include <algorithm>

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

#ifndef M_PI
#define M_PI 3.14159265358979323846264338327950288
#endif

class doublegyre2d
{
public:
	doublegyre2d(double a = 0.1, double eps = 0.25, double omega = M_PI / 5) : A(a), EPS(eps), OMEGA(omega) {}

	// Samples the vector field of the Double Gyre flow
	Vec2d sample(double x, double y, double t)
	{
		return Vec2d(
			-A * M_PI*sin(M_PI*(EPS*sin(OMEGA*t)*x*x + (1 - 2 * EPS*sin(OMEGA*t))*x))*cos(M_PI*y),
			A * M_PI*(2 * EPS*sin(OMEGA*t)*x - 2 * EPS*sin(OMEGA*t) + 1)* cos(M_PI*(EPS*sin(OMEGA*t)*x*x + (1 - 2 * EPS*sin(OMEGA*t))*x))*sin(M_PI*y));
	}

	// Samples the x-partial of the vector field of the Double Gyre flow
	Vec2d sample_dx(double x, double y, double t)
	{
		return Vec2d(
			-A * M_PI*M_PI * (2 * EPS*sin(OMEGA*t)*x - 2 * EPS*sin(OMEGA*t) + 1)*cos(M_PI*(EPS*sin(OMEGA*t)*x*x + (1 - 2 * EPS*sin(OMEGA*t))*x))*cos(M_PI*y),
			2 * A*EPS*M_PI* sin(OMEGA*t)*cos(M_PI*(EPS*sin(OMEGA*t)*x*x + (1 - 2 * EPS*sin(OMEGA*t))*x))*sin(M_PI*y) - A * M_PI*M_PI * (2 * EPS*sin(OMEGA*t)*x - 2 * EPS*sin(OMEGA*t) + 1) * (2 * EPS*sin(OMEGA*t)*x - 2 * EPS*sin(OMEGA*t) + 1) * sin(M_PI*(EPS*sin(OMEGA*t)*x*x + (1 - 2 * EPS*sin(OMEGA*t))*x))*sin(M_PI*y));
	}

	// Samples the y-partial of the vector field of the Double Gyre flow
	Vec2d sample_dy(double x, double y, double t)
	{
		return Vec2d(
			A * M_PI*M_PI * sin(M_PI*(EPS*sin(OMEGA*t)*x*x + (1 - 2 * EPS*sin(OMEGA*t))*x))*sin(M_PI*y),
			A * M_PI*M_PI * (2 * EPS*sin(OMEGA*t)*x - 2 * EPS*sin(OMEGA*t) + 1) * cos(M_PI*(EPS*sin(OMEGA*t)*x*x + (1 - 2 * EPS*sin(OMEGA*t))*x))*cos(M_PI*y));
	}

	// Samples the t-partial of the vector field of the Double Gyre flow
	Vec2d sample_dt(double x, double y, double t)
	{
		return Vec2d(
			-A * M_PI*M_PI * (EPS*OMEGA*cos(OMEGA*t)*x*x - 2 * EPS*OMEGA*cos(OMEGA*t)*x)*cos(M_PI*(EPS*sin(OMEGA*t)*x*x + (1 - 2 * EPS*sin(OMEGA*t))*x))*cos(M_PI*y),
			A * M_PI *	(2 * EPS*OMEGA*cos(OMEGA*t)*x - 2 * EPS*OMEGA*cos(OMEGA*t))*cos(M_PI*(EPS*sin(OMEGA*t)*x*x + (1 - 2 * EPS*sin(OMEGA*t))*x))*sin(M_PI*y) - A * M_PI*M_PI * (2 * EPS*sin(OMEGA*t)*x - 2 * EPS*sin(OMEGA*t) + 1)*(EPS*OMEGA*cos(OMEGA*t)*x*x - 2 * EPS*OMEGA*cos(OMEGA*t)*x) * sin(M_PI*(EPS*sin(OMEGA*t)*x*x + (1 - 2 * EPS*sin(OMEGA*t))*x))*sin(M_PI*y));
	}

private:
	double A; 
	double EPS;
	double OMEGA;
};

#endif