// Vel2dEx.cpp
//
// J.Liu, R.Cali
// ColoState, Jan.-Jun. 2007


#include <cmath>
#include "PtVec2d.h"

#ifndef M_PI
#define M_PI 3.141592653589793
#endif  // M_PI


// A band of 4 eddies on the unit square

PtVec2d vel2dExEddyBand(PtVec2d P)
{
   double x = P.xCrd()-0.5;
   double y = P.yCrd()-0.5;
   double w = 1-4*x*x;
   double v1 = pow(w,4)*cos(4*M_PI*y);
   double v2 = 8*x*pow(w,3)*sin(4*M_PI*y)/M_PI;
   return PtVec2d(v1, v2);
}


// Four eddies on the unit square

PtVec2d vel2dExEddyFour(PtVec2d P)
{
   double x = P.xCrd()-0.5;
   double y = P.yCrd()-0.5;
   double v1 = -sin(2*M_PI*x)*cos(2*M_PI*y);
   double v2 = cos(2*M_PI*x)*sin(2*M_PI*y);
   return PtVec2d(v1, v2);
}


// A single eddy on the unit square

PtVec2d vel2dExEddyOne(PtVec2d P) 
{
   double x = P.xCrd() - 0.5;
   double y = P.yCrd() - 0.5;
   double w = 1-4*y*y;
   double v1 = cos(M_PI*x)/M_PI*32*y*pow(w,3);
   double v2 = -sin(M_PI*x)*pow(w,4);
   return PtVec2d(v1, v2);
}


// A rotating velocity field on the unit square 

PtVec2d vel2dExRotating(PtVec2d P) 
{
   double x = P.xCrd()-0.5;
   double y = P.yCrd()-0.5;
   return PtVec2d(y, -x);
}


// A swirling velocity field on the unit square 
// An example proposed by R.LeVeque

PtVec2d vel2dExSwirling(PtVec2d P) 
{
   double x = P.xCrd();
   double y = P.yCrd();
   double wx = sin(M_PI*x);
   double wy = sin(M_PI*y);
   double v1 = wx*wx*sin(2*M_PI*y);
   double v2 = -sin(2*M_PI*x)*wy*wy;
   return PtVec2d(v1, v2);
}