/*
This file is part of EZRide.
Copyright 2010 Vadim Biktashev & Andrew Foulkes.

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published
by the Free Software Foundation, either version 3 of the License,
or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include "ezride.h"
#include <stdio.h>
#include <math.h>

/* ========================================================================= */

Real dt, dt_h, dt_o_2h, Au, Av, cx, cy, omega, a1;
Real X1, Y1, TH1, grid_h;
int  NX, NY, NX2, NY2, istep;

/* ========================================================================= */

void Perturb (void)
{

  /* "Electrophoresis" perturbation */

int  i,j;
Real Aus, Auc, Avs, Avc;
Real Ux, Uy, Vx, Vy;
Real dt2hAuc, dt2hAus, dt2hAvc, dt2hAvs;

  dt2hAuc = dt_o_2h*Au*cos(TH1);
  dt2hAus = dt_o_2h*Au*sin(TH1);
  dt2hAvc = dt_o_2h*Av*cos(TH1);
  dt2hAvs = dt_o_2h*Av*sin(TH1);

  for (i=1;i<=NX;i++){
   for (j=1;j<=NY;j++){

/* transformed perturbation per analysis is;     *
 *                                               *
 *    eps h = A( cos(th)*dv/dx - sin(th)*dv/dy ) *
 *                                               */



	if(Auc>0)     {
	  if (i==NX) {
	    Ux=-U(NX-2,j)+(4*U(NX-1,j))-(3*U(NX,j));
	     }
          else if (i==NX-1) {
	    Ux=-U(NX-1,j)+(4*U(NX,j))-(3*U(NX-1,j));
	     } 
          else {
	    Ux=-U(i+2,j)+(4*U(i+1,j))-(3*U(i,j));
	     }
          }

	else {
          if(i==1){
            Ux=U(3,j)-(4*U(2,j))+(3*U(1,j));
             }
          else if(i==2){
            Ux=U(2,j)-(4*U(1,j))+(3*U(2,j));
             }
          else{
            Ux=U(i-2,j)-(4*U(i-1,j))+(3*U(i,j));
             }
          }

	if(Avc>0)     {
	  if (i==NX) {
	    Vx=-V(NX-2,j)+(4*V(NX-1,j))-(3*V(NX,j));
	     }
          else if (i==NX-1) {
	    Vx=-V(NX-1,j)+(4*V(NX,j))-(3*V(NX-1,j));
	     } 
          else {
	    Vx=-V(i+2,j)+(4*V(i+1,j))-(3*V(i,j));
	     }
          }

	else {
          if(i==1){
            Vx=V(3,j)-(4*V(2,j))+(3*V(1,j));
             }
          else if(i==2){
            Vx=V(2,j)-(4*V(1,j))+(3*V(2,j));
             }
          else{
            Vx=V(i-2,j)-(4*V(i-1,j))+(3*V(i,j));
             }
          }

        if(Aus>0)     {
	  if (j==NY) {
	    Uy=-U(i,NY-2)+(4*U(i,NY-1))-(3*U(i,NY));
	     }
          else if (j==NY-1) {
	    Uy=-U(i,NY-1)+(4*U(i,NY))-(3*U(i,NY-1));
	     } 
          else {
	    Uy=-U(i,j+2)+(4*U(i,j+1))-(3*U(i,j));
	     }
          }

	else {
          if(j==1){
            Uy=U(i,3)-(4*U(i,2))+(3*U(i,1));
             }
          else if(j==2){
            Uy=U(i,2)-(4*U(i,1))+(3*U(i,2));
             }
          else{
            Uy=U(i,j-2)-(4*U(i,j-1))+(3*U(i,j));
             }
          }
	if(Avs>0)     {
	  if (j==NY) {
	    Vy=-V(i,NY-2)+(4*V(i,NY-1))-(3*V(i,NY));
	     }
          else if (j==NY-1) {
	    Vy=-V(i,NY-1)+(4*V(i,NY))-(3*V(i,NY-1));
	     } 
          else {
            Vy=-V(i,j+2)+(4*V(i,j+1))-(3*V(i,j));
	     }
          }

	else {
          if(j==1){
            Vy=V(i,3)-(4*V(i,2))+(3*V(i,1));
             }
          else if(j==2){
            Vy=V(i,2)+(4*V(i,1))+(3*V(i,2));
             }
          else{
            Vy=V(i,j-2)-(4*V(i,j-1))+(3*V(i,j));
             }
          }

 U1(i,j) = U(i,j) + ( (dt2hAuc*Ux) - (dt2hAus*Uy) );
 V1(i,j) = V(i,j) + ( (dt2hAvc*Vx) - (dt2hAvs*Vy) );
  }
 }
  
 for (i=1;i<=NX;i++){
  for (j=1;j<=NY;j++){
    U(i,j) = U1(i,j);
    V(i,j) = V1(i,j);
  }
 }
}

/* ========================================================================= */