#ifndef INCLUDE_ONCE_VECTOR_FIELD_DOUBLE_GYRE_2D #define INCLUDE_ONCE_VECTOR_FIELD_DOUBLE_GYRE_2D #include #include typedef std::pair 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