diff --git a/fdelmodc3D/acoustic4_3D.c b/fdelmodc3D/acoustic4_3D.c
index f3c2c346f37b0b935d0f32ab3d8c466bf3b16ade..36316d652c7b0bbb1c4006368d18d0ad13003a1b 100644
--- a/fdelmodc3D/acoustic4_3D.c
+++ b/fdelmodc3D/acoustic4_3D.c
@@ -6,22 +6,26 @@
 
 #define MIN(x,y) ((x) < (y) ? (x) : (y))
 
-long applySource3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, float *rox, float *roy, float *roz, float *l2m, float **src_nwav, long verbose);
-//long applySource(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long izsrc, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float **src_nwav, long verbose);
+long applySource3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, 
+    float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz,
+    float *rox, float *roy, float *roz, float *l2m, float **src_nwav, long verbose);
 
-long storeSourceOnSurface3D(modPar mod, srcPar src, bndPar bnd, long ixsrc, long iysrc, long izsrc, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, long verbose);
-//long storeSourceOnSurface(modPar mod, srcPar src, bndPar bnd, long ixsrc, long izsrc, float *vx, float *vz, float *tzz, float *txx, float *txz, long verbose);
+long storeSourceOnSurface3D(modPar mod, srcPar src, bndPar bnd, long ixsrc, long iysrc, long izsrc,
+    float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, long verbose);
 
-long reStoreSourceOnSurface3D(modPar mod, srcPar src, bndPar bnd, long ixsrc, long iysrc, long izsrc, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, long verbose);
-//long reStoreSourceOnSurface(modPar mod, srcPar src, bndPar bnd, long ixsrc, long izsrc, float *vx, float *vz, float *tzz, float *txx, float *txz, long verbose);
+long reStoreSourceOnSurface3D(modPar mod, srcPar src, bndPar bnd, long ixsrc, long iysrc, long izsrc,
+    float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, long verbose);
 
-long boundariesP3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, float *rox, float *roy, float *roz, float *l2m, float *lam, float *mul, long itime, long verbose);
-//long boundariesP(modPar mod, bndPar bnd, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float *lam, float *mul, long itime, long verbose);
+long boundariesP3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz,
+    float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, float *rox, float *roy, float *roz,
+    float *l2m, float *lam, float *mul, long itime, long verbose);
 
-long boundariesV3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, float *rox, float *roy, float *roz, float *l2m, float *lam, float *mul, long itime, long verbose);
-//long boundariesV(modPar mod, bndPar bnd, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float *lam, float *mul, long itime, long verbose);
+long boundariesV3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz,
+    float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, float *rox, float *roy, float *roz,
+    float *l2m, float *lam, float *mul, long itime, long verbose);
 
