/*------------------------------------------------------------------------
 * Copyright (C) 2011 For the list of authors, see file AUTHORS.
 *
 * This file is part of SOFI3D.
 * 
 * SOFI3D 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, version 2.0 of the License only.
 * 
 * SOFI3D 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 SOFI3D. See file COPYING and/or 
  * <http://www.gnu.org/licenses/gpl-2.0.html>.
--------------------------------------------------------------------------*/
/*------------------------------------------------------------------------
 *   updating particle velocities at gridpoints [nx1...nx2][ny1...ny2][nz1...nz2]
 *   by a staggered grid finite difference scheme of FDORDER order accuracy in space
 *   and second order accuracy in time
 *
 *  ----------------------------------------------------------------------*/

#include "fd.h"

double update_v_acoustic_PML(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2,
		int nt, float *** vx, float *** vy, float *** vz, float *** sxx,
		float *** vx1, float *** vy2, float *** vz3, float  ***  rho, float **  srcpos_loc, float ** signals,
		int nsrc, float *** absorb_coeffx, float *** absorb_coeffy, float *** absorb_coeffz, int * stype){


	extern float DT, DX, DY, DZ;
	extern int MYID, FDORDER, LOG, FDCOEFF;
	extern FILE *FP;
	extern int OUTNTIMESTEPINFO;

	register int i, j, k, l;
	float  amp, rjp, rkp, rip;
	double time=0.0, time1=0.0, time2=0.0;
	register float b1, b2, b3, b4, b5, b6, dx, dy, dz; 
	register float sxx_x, syy_y, szz_z, PML1, PML2;

	dx=1.0/DX;
	dy=1.0/DY;
	dz=1.0/DZ;	

	if (LOG)
		if ((MYID==0) && ((nt+(OUTNTIMESTEPINFO-1))%OUTNTIMESTEPINFO)==0) time1=MPI_Wtime();

	switch (FDORDER){
	case 2 :


		for (j=ny1;j<=ny2;j++){
			for (i=nx1;i<=nx2;i++){
				for (k=nz1;k<=nz2;k++){

					rjp=0.5*(rho[j][i][k]+rho[j+1][i][k+1]);
					rkp=0.5*(rho[j][i][k]+rho[j][i][k+1]);
					rip=0.5*(rho[j][i][k]+rho[j][i+1][k]);

					PML1=((((1.0/DT)-(absorb_coeffx[j][i][k]/2.0)))/(((1.0/DT)+(absorb_coeffx[j][i][k]/2.0))));
					PML2=(1.0/(((1.0/DT)+(absorb_coeffx[j][i][k]/2.0))));

					vx1[j][i][k] = (PML1*vx1[j][i][k]) + PML2*dx*(sxx[j][i+1][k]-sxx[j][i][k])/rip;

					PML1=((((1.0/DT)-(absorb_coeffy[j][i][k]/2.0)))/(((1.0/DT)+(absorb_coeffy[j][i][k]/2.0))));
					PML2=(1.0/(((1.0/DT)+(absorb_coeffy[j][i][k]/2.0))));

					vy2[j][i][k] = (PML1*vy2[j][i][k]) + PML2*dy*(sxx[j+1][i][k]-sxx[j][i][k])/rjp;

					PML1=((((1.0/DT)-(absorb_coeffz[j][i][k]/2.0)))/(((1.0/DT)+(absorb_coeffz[j][i][k]/2.0))));
					PML2=(1.0/(((1.0/DT)+(absorb_coeffz[j][i][k]/2.0))));

					vz3[j][i][k] = (PML1*vz3[j][i][k]) + PML2*dz*(sxx[j][i][k+1]-sxx[j][i][k])/rkp;

					vx[j][i][k] = vx1[j][i][k];
					vy[j][i][k] = vy2[j][i][k];
					vz[j][i][k] = vz3[j][i][k];

				}
			}
		}
		break;

	case 4 :
		b1=9.0/8.0; b2=-1.0/24.0; /* Taylor coefficients */
		if(FDCOEFF==2){
			b1=1.1382; b2=-0.046414;} /* Holberg coefficients E=0.1 %*/

			for (j=ny1;j<=ny2;j++){
				for (i=nx1;i<=nx2;i++){
					for (k=nz1;k<=nz2;k++){

						rjp=0.5*(rho[j][i][k]+rho[j+1][i][k+1]);
						rkp=0.5*(rho[j][i][k]+rho[j][i][k+1]);
						rip=0.5*(rho[j][i][k]+rho[j][i+1][k]);

						sxx_x = b1*(sxx[j][i+1][k]-sxx[j][i][k])+
								b2*(sxx[j][i+2][k]-sxx[j][i-1][k]);
						syy_y =	b1*(sxx[j+1][i][k]-sxx[j][i][k])+
								b2*(sxx[j+2][i][k]-sxx[j-1][i][k]);
						szz_z = b1*(sxx[j][i][k+1]-sxx[j][i][k])+
								b2*(sxx[j][i][k+2]-sxx[j][i][k-1]);

						PML1=((((1.0/DT)-(absorb_coeffx[j][i][k]/2.0)))/(((1.0/DT)+(absorb_coeffx[j][i][k]/2.0))));
						PML2=(1.0/(((1.0/DT)+(absorb_coeffx[j][i][k]/2.0))));

						vx1[j][i][k] = (PML1*vx1[j][i][k]) + PML2*dx*sxx_x/rip;

						PML1=((((1.0/DT)-(absorb_coeffy[j][i][k]/2.0)))/(((1.0/DT)+(absorb_coeffy[j][i][k]/2.0))));
						PML2=(1.0/(((1.0/DT)+(absorb_coeffy[j][i][k]/2.0))));

						vy2[j][i][k] = (PML1*vy2[j][i][k]) + PML2*dy*syy_y/rjp;

						PML1=((((1.0/DT)-(absorb_coeffz[j][i][k]/2.0)))/(((1.0/DT)+(absorb_coeffz[j][i][k]/2.0))));
						PML2=(1.0/(((1.0/DT)+(absorb_coeffz[j][i][k]/2.0))));

						vz3[j][i][k] = (PML1*vz3[j][i][k]) + PML2*dz*szz_z/rkp;

						vx[j][i][k] = vx1[j][i][k];
						vy[j][i][k] = vy2[j][i][k];
						vz[j][i][k] = vz3[j][i][k];


					}
				}
			}
			break;

	case 6 :
		b1=75.0/64.0; b2=-25.0/384.0; b3=3.0/640.0; /* Taylor coefficients */
		if(FDCOEFF==2){
			b1=1.1965; b2=-0.078804; b3=0.0081781;}   /* Holberg coefficients E=0.1 %*/

		for (j=ny1;j<=ny2;j++){
			for (i=nx1;i<=nx2;i++){
				for (k=nz1;k<=nz2;k++){

					rjp=0.5*(rho[j][i][k]+rho[j+1][i][k+1]);
					rkp=0.5*(rho[j][i][k]+rho[j][i][k+1]);
					rip=0.5*(rho[j][i][k]+rho[j][i+1][k]);

					sxx_x = b1*(sxx[j][i+1][k]-sxx[j][i][k])+
							b2*(sxx[j][i+2][k]-sxx[j][i-1][k])+
							b3*(sxx[j][i+3][k]-sxx[j][i-2][k]);
					syy_y =	b1*(sxx[j+1][i][k]-sxx[j][i][k])+
							b2*(sxx[j+2][i][k]-sxx[j-1][i][k])+
							b3*(sxx[j+3][i][k]-sxx[j-2][i][k]);
					szz_z = b1*(sxx[j][i][k+1]-sxx[j][i][k])+
							b2*(sxx[j][i][k+2]-sxx[j][i][k-1])+
							b3*(sxx[j][i][k+3]-sxx[j][i][k-2]);

					PML1=((((1.0/DT)-(absorb_coeffx[j][i][k]/2.0)))/(((1.0/DT)+(absorb_coeffx[j][i][k]/2.0))));
					PML2=(1.0/(((1.0/DT)+(absorb_coeffx[j][i][k]/2.0))));

					vx1[j][i][k] = (PML1*vx1[j][i][k]) + PML2*dx*sxx_x/rip;

					PML1=((((1.0/DT)-(absorb_coeffy[j][i][k]/2.0)))/(((1.0/DT)+(absorb_coeffy[j][i][k]/2.0))));
					PML2=(1.0/(((1.0/DT)+(absorb_coeffy[j][i][k]/2.0))));

					vy2[j][i][k] = (PML1*vy2[j][i][k]) + PML2*dy*syy_y/rjp;

					PML1=((((1.0/DT)-(absorb_coeffz[j][i][k]/2.0)))/(((1.0/DT)+(absorb_coeffz[j][i][k]/2.0))));
					PML2=(1.0/(((1.0/DT)+(absorb_coeffz[j][i][k]/2.0))));

					vz3[j][i][k] = (PML1*vz3[j][i][k]) + PML2*dz*szz_z/rkp;

					vx[j][i][k] = vx1[j][i][k];
					vy[j][i][k] = vy2[j][i][k];
					vz[j][i][k] = vz3[j][i][k];

				}
			}
		}
		break;

	case 8 :
		b1=1225.0/1024.0; b2=-245.0/3072.0; b3=49.0/5120.0; b4=-5.0/7168.0; /* Taylor coefficients */
		if(FDCOEFF==2){
			b1=1.2257; b2=-0.099537; b3=0.018063; b4=-0.0026274;} /* Holberg coefficients E=0.1 %*/

		for (j=ny1;j<=ny2;j++){
			for (i=nx1;i<=nx2;i++){
				for (k=nz1;k<=nz2;k++){

					rjp=0.5*(rho[j][i][k]+rho[j+1][i][k+1]);
					rkp=0.5*(rho[j][i][k]+rho[j][i][k+1]);
					rip=0.5*(rho[j][i][k]+rho[j][i+1][k]);

					sxx_x = b1*(sxx[j][i+1][k]-sxx[j][i][k])+
							b2*(sxx[j][i+2][k]-sxx[j][i-1][k])+
							b3*(sxx[j][i+3][k]-sxx[j][i-2][k])+
							b4*(sxx[j][i+4][k]-sxx[j][i-3][k]);
					syy_y =	b1*(sxx[j+1][i][k]-sxx[j][i][k])+
							b2*(sxx[j+2][i][k]-sxx[j-1][i][k])+
							b3*(sxx[j+3][i][k]-sxx[j-2][i][k])+
							b4*(sxx[j+4][i][k]-sxx[j-3][i][k]);
					szz_z = b1*(sxx[j][i][k+1]-sxx[j][i][k])+
							b2*(sxx[j][i][k+2]-sxx[j][i][k-1])+
							b3*(sxx[j][i][k+3]-sxx[j][i][k-2])+
							b4*(sxx[j][i][k+4]-sxx[j][i][k-3]);

					PML1=((((1.0/DT)-(absorb_coeffx[j][i][k]/2.0)))/(((1.0/DT)+(absorb_coeffx[j][i][k]/2.0))));
					PML2=(1.0/(((1.0/DT)+(absorb_coeffx[j][i][k]/2.0))));

					vx1[j][i][k] = (PML1*vx1[j][i][k]) + PML2*dx*sxx_x/rip;

					PML1=((((1.0/DT)-(absorb_coeffy[j][i][k]/2.0)))/(((1.0/DT)+(absorb_coeffy[j][i][k]/2.0))));
					PML2=(1.0/(((1.0/DT)+(absorb_coeffy[j][i][k]/2.0))));

					vy2[j][i][k] = (PML1*vy2[j][i][k]) + PML2*dy*syy_y/rjp;

					PML1=((((1.0/DT)-(absorb_coeffz[j][i][k]/2.0)))/(((1.0/DT)+(absorb_coeffz[j][i][k]/2.0))));
					PML2=(1.0/(((1.0/DT)+(absorb_coeffz[j][i][k]/2.0))));

					vz3[j][i][k] = (PML1*vz3[j][i][k]) + PML2*dz*szz_z/rkp;

					vx[j][i][k] = vx1[j][i][k];
					vy[j][i][k] = vy2[j][i][k];
					vz[j][i][k] = vz3[j][i][k];

				}
			}
		}
		break;

	case 10 :
		b1=19845.0/16384.0; b2=-735.0/8192.0; b3=567.0/40960.0; b4=-405.0/229376.0; b5=35.0/294912.0; /* Taylor coefficients */
		if(FDCOEFF==2){
			b1=1.2415; b2=-0.11231; b3=0.026191; b4=-0.0064682; b5=0.001191;} /* Holberg coefficients E=0.1 %*/

		for (j=ny1;j<=ny2;j++){
			for (i=nx1;i<=nx2;i++){
				for (k=nz1;k<=nz2;k++){

					rjp=0.5*(rho[j][i][k]+rho[j+1][i][k+1]);
					rkp=0.5*(rho[j][i][k]+rho[j][i][k+1]);
					rip=0.5*(rho[j][i][k]+rho[j][i+1][k]);

					sxx_x = b1*(sxx[j][i+1][k]-sxx[j][i][k])+
							b2*(sxx[j][i+2][k]-sxx[j][i-1][k])+
							b3*(sxx[j][i+3][k]-sxx[j][i-2][k])+
							b4*(sxx[j][i+4][k]-sxx[j][i-3][k])+
							b5*(sxx[j][i+5][k]-sxx[j][i-4][k]);
					syy_y =	b1*(sxx[j+1][i][k]-sxx[j][i][k])+
							b2*(sxx[j+2][i][k]-sxx[j-1][i][k])+
							b3*(sxx[j+3][i][k]-sxx[j-2][i][k])+
							b4*(sxx[j+4][i][k]-sxx[j-3][i][k])+
							b5*(sxx[j+5][i][k]-sxx[j-4][i][k]);
					szz_z = b1*(sxx[j][i][k+1]-sxx[j][i][k])+
							b2*(sxx[j][i][k+2]-sxx[j][i][k-1])+
							b3*(sxx[j][i][k+3]-sxx[j][i][k-2])+
							b4*(sxx[j][i][k+4]-sxx[j][i][k-3])+
							b5*(sxx[j][i][k+5]-sxx[j][i][k-4]);

					PML1=((((1.0/DT)-(absorb_coeffx[j][i][k]/2.0)))/(((1.0/DT)+(absorb_coeffx[j][i][k]/2.0))));
					PML2=(1.0/(((1.0/DT)+(absorb_coeffx[j][i][k]/2.0))));

					vx1[j][i][k] = (PML1*vx1[j][i][k]) + PML2*dx*sxx_x/rip;

					PML1=((((1.0/DT)-(absorb_coeffy[j][i][k]/2.0)))/(((1.0/DT)+(absorb_coeffy[j][i][k]/2.0))));
					PML2=(1.0/(((1.0/DT)+(absorb_coeffy[j][i][k]/2.0))));

					vy2[j][i][k] = (PML1*vy2[j][i][k]) + PML2*dy*syy_y/rjp;

					PML1=((((1.0/DT)-(absorb_coeffz[j][i][k]/2.0)))/(((1.0/DT)+(absorb_coeffz[j][i][k]/2.0))));
					PML2=(1.0/(((1.0/DT)+(absorb_coeffz[j][i][k]/2.0))));

					vz3[j][i][k] = (PML1*vz3[j][i][k]) + PML2*dz*szz_z/rkp;

					vx[j][i][k] = vx1[j][i][k];
					vy[j][i][k] = vy2[j][i][k];
					vz[j][i][k] = vz3[j][i][k];

				}
			}
		}
		break;

	case 12 :

		/* Taylor coefficients */
		b1=160083.0/131072.0; b2=-12705.0/131072.0; b3=22869.0/1310720.0; 
		b4=-5445.0/1835008.0; b5=847.0/2359296.0; b6=-63.0/2883584;

		/* Holberg coefficients E=0.1 %*/
		if(FDCOEFF==2){
			b1=1.2508; b2=-0.12034; b3=0.032131; b4=-0.010142; b5=0.0029857; b6=-0.00066667;}

		for (j=ny1;j<=ny2;j++){
			for (i=nx1;i<=nx2;i++){
				for (k=nz1;k<=nz2;k++){

					rjp=0.5*(rho[j][i][k]+rho[j+1][i][k+1]);
					rkp=0.5*(rho[j][i][k]+rho[j][i][k+1]);
					rip=0.5*(rho[j][i][k]+rho[j][i+1][k]);

					sxx_x = b1*(sxx[j][i+1][k]-sxx[j][i][k])+
							b2*(sxx[j][i+2][k]-sxx[j][i-1][k])+
							b3*(sxx[j][i+3][k]-sxx[j][i-2][k])+
							b4*(sxx[j][i+4][k]-sxx[j][i-3][k])+
							b5*(sxx[j][i+5][k]-sxx[j][i-4][k])+
							b6*(sxx[j][i+6][k]-sxx[j][i-5][k]);
					syy_y =	b1*(sxx[j+1][i][k]-sxx[j][i][k])+
							b2*(sxx[j+2][i][k]-sxx[j-1][i][k])+
							b3*(sxx[j+3][i][k]-sxx[j-2][i][k])+
							b4*(sxx[j+4][i][k]-sxx[j-3][i][k])+
							b5*(sxx[j+5][i][k]-sxx[j-4][i][k])+
							b6*(sxx[j+6][i][k]-sxx[j-5][i][k]);
					szz_z = b1*(sxx[j][i][k+1]-sxx[j][i][k])+
							b2*(sxx[j][i][k+2]-sxx[j][i][k-1])+
							b3*(sxx[j][i][k+3]-sxx[j][i][k-2])+
							b4*(sxx[j][i][k+4]-sxx[j][i][k-3])+
							b5*(sxx[j][i][k+5]-sxx[j][i][k-4])+
							b6*(sxx[j][i][k+6]-sxx[j][i][k-5]);

					PML1=((((1.0/DT)-(absorb_coeffx[j][i][k]/2.0)))/(((1.0/DT)+(absorb_coeffx[j][i][k]/2.0))));
					PML2=(1.0/(((1.0/DT)+(absorb_coeffx[j][i][k]/2.0))));

					vx1[j][i][k] = (PML1*vx1[j][i][k]) + PML2*dx*sxx_x/rip;

					PML1=((((1.0/DT)-(absorb_coeffy[j][i][k]/2.0)))/(((1.0/DT)+(absorb_coeffy[j][i][k]/2.0))));
					PML2=(1.0/(((1.0/DT)+(absorb_coeffy[j][i][k]/2.0))));

					vy2[j][i][k] = (PML1*vy2[j][i][k]) + PML2*dy*syy_y/rjp;

					PML1=((((1.0/DT)-(absorb_coeffz[j][i][k]/2.0)))/(((1.0/DT)+(absorb_coeffz[j][i][k]/2.0))));
					PML2=(1.0/(((1.0/DT)+(absorb_coeffz[j][i][k]/2.0))));

					vz3[j][i][k] = (PML1*vz3[j][i][k]) + PML2*dz*szz_z/rkp;

					vx[j][i][k] = vx1[j][i][k];
					vy[j][i][k] = vy2[j][i][k];
					vz[j][i][k] = vz3[j][i][k];

				}
			}
		}
		break;


	default:
		err(" error in particle velocity update: wrong FDORDER. ");
		break;
	} /* end of switch (FDORDER) */


	/* Adding body force components to corresponding particle velocities */
	for (l=1;l<=nsrc;l++) {
		i=(int)srcpos_loc[1][l];
		j=(int)srcpos_loc[2][l];
		k=(int)srcpos_loc[3][l];
		amp=signals[l][nt];

		switch (stype[l]){
		case 2 : 
			vx[j][i][k]+=amp;  /* single force in x */
			break;
		case 3 : 
			vy[j][i][k]+=amp;  /* single force in y */
			break;
		case 4 : 
			vz[j][i][k]+=amp;  /* single force in z */
			break;
		}
	}


	if (LOG)
		if ((MYID==0) && ((nt+(OUTNTIMESTEPINFO-1))%OUTNTIMESTEPINFO)==0){
			time2=MPI_Wtime();
			time=time2-time1;
			fprintf(FP," Real time for particle velocity update: \t %4.2f s.\n",time);
		}
	return time;

}