-long acoustic4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, float **src_nwav, float *vx, float *vy, float *vz, float *p, float *rox, float *roy, float *roz, float *l2m, long verbose)
+long acoustic4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, float **src_nwav,
+    float *vx, float *vy, float *vz, float *p, float *rox, float *roy, float *roz, float *l2m, long verbose)
 {
 /*********************************************************************
        COMPUTATIONAL OVERVIEW OF THE 4th ORDER STAGGERED GRID: 
@@ -82,9 +86,8 @@ long acoustic4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, lo
 	/* calculate vx for all grid points except on the virtual boundary*/
 #pragma omp for private (ix, iy, iz) nowait schedule(guided,1)
 	for (iy=mod.ioXy; iy<mod.ieXy; iy++) {
-//#pragma ivdep
         for (ix=mod.ioXx; ix<mod.ieXx; ix++) {
-       #pragma ivdep //deletar 
+       #pragma ivdep 
             for (iz=mod.ioXz; iz<mod.ieXz; iz++) {
                 vx[iy*n2*n1+ix*n1+iz] -= rox[iy*n2*n1+ix*n1+iz]*(
                     c1*(p[iy*n2*n1+ix*n1+iz]     - p[iy*n2*n1+(ix-1)*n1+iz]) +
@@ -96,9 +99,8 @@ long acoustic4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, lo
 	/* calculate vy for all grid points except on the virtual boundary*/
 #pragma omp for private (ix, iy, iz) nowait schedule(guided,1)
 	for (iy=mod.ioYy; iy<mod.ieYy; iy++) {
-//#pragma ivdep
         for (ix=mod.ioYx; ix<mod.ieYx; ix++) {
-            #pragma ivdep //deletar
+            #pragma ivdep
             for (iz=mod.ioYz; iz<mod.ieYz; iz++) {
                 vy[iy*n2*n1+ix*n1+iz] -= roy[iy*n2*n1+ix*n1+iz]*(
                     c1*(p[iy*n2*n1+ix*n1+iz]     - p[(iy-1)*n2*n1+ix*n1+iz]) +
@@ -110,9 +112,8 @@ long acoustic4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, lo
 	/* calculate vz for all grid points except on the virtual boundary */
 #pragma omp for private (ix, iy, iz) schedule(guided,1) 
 	for (iy=mod.ioZy; iy<mod.ieZy; iy++) {
-//#pragma ivdep
         for (ix=mod.ioZx; ix<mod.ieZx; ix++) {
-            #pragma ivdep //deletar
+            #pragma ivdep 
             for (iz=mod.ioZz; iz<mod.ieZz; iz++) {
                 vz[iy*n2*n1+ix*n1+iz] -= roz[iy*n2*n1+ix*n1+iz]*(
                     c1*(p[iy*n2*n1+ix*n1+iz]   - p[iy*n2*n1+ix*n1+iz-1]) +
@@ -142,12 +143,10 @@ long acoustic4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, lo
 
 	/* calculate p/tzz for all grid points except on the virtual boundary */
 #pragma omp for private (ix, iy, iz) schedule(guided,1) 
-//#pragma omp for private (ix, iz) schedule(dynamic) 
 #pragma ivdep
 	for (ix=mod.ioPx; ix<mod.iePx; ix++) {
-//#pragma ivdep
 		for (iy=mod.ioPy; iy<mod.iePy; iy++) {
-            #pragma ivdep //deletar
+            #pragma ivdep
             for (iz=mod.ioPz; iz<mod.iePz; iz++) {
                 p[iy*n1*n2+ix*n1+iz] -= l2m[iy*n1*n2+ix*n1+iz]*(
                             c1*(vx[iy*n1*n2+(ix+1)*n1+iz] - vx[iy*n1*n2+ix*n1+iz]) +
@@ -156,11 +155,6 @@ long acoustic4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, lo
                             c2*(vy[(iy+2)*n1*n2+ix*n1+iz] - vy[(iy-1)*n1*n2+ix*n1+iz]) +
                             c1*(vz[iy*n1*n2+ix*n1+iz+1]   - vz[iy*n1*n2+ix*n1+iz]) +
                             c2*(vz[iy*n1*n2+ix*n1+iz+2]   - vz[iy*n1*n2+ix*n1+iz-1]));
-
-                                //checking //deletar   
-                if(iz==(izsrc+mod.ioPz+bnd.ntap) && ix==(ixsrc+mod.ioPx+bnd.ntap) && iy==(iysrc+mod.ioPy+bnd.ntap) ){
-                    //vmess("inloop--p[ix=%li,iy=%li,iz=%li] at itime %li has value %e", ix, iy, iz, itime, p[iy*n1*n2+ix*n1+iz]); //deletar
-                }
             }
         }
 	}
diff --git a/fdelmodc3D/applySource3D.c b/fdelmodc3D/applySource3D.c
index fa6437dab6aea1d48056c02a62228145728b9816..0a7d1f1165a4eb2d7c7f4c3bdcaa374fd205c1c0 100644
--- a/fdelmodc3D/applySource3D.c
+++ b/fdelmodc3D/applySource3D.c
@@ -20,14 +20,15 @@ void vmess(char *fmt, ...);
  *           The Netherlands 
  *
  **********************************************************************/
-//      applySource3D(mod, src,            wav,              bnd,        itime,      ixsrc,    iysrc,    izsrc,            vx,     vy,           vz,       p,        NULL,         NULL,         NULL,    NULL,       NULL,         rox,       roy,        roz,        l2m,          src_nwav,     verbose);
-long applySource3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, float *rox, float *roy, float *roz, float *l2m, float **src_nwav, long verbose)
+long applySource3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc,
+	float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz,
+	float *rox, float *roy, float *roz, float *l2m, float **src_nwav, long verbose)
 {
 	long is0, ibndz, ibndy, ibndx;
 	long isrc, ix, iy, iz, n1, n2;
 	long id1, id2, id3;
 	float src_ampl, time, scl, dt, sdx;
-	float Mxx, Myy, Mzz, Mxz, Myz, Mxy, rake;
+	float Mxx, Myy, Mzz, Mxz, Myz, Mxy;
 	static long first=1;
 
 	if (src.type==6) {
@@ -350,19 +351,22 @@ long applySource3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, l
                 }
 			}
             else if(src.type == 9) {
-				rake = 0.5*M_PI;
-				Mxx = -1.0*(sin(src.dip)*cos(rake)*sin(2.0*src.strike)+sin(src.dip*2.0)*sin(rake)*sin(src.strike)*sin(src.strike));
-				Mxz = -1.0*(cos(src.dip)*cos(rake)*cos(src.strike)+cos(src.dip*2.0)*sin(rake)*sin(src.strike));
-				Mzz = sin(src.dip*2.0)*sin(rake);
+				Mxx = -1.0*(sin(src.dip)*cos(src.rake)*sin(2.0*src.strike)+sin(src.dip*2.0)*sin(src.rake)*sin(src.strike)*sin(src.strike));
+				Myy = sin(src.dip)*cos(src.rake)*sin(2.0*src.strike)-sin(src.dip*2.0)*sin(src.rake)*cos(src.strike)*cos(src.strike);
+				Mzz = sin(src.dip*2.0)*sin(src.rake);
+				Mxz = -1.0*(cos(src.dip)*cos(src.rake)*cos(src.strike)+cos(src.dip*2.0)*sin(src.rake)*sin(src.strike));
+				Mxy = sin(src.dip)*cos(src.rake)*cos(src.strike*2.0)+0.5*(sin(src.dip*2.0)*sin(src.rake)*sin(src.strike*2.0));
+				Myz = -1.0*(cos(src.dip)*cos(src.rake)*sin(src.strike)-cos(src.dip*2.0)*sin(src.rake)*cos(src.strike));
 
 				txx[iy*n1*n2+ix*n1+iz] -= Mxx*src_ampl;
+				tyy[iy*n1*n2+ix*n1+iz] -= Myy*src_ampl;
 				tzz[iy*n1*n2+ix*n1+iz] -= Mzz*src_ampl;
 				txz[iy*n1*n2+ix*n1+iz] -= Mxz*src_ampl;
+				txy[iy*n1*n2+ix*n1+iz] -= Mxy*src_ampl;
+				tyz[iy*n1*n2+ix*n1+iz] -= Myz*src_ampl;
 			} /* src.type */
 		} /* ischeme */
 	} /* loop over isrc */
 
-
-
 	return 0;
-}
+}
\ No newline at end of file
diff --git a/fdelmodc3D/boundaries3D.c b/fdelmodc3D/boundaries3D.c
index 180e774d6ae6be26f7d07f14bfd3aef92eaeb617..3a3f2cb2f7a7cc245c654ea4163175ab5b0f38d5 100644
--- a/fdelmodc3D/boundaries3D.c
+++ b/fdelmodc3D/boundaries3D.c
@@ -6,7 +6,9 @@
 
 void vmess(char *fmt, ...);
 
-long boundariesP3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, float *rox, float *roy, float *roz, float *l2m, float *lam, float *mul, long itime, long verbose)
+long boundariesP3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz,
+	float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz,
+	float *rox, float *roy, float *roz, float *l2m, float *lam, float *mul, long itime, long verbose)
 {
 /*********************************************************************
 Joeri's original filling order
diff --git a/fdelmodc3D/fdelmodc3D.h b/fdelmodc3D/fdelmodc3D.h
index 158899df6305898453f54fc4344e21bd645b7fe8..1374734c6a2ebbd2c2146d0c3d1ad7d114dcca53 100644
--- a/fdelmodc3D/fdelmodc3D.h
+++ b/fdelmodc3D/fdelmodc3D.h
@@ -176,6 +176,7 @@ typedef struct _sourcePar { /* Source Array Parameters */
 	float amplitude;
 	float dip;
 	float strike;
+	float rake;
 	long distribution;
 	long window;
     long injectionrate;
diff --git a/fdelmodc3D/getParameters3D.c b/fdelmodc3D/getParameters3D.c
index 8771c6bacae14eeafa0e16f643d681555dc454ba..caf7fbb000dc704762ce7c270a4789d02949e670 100644
--- a/fdelmodc3D/getParameters3D.c
+++ b/fdelmodc3D/getParameters3D.c
@@ -540,10 +540,11 @@ criteria we have imposed.*/
 	if (!getparfloat("dyshot",&dyshot)) dyshot=dy;
 	if (!getparfloat("dzshot",&dzshot)) dzshot=0.0;
 	if (!getparfloat("dip",&src->dip)) src->dip=0.0;
-	if (!getparfloat("strike",&src->strike)) src->strike=1.0;
-	if (src->strike>=0) src->strike=0.5*M_PI;
-	else src->strike = -0.5*M_PI;
+	if (!getparfloat("strike",&src->strike)) src->strike=0.0;
+	if (!getparfloat("rake",&src->rake)) src->rake=0.0;
 	src->dip = M_PI*(src->dip/180.0);
+	src->strike = M_PI*(src->strike/180.0);
+	src->rake = M_PI*(src->rake/180.0);
 
 	if (shot->n>1) {
 		idxshot=MAX(0,NINT(dxshot/dx));