diff --git a/fdelmodc3D/Makefile b/fdelmodc3D/Makefile
index 6e103682a7d1295b5c5659b130004210b2e6a498..051ebfb2b300a11d4f05cfb0dc00f7fd00889210 100644
--- a/fdelmodc3D/Makefile
+++ b/fdelmodc3D/Makefile
@@ -25,6 +25,7 @@ SRCC = $(PRG).c \
acoustic6_3D.c \
acousticSH4_3D.c \
acoustic4_qr_3D.c \
+ elastic4dc_3D.c \
defineSource3D.c \
getParameters3D.c \
getWaveletInfo3D.c \
diff --git a/fdelmodc3D/boundaries3D.c b/fdelmodc3D/boundaries3D.c
index 1496fbb7388ce82d83eb4a0913ee289ac8c1b417..8281c66e6db3dc7b1e34c946381aeec9c3711839 100644
--- a/fdelmodc3D/boundaries3D.c
+++ b/fdelmodc3D/boundaries3D.c
@@ -3441,8 +3441,7 @@ MID left mid mid
}
}
-
- if (mod.ischeme <= 2) { /* Elastic scheme */
+ else { /* Elastic scheme */
/* Vx field */
/* right mid mid vx */
diff --git a/fdelmodc3D/fdelmodc3D b/fdelmodc3D/fdelmodc3D
index e8d703ce3c5944cb85865e241beb62f5fdc3123f..ca2229cabd8cf69cb8d1e9832c8f7155b652bae7 100755
Binary files a/fdelmodc3D/fdelmodc3D and b/fdelmodc3D/fdelmodc3D differ
diff --git a/fdelmodc3D/fdelmodc3D.c b/fdelmodc3D/fdelmodc3D.c
index 7cc4e25a275328fa6cf8ebece1facf9bcf041506..705ee6ab9717937cde57cd0ae6de111d94c869eb 100644
--- a/fdelmodc3D/fdelmodc3D.c
+++ b/fdelmodc3D/fdelmodc3D.c
@@ -48,6 +48,11 @@ long acousticSH4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime,
long ixsrc, long iysrc, long izsrc, float **src_nwav, float *tx, float *ty, float *tz,
float *vz, float *rox, float *roy, float *roz, float *mul, long verbose);
+long elastic4dc_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 *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz,
+ float *rox, float *roy, float *roz, float *l2m, float *lam, float *mul, long verbose);
+
long getRecTimes3D(modPar mod, recPar rec, bndPar bnd, long itime, long isam,
float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx,
float *txz, float *txy, float *tyz, float *l2m, float *rox, float *roy, float *roz,
@@ -481,7 +486,6 @@ int main(int argc, char **argv)
shot.z[ishot]=iz+src.sinkdepth;
}
-
/* scan for free surface boundary in case it has a topography */
for (iy=0; iy<mod.ny; iy++) {
for (ix=0; ix<mod.nx; ix++) {
@@ -631,7 +635,8 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, iysrc, izsrc, it, src_nwav, verbos
vmess("Visco-Elastic order 4 not yet available");
break;
case 5 : /* Elastic FD kernel with S-velocity set to zero*/
- vmess("DC-Elastic order 4 not yet available");
+ elastic4dc_3D(mod, src, wav, bnd, it, ixsrc, iysrc, izsrc, src_nwav,
+ vx, vy, vz, tzz, tyy, txx, txz, txy, tyz, rox, roy, roz, l2m, lam, mul, verbose);
break;
}
diff --git a/fdelmodc3D/getParameters3D.c b/fdelmodc3D/getParameters3D.c
index caf7fbb000dc704762ce7c270a4789d02949e670..bce77f757574fdb37bef38104cf88b27967af84b 100644
--- a/fdelmodc3D/getParameters3D.c
+++ b/fdelmodc3D/getParameters3D.c
@@ -959,6 +959,7 @@ criteria we have imposed.*/
vmess("*******************************************");
vmess("wav_nt = %6li wav_nx = %li", wav->ns, wav->nx);
vmess("src_type = %6li src_orient = %li", src->type, src->orient);
+ vmess("number of sources = %li", shot->n);
vmess("fmax = %8.2f", fmax);
fprintf(stderr," %s: Source type : ",xargv[0]);
switch ( src->type ) {
diff --git a/fdelmodc3D/getRecTimes3D.c b/fdelmodc3D/getRecTimes3D.c
index fa007f88322b1cb038b7d2d820218d118ffe7835..4c19786a8028eb10c0bc3183fe3c31c72e3603e8 100644
--- a/fdelmodc3D/getRecTimes3D.c
+++ b/fdelmodc3D/getRecTimes3D.c
@@ -402,6 +402,7 @@ long getRecTimes3D(modPar mod, recPar rec, bndPar bnd, long itime, long isam,
}
}
if (rec.type.txx) rec_txx[irec*rec.nt+isam] = txx[iy*n1*n2+ix*n1+iz];
+ if (rec.type.tyy) rec_tyy[irec*rec.nt+isam] = tyy[iy*n1*n2+ix*n1+iz];
if (rec.type.tzz) rec_tzz[irec*rec.nt+isam] = tzz[iy*n1*n2+ix*n1+iz];
if (rec.type.txz) { /* time interpolation to be done */
if (rec.int_vz == 2 || rec.int_vx == 2) {
diff --git a/fdelmodc3D/readModel3D.c b/fdelmodc3D/readModel3D.c
index 024e279289a74f6ac8a68960a4058fb0f28539ac..157f7a1790d4bcb1c8539aebf37385b170679e63 100644
--- a/fdelmodc3D/readModel3D.c
+++ b/fdelmodc3D/readModel3D.c
@@ -27,7 +27,7 @@
long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz,
- float *l2m, float *lam, float *muu, float *tss, float *tes, float *tep)
+ float *l2m, float *lam, float *muxz, float *tss, float *tes, float *tep)
{
FILE *fpcp, *fpcs, *fpro;
FILE *fpqp=NULL, *fpqs=NULL;
@@ -36,7 +36,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz,
long n1, n2, n3, ix, iy, iz, nz, ny, nx;
long ixo, iyo, izo, ixe, iye, ize;
long ioXx, ioXy, ioXz, ioYx, ioYy, ioYz, ioZz, ioZy, ioZx, ioPx, ioPy, ioPz, ioTx, ioTy, ioTz;
- float cp2, cs2, cs11, cs12, cs21, cs22, mul, mu, lamda2mu, lamda;
+ float cp2, cs2, cs111, cs112, cs121, cs211, cs122, cs212, cs221, mul, mu, lamda2mu, lamda;
float cs2c, cs2b, cs2a, cpx, cpy, cpz, bx, by, bz, fac;
float *cp, *cs, *ro, *qp, *qs;
float a, b;
@@ -69,7 +69,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz,
ioPx=mod.ioPx;
ioPy=mod.ioPy;
ioPz=mod.ioPz;
- /* Txz, Txy, Tyz,: muu */
+ /* Txz, Txy, Tyz,: muxz */
ioTx=mod.ioTx;
ioTy=mod.ioTy;
ioTz=mod.ioTz;
@@ -231,23 +231,263 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz,
/* the edges of the model */
if (mod.ischeme>2) { /* Elastic Scheme */
+ iz = nz-1;
+ for (iy=0;iy<ny-1;iy++) {
+ for (ix=0;ix<nx-1;ix++) {
+ /* for muxz field */
+ /* csxyz */
+ cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz];
+ cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz];
+ cs2a = cs[iy*nx*nz+(ix+1)*nz+iz]*cs[iy*nx*nz+(ix+1)*nz+iz];
+ cs111 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ cs112 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ cs211 = cs2a*ro[iy*nx*nz+(ix+1)*nz+iz];
+ cs212 = cs2a*ro[iy*nx*nz+(ix+1)*nz+iz];
+ if (cs111 > 0.0) {
+ mul = 4.0/(1.0/cs111+1.0/cs112+1.0/cs211+1.0/cs212);
+ }
+ else {
+ mul = 0.0;
+ }
+
+ /* for muyz field IN PROGRESS!!!!!!!!!!!!!!!!! */
+ /* csxyz */
+ // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz];
+ // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz];
+ // cs2a = cs[(iy+1)*nx*nz+ix*nz+iz]*cs[(iy+1)*nx*nz+ix*nz+iz];
+ // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // cs121 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // cs112 = cs2a*ro[(iy+1)*nx*nz+ix*nz+iz];
+ // cs122 = cs2a*ro[(iy+1)*nx*nz+ix*nz+iz];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs121+1.0/cs112+1.0/cs122);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
+ /* for muxy field IN PROGRESS!!!!!!!!!!!!!!!!! */
+ /* csxyz */
+ // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz];
+ // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz];
+ // cs2a = cs[(iy+1)*nx*nz+ix*nz+iz]*cs[(iy+1)*nx*nz+ix*nz+iz];
+ // cs2b = cs[iy*nx*nz+(ix+1)*nz+iz]*cs[iy*nx*nz+(ix+1)*nz+iz];
+ // cs2c = cs[(iy+1)*nx*nz+(ix+1)*nz+iz]*cs[(iy+1)*nx*nz+(ix+1)*nz+iz];
+ // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // cs121 = cs2a*ro[(iy+1)*nx*nz+ix*nz+iz];
+ // cs211 = cs2b*ro[iy*nx*nz+(ix+1)*nz+iz];
+ // cs221 = cs2c*ro[(iy+1)*nx*nz+(ix+1)*nz+iz];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs211+1.0/cs121+1.0/cs221);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
+ mu = cs2*ro[iy*nx*nz+ix*nz+iz];
+ lamda2mu = cp2*ro[iy*nx*nz+ix*nz+iz];
+ lamda = lamda2mu - 2*mu;
+
+ bx = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[iy*nx*nz+(ix+1)*nz+iz]);
+ by = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[(iy+1)*nx*nz+ix*nz+iz]);
+ bz = ro[iy*nx*nz+ix*nz+iz];
+ rox[(iy+ioXy)*n2*n1+(ix+ioXx)*n1+iz+ioXz]=fac/bx;
+ roy[(iy+ioYy)*n2*n1+(ix+ioYx)*n1+iz+ioYz]=fac/by;
+ roz[(iy+ioZy)*n2*n1+(ix+ioZx)*n1+iz+ioZz]=fac/bz;
+ l2m[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu;
+ lam[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda;
+ muxz[(iy+ioTy)*n2*n1+(ix+ioTx)*n1+iz+ioTz]=fac*mul;
+ }
+ }
+
+ iy = ny-1;
+ for (ix=0;ix<nx-1;ix++) {
+ for (iz=0;iz<nz-1;iz++) {
+ /* for muxz field */
+ /* csxyz */
+ cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz];
+ cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz];
+ cs2a = cs[iy*nx*nz+(ix+1)*nz+iz]*cs[iy*nx*nz+(ix+1)*nz+iz];
+ cs2b = cs[iy*nx*nz+ix*nz+iz+1]*cs[iy*nx*nz+ix*nz+iz+1];
+ cs2c = cs[iy*nx*nz+(ix+1)*nz+iz+1]*cs[iy*nx*nz+(ix+1)*nz+iz+1];
+ cs111 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ cs112 = cs2b*ro[iy*nx*nz+ix*nz+iz+1];
+ cs211 = cs2a*ro[iy*nx*nz+(ix+1)*nz+iz];
+ cs212 = cs2c*ro[iy*nx*nz+(ix+1)*nz+iz+1];
+ if (cs111 > 0.0) {
+ mul = 4.0/(1.0/cs111+1.0/cs112+1.0/cs211+1.0/cs212);
+ }
+ else {
+ mul = 0.0;
+ }
+
+ /* for muyz field IN PROGRESS!!!!!!!!!!!!!!!!! */
+ /* csxyz */
+ // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz];
+ // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz];
+ // cs2a = cs[iy*nx*nz+ix*nz+iz+1]*cs[iy*nx*nz+ix*nz+iz+1];
+ // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // cs121 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // cs112 = cs2a*ro[iy*nx*nz+ix*nz+iz+1];
+ // cs122 = cs2a*ro[iy*nx*nz+ix*nz+iz+1];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs121+1.0/cs112+1.0/cs122);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
+ /* for muxy field IN PROGRESS!!!!!!!!!!!!!!!!! */
+ /* csxyz */
+ // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz];
+ // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz];
+ // cs2a = cs[iy*nx*nz+(ix+1)*nz+iz]*cs[iy*nx*nz+(ix+1)*nz+iz];
+ // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // cs121 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // cs211 = cs2a*ro[iy*nx*nz+(ix+1)*nz+iz];
+ // cs221 = cs2a*ro[iy*nx*nz+(ix+1)*nz+iz];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs211+1.0/cs121+1.0/cs221);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
+ mu = cs2*ro[iy*nx*nz+ix*nz+iz];
+ lamda2mu = cp2*ro[iy*nx*nz+ix*nz+iz];
+ lamda = lamda2mu - 2*mu;
+
+ bx = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[iy*nx*nz+(ix+1)*nz+iz]);
+ by = ro[iy*nx*nz+ix*nz+iz];
+ bz = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[iy*nx*nz+ix*nz+iz+1]);
+ rox[(iy+ioXy)*n2*n1+(ix+ioXx)*n1+iz+ioXz]=fac/bx;
+ roy[(iy+ioYy)*n2*n1+(ix+ioYx)*n1+iz+ioYz]=fac/by;
+ roz[(iy+ioZy)*n2*n1+(ix+ioZx)*n1+iz+ioZz]=fac/bz;
+ l2m[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu;
+ lam[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda;
+ muxz[(iy+ioTy)*n2*n1+(ix+ioTx)*n1+iz+ioTz]=fac*mul;
+ }
+ }
+
+ ix = nx-1;
+ for (iy=0;iy<ny-1;iy++) {
+ for (iz=0;iz<nz-1;iz++) {
+ /* for muxz field */
+ /* csxyz */
+ cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz];
+ cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz];
+ cs2a = cs[iy*nx*nz+ix*nz+iz+1]*cs[iy*nx*nz+ix*nz+iz+1];
+ cs111 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ cs112 = cs2a*ro[iy*nx*nz+ix*nz+iz+1];
+ cs211 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ cs212 = cs2a*ro[iy*nx*nz+ix*nz+iz+1];
+ if (cs111 > 0.0) {
+ mul = 4.0/(1.0/cs111+1.0/cs112+1.0/cs211+1.0/cs212);
+ }
+ else {
+ mul = 0.0;
+ }
+
+ /* for muyz field IN PROGRESS!!!!!!!!!!!!!!!!! */
+ /* csxyz */
+ // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz];
+ // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz];
+ // cs2a = cs[iy*nx*nz+ix*nz+iz+1]*cs[iy*nx*nz+ix*nz+iz+1];
+ // cs2b = cs[(iy+1)*nx*nz+ix*nz+iz]*cs[(iy+1)*nx*nz+ix*nz+iz];
+ // cs2c = cs[(iy+1)*nx*nz+ix*nz+iz+1]*cs[(iy+1)*nx*nz+ix*nz+iz+1];
+ // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // cs121 = cs2b*ro[(iy+1)*nx*nz+ix*nz+iz];
+ // cs112 = cs2a*ro[iy*nx*nz+ix*nz+iz+1];
+ // cs122 = cs2c*ro[(iy+1)*nx*nz+ix*nz+iz+1];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs121+1.0/cs112+1.0/cs122);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
+ /* for muxy field IN PROGRESS!!!!!!!!!!!!!!!!! */
+ /* csxyz */
+ // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz];
+ // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz];
+ // cs2a = cs[(iy+1)*nx*nz+ix*nz+iz]*cs[(iy+1)*nx*nz+ix*nz+iz];
+ // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // cs121 = cs2a*ro[(iy+1)*nx*nz+ix*nz+iz];
+ // cs211 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // cs221 = cs2a*ro[(iy+1)*nx*nz+ix*nz+iz];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs211+1.0/cs121+1.0/cs221);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
+ mu = cs2*ro[iy*nx*nz+ix*nz+iz];
+ lamda2mu = cp2*ro[iy*nx*nz+ix*nz+iz];
+ lamda = lamda2mu - 2*mu;
+
+ bx = ro[iy*nx*nz+ix*nz+iz];
+ by = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[(iy+1)*nx*nz+ix*nz+iz]);
+ bz = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[iy*nx*nz+ix*nz+iz+1]);
+ rox[(iy+ioXy)*n2*n1+(ix+ioXx)*n1+iz+ioXz]=fac/bx;
+ roy[(iy+ioYy)*n2*n1+(ix+ioYx)*n1+iz+ioYz]=fac/by;
+ roz[(iy+ioZy)*n2*n1+(ix+ioZx)*n1+iz+ioZz]=fac/bz;
+ l2m[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu;
+ lam[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda;
+ muxz[(iy+ioTy)*n2*n1+(ix+ioTx)*n1+iz+ioTz]=fac*mul;
+ }
+ }
+
iz = nz-1;
iy = ny-1;
for (ix=0;ix<nx-1;ix++) {
+ /* for muxz field */
+ /* csxyz */
cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz];
cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz];
cs2a = cs[iy*nx*nz+(ix+1)*nz+iz]*cs[iy*nx*nz+(ix+1)*nz+iz];
- cs11 = cs2*ro[iy*nx*nz+ix*nz+iz];
- cs12 = cs2*ro[iy*nx*nz+ix*nz+iz];
- cs21 = cs2a*ro[iy*nx*nz+(ix+1)*nz+iz];
- cs22 = cs2a*ro[iy*nx*nz+(ix+1)*nz+iz];
-
- if (cs11 > 0.0) {
- mul = 4.0/(1.0/cs11+1.0/cs12+1.0/cs21+1.0/cs22);
+ cs111 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ cs112 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ cs211 = cs2a*ro[iy*nx*nz+(ix+1)*nz+iz];
+ cs212 = cs2a*ro[iy*nx*nz+(ix+1)*nz+iz];
+ if (cs111 > 0.0) {
+ mul = 4.0/(1.0/cs111+1.0/cs112+1.0/cs211+1.0/cs212);
}
else {
mul = 0.0;
}
+
+ /* for muyz field IN PROGRESS!!!!!!!!!!!!!!!!! */
+ /* csxyz */
+ // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz];
+ // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz];
+ // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // cs121 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // cs112 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // cs122 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs121+1.0/cs112+1.0/cs122);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
+ /* for muxy field IN PROGRESS!!!!!!!!!!!!!!!!! */
+ /* csxyz */
+ // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz];
+ // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz];
+ // cs2a = cs[iy*nx*nz+(ix+1)*nz+iz]*cs[iy*nx*nz+(ix+1)*nz+iz];
+ // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // cs121 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // cs211 = cs2a*ro[iy*nx*nz+(ix+1)*nz+iz];
+ // cs221 = cs2a*ro[iy*nx*nz+(ix+1)*nz+iz];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs211+1.0/cs121+1.0/cs221);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
mu = cs2*ro[iy*nx*nz+ix*nz+iz];
lamda2mu = cp2*ro[iy*nx*nz+ix*nz+iz];
lamda = lamda2mu - 2*mu;
@@ -260,58 +500,124 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz,
roz[(iy+ioZy)*n2*n1+(ix+ioZx)*n1+iz+ioZz]=fac/bz;
l2m[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu;
lam[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda;
- muu[(iy+ioTy)*n2*n1+(ix+ioTx)*n1+iz+ioTz]=fac*mul;
+ muxz[(iy+ioTy)*n2*n1+(ix+ioTx)*n1+iz+ioTz]=fac*mul;
}
ix = nx-1;
iz = nz-1;
for (iy=0;iy<ny-1;iy++) {
+ /* for muxz field */
+ /* csxyz */
cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz];
cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz];
- cs2b = cs[(iy+1)*nx*nz+ix*nz+iz]*cs[(iy+1)*nx*nz+ix*nz+iz];
- cs11 = cs2*ro[iy*nx*nz+ix*nz+iz];
- cs12 = cs2b*ro[iy*nx*nz+ix*nz+iz];
- cs21 = cs2*ro[iy*nx*nz+ix*nz+iz];
- cs22 = cs2b*ro[iy*nx*nz+ix*nz+iz];
-
- if (cs11 > 0.0) {
- mul = 4.0/(1.0/cs11+1.0/cs12+1.0/cs21+1.0/cs22);
+ cs111 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ cs112 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ cs211 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ cs212 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ if (cs111 > 0.0) {
+ mul = 4.0/(1.0/cs111+1.0/cs112+1.0/cs211+1.0/cs212);
}
else {
mul = 0.0;
}
+
+ /* for muyz field IN PROGRESS!!!!!!!!!!!!!!!!! */
+ /* csxyz */
+ // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz];
+ // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz];
+ // cs2a = cs[(iy+1)*nx*nz+ix*nz+iz]*cs[(iy+1)*nx*nz+ix*nz+iz];
+ // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // cs121 = cs2a*ro[(iy+1)*nx*nz+ix*nz+iz];
+ // cs112 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // cs122 = cs2a*ro[(iy+1)*nx*nz+ix*nz+iz];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs121+1.0/cs112+1.0/cs122);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
+ /* for muxy field IN PROGRESS!!!!!!!!!!!!!!!!! */
+ /* csxyz */
+ // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz];
+ // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz];
+ // cs2a = cs[(iy+1)*nx*nz+ix*nz+iz]*cs[(iy+1)*nx*nz+ix*nz+iz];
+ // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // cs121 = cs2a*ro[(iy+1)*nx*nz+ix*nz+iz];
+ // cs211 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // cs221 = cs2a*ro[(iy+1)*nx*nz+ix*nz+iz];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs211+1.0/cs121+1.0/cs221);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
mu = cs2*ro[iy*nx*nz+ix*nz+iz];
lamda2mu = cp2*ro[iy*nx*nz+ix*nz+iz];
lamda = lamda2mu - 2*mu;
bx = ro[iy*nx*nz+ix*nz+iz];
- by = ro[iy*nx*nz+ix*nz+iz];
- bz = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[iy*nx*nz+ix*nz+iz+1]);
+ by = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[(iy+1)*nx*nz+ix*nz+iz]);
+ bz = ro[iy*nx*nz+ix*nz+iz];
rox[(iy+ioXy)*n2*n1+(ix+ioXx)*n1+iz+ioXz]=fac/bx;
roy[(iy+ioYy)*n2*n1+(ix+ioYx)*n1+iz+ioYz]=fac/bx;
roz[(iy+ioZy)*n2*n1+(ix+ioZx)*n1+iz+ioZz]=fac/bz;
l2m[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu;
lam[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda;
- muu[(iy+ioTy)*n2*n1+(ix+ioTx)*n1+iz+ioTz]=fac*mul;
+ muxz[(iy+ioTy)*n2*n1+(ix+ioTx)*n1+iz+ioTz]=fac*mul;
}
ix = nx-1;
iy = ny-1;
for (iz=0;iz<nz-1;iz++) {
+ /* for muxz field */
+ /* csxyz */
cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz];
cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz];
- cs2b = cs[iy*nx*nz+ix*nz+iz+1]*cs[iy*nx*nz+ix*nz+iz+1];
- cs11 = cs2*ro[iy*nx*nz+ix*nz+iz];
- cs12 = cs2b*ro[iy*nx*nz+ix*nz+iz+1];
- cs21 = cs2*ro[iy*nx*nz+ix*nz+iz];
- cs22 = cs2b*ro[iy*nx*nz+ix*nz+iz+1];
-
- if (cs11 > 0.0) {
- mul = 4.0/(1.0/cs11+1.0/cs12+1.0/cs21+1.0/cs22);
+ cs2a = cs[iy*nx*nz+ix*nz+iz+1]*cs[iy*nx*nz+ix*nz+iz+1];
+ cs111 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ cs112 = cs2a*ro[iy*nx*nz+ix*nz+iz+1];
+ cs211 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ cs212 = cs2a*ro[iy*nx*nz+ix*nz+iz+1];
+ if (cs111 > 0.0) {
+ mul = 4.0/(1.0/cs111+1.0/cs112+1.0/cs211+1.0/cs212);
}
else {
mul = 0.0;
}
+
+ /* for muyz field IN PROGRESS!!!!!!!!!!!!!!!!! */
+ /* csxyz */
+ // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz];
+ // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz];
+ // cs2a = cs[iy*nx*nz+ix*nz+iz+1]*cs[iy*nx*nz+ix*nz+iz+1];
+ // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // cs121 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // cs112 = cs2a*ro[iy*nx*nz+ix*nz+iz+1];
+ // cs122 = cs2a*ro[iy*nx*nz+ix*nz+iz+1];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs121+1.0/cs112+1.0/cs122);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
+ /* for muxy field IN PROGRESS!!!!!!!!!!!!!!!!! */
+ /* csxyz */
+ // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz];
+ // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz];
+ // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // cs121 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // cs211 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // cs221 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs211+1.0/cs121+1.0/cs221);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
mu = cs2*ro[iy*nx*nz+ix*nz+iz];
lamda2mu = cp2*ro[iy*nx*nz+ix*nz+iz];
lamda = lamda2mu - 2*mu;
@@ -324,7 +630,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz,
roz[(iy+ioZy)*n2*n1+(ix+ioZx)*n1+iz+ioZz]=fac/bz;
l2m[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu;
lam[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda;
- muu[(iy+ioTy)*n2*n1+(ix+ioTx)*n1+iz+ioTz]=fac*mul;
+ muxz[(iy+ioTy)*n2*n1+(ix+ioTx)*n1+iz+ioTz]=fac*mul;
}
ix=nx-1;
@@ -339,52 +645,85 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz,
by = ro[iy*nx*nz+ix*nz+iz];
bz = ro[iy*nx*nz+ix*nz+iz];
rox[(iy+ioXy)*n2*n1+(ix+ioXx)*n1+iz+ioXz]=fac/bx;
- roz[(iy+ioYy)*n2*n1+(ix+ioYx)*n1+iz+ioYz]=fac/by;
+ roy[(iy+ioYy)*n2*n1+(ix+ioYx)*n1+iz+ioYz]=fac/by;
roz[(iy+ioZy)*n2*n1+(ix+ioZx)*n1+iz+ioZz]=fac/bz;
- l2m[(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu;
- lam[(ix+ioPx)*n1+iz+ioPz]=fac*lamda;
- muu[(ix+ioTx)*n1+iz+ioTz]=fac*mu;
+ l2m[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu;
+ lam[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda;
+ muxz[(iy+ioTy)*n2*n1+(ix+ioTx)*n1+iz+ioTz]=fac*mu;
+
+ for (iy=0;iy<ny-1;iy++) {
+ for (ix=0;ix<nx-1;ix++) {
+ for (iz=0;iz<nz-1;iz++) {
+ /* for muxz field */
+ /* csxyz */
+ cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz];
+ cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz];
+ cs2a = cs[iy*nx*nz+(ix+1)*nz+iz]*cs[iy*nx*nz+(ix+1)*nz+iz];
+ cs2b = cs[iy*nx*nz+ix*nz+iz+1]*cs[iy*nx*nz+ix*nz+iz+1];
+ cs2c = cs[iy*nx*nz+(ix+1)*nz+iz+1]*cs[iy*nx*nz+(ix+1)*nz+iz+1];
+ cs111 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ cs112 = cs2b*ro[iy*nx*nz+ix*nz+iz+1];
+ cs211 = cs2a*ro[iy*nx*nz+ix*nz+iz];
+ cs212 = cs2c*ro[iy*nx*nz+ix*nz+iz+1];
+ if (cs111 > 0.0) {
+ mul = 4.0/(1.0/cs111+1.0/cs112+1.0/cs211+1.0/cs212);
+ }
+ else {
+ mul = 0.0;
+ }
- for (ix=0;ix<nx-1;ix++) {
- for (iz=0;iz<nz-1;iz++) {
- cp2 = cp[ix*nz+iz]*cp[ix*nz+iz];
- cs2 = cs[ix*nz+iz]*cs[ix*nz+iz];
- cs2a = cs[(ix+1)*nz+iz]*cs[(ix+1)*nz+iz];
- cs2b = cs[ix*nz+iz+1]*cs[ix*nz+iz+1];
- cs2c = cs[(ix+1)*nz+iz+1]*cs[(ix+1)*nz+iz+1];
-
-/*
-Compute harmonic average of mul for accurate and stable fluid-solid interface
-see Finite-difference modeling of wave propagation in a fluid-solid configuration
-Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
-*/
-
- cs11 = cs2*ro[ix*nz+iz];
- cs12 = cs2b*ro[ix*nz+iz+1];
- cs21 = cs2a*ro[ix*nz+iz];
- cs22 = cs2c*ro[ix*nz+iz+1];
-// cpx = 0.5*(cp[ix*nz+iz]+cp[(ix+1)*nz+iz])
-// cpz = 0.5*(cp[ix*nz+iz]+cp[ix*nz+iz+1])
-
- if (cs11 > 0.0) {
- mul = 4.0/(1.0/cs11+1.0/cs12+1.0/cs21+1.0/cs22);
- }
- else {
- mul = 0.0;
- }
- mu = cs2*ro[ix*nz+iz];
- lamda2mu = cp2*ro[ix*nz+iz];
- lamda = lamda2mu - 2*mu; /* could also use mul to calculate lambda, but that might not be correct: question from Chaoshun Hu. Note use mu or mul as well on boundaries */
-
- bx = 0.5*(ro[ix*nz+iz]+ro[(ix+1)*nz+iz]);
- bz = 0.5*(ro[ix*nz+iz]+ro[ix*nz+iz+1]);
- rox[(ix+ioXx)*n1+iz+ioXz]=fac/bx;
- roz[(ix+ioZx)*n1+iz+ioZz]=fac/bz;
- l2m[(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu;
- lam[(ix+ioPx)*n1+iz+ioPz]=fac*lamda;
- muu[(ix+ioTx)*n1+iz+ioTz]=fac*mul;
- }
- }
+ /* for muyz field IN PROGRESS!!!!!!!!!!!!!!!!! */
+ /* csxyz */
+ // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz];
+ // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz];
+ // cs2a = cs[iy*nx*nz+ix*nz+iz+1]*cs[iy*nx*nz+ix*nz+iz+1];
+ // cs2b = cs[(iy+1)*nx*nz+ix*nz+iz]*cs[(iy+1)*nx*nz+ix*nz+iz];
+ // cs2c = cs[(iy+1)*nx*nz+ix*nz+iz+1]*cs[(iy+1)*nx*nz+ix*nz+iz+1];
+ // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // cs121 = cs2b*ro[(iy+1)*nx*nz+ix*nz+iz];
+ // cs112 = cs2a*ro[iy*nx*nz+ix*nz+iz+1];
+ // cs122 = cs2c*ro[(iy+1)*nx*nz+ix*nz+iz+1];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs121+1.0/cs112+1.0/cs122);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
+ /* for muxy field IN PROGRESS!!!!!!!!!!!!!!!!! */
+ /* csxyz */
+ // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz];
+ // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz];
+ // cs2a = cs[iy*nx*nz+(ix+1)*nz+iz]*cs[iy*nx*nz+(ix+1)*nz+iz];
+ // cs2b = cs[(iy+1)*nx*nz+ix*nz+iz]*cs[(iy+1)*nx*nz+ix*nz+iz];
+ // cs2c = cs[(iy+1)*nx*nz+(ix+1)*nz+iz]*cs[(iy+1)*nx*nz+(ix+1)*nz+iz];
+ // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz];
+ // cs121 = cs2b*ro[(iy+1)*nx*nz+ix*nz+iz];
+ // cs211 = cs2a*ro[iy*nx*nz+(ix+1)*nz+iz];
+ // cs221 = cs2c*ro[(iy+1)*nx*nz+(ix+1)*nz+iz];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs211+1.0/cs121+1.0/cs221);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
+ mu = cs2*ro[iy*nx*nz+ix*nz+iz];
+ lamda2mu = cp2*ro[iy*nx*nz+ix*nz+iz];
+ lamda = lamda2mu - 2*mu;
+
+ bx = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[iy*nx*nz+(ix+1)*nz+iz]);
+ by = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[(iy+1)*nx*nz+ix*nz+iz]);
+ bz = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[iy*nx*nz+ix*nz+iz+1]);
+ rox[(iy+ioXy)*n2*n1+(ix+ioXx)*n1+iz+ioXz]=fac/bx;
+ roy[(iy+ioYy)*n2*n1+(ix+ioYx)*n1+iz+ioYz]=fac/by;
+ roz[(iy+ioZy)*n2*n1+(ix+ioZx)*n1+iz+ioZz]=fac/bz;
+ l2m[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu;
+ lam[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda;
+ muxz[(iy+ioTy)*n2*n1+(ix+ioTx)*n1+iz+ioTz]=fac*mul;
+ }
+ }
+ }
}
else { /* Acoustic Scheme */
@@ -610,7 +949,7 @@ Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
}
}
}
- /* muu field */
+ /* muxz field */
ixo = mod.ioTx;
ixe = mod.ioTx+bnd.ntap;
iyo = mod.ioTy;
@@ -620,7 +959,7 @@ Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
for (iy=iyo; iy<iye; iy++) {
for (ix=ixo; ix<ixe; ix++) {
for (iz=izo; iz<ize; iz++) {
- muu[iy*n2*n1+ix*n1+iz] = muu[iy*n2*n1+ixe*n1+iz];
+ muxz[iy*n2*n1+ix*n1+iz] = muxz[iy*n2*n1+ixe*n1+iz];
}
}
}
@@ -741,7 +1080,7 @@ Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
}
}
- /* muu field */
+ /* muxz field */
ixo = mod.ieTx-bnd.ntap;
ixe = mod.ieTx;
iyo = mod.ioTy;
@@ -751,7 +1090,7 @@ Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
for (iy=iyo; iy<iye; iy++) {
for (ix=ixo; ix<ixe; ix++) {
for (iz=izo; iz<ize; iz++) {
- muu[iy*n2*n1+ix*n1+iz] = muu[iy*n2*n1+(ixo-1)*n1+iz];
+ muxz[iy*n2*n1+ix*n1+iz] = muxz[iy*n2*n1+(ixo-1)*n1+iz];
}
}
}
@@ -883,7 +1222,7 @@ Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
}
}
}
- /* muu field */
+ /* muxz field */
ixo = mod.ioTx;
ixe = mod.ieTx;
if (bnd.lef==4 || bnd.lef==2) ixo -= bnd.ntap;
@@ -895,7 +1234,7 @@ Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
for (iy=iyo; iy<iye; iy++) {
for (ix=ixo; ix<ixe; ix++) {
for (iz=izo; iz<ize; iz++) {
- muu[iy*n2*n1+ix*n1+iz] = muu[iye*n2*n1+ix*n1+iz];
+ muxz[iy*n2*n1+ix*n1+iz] = muxz[iye*n2*n1+ix*n1+iz];
}
}
}
@@ -1032,7 +1371,7 @@ Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
}
}
- /* muu field */
+ /* muxz field */
ixo = mod.ioTx;
ixe = mod.ieTx;
if (bnd.lef==4 || bnd.lef==2) ixo -= bnd.ntap;
@@ -1044,7 +1383,7 @@ Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
for (iy=iyo; iy<iye; iy++) {
for (ix=ixo; ix<ixe; ix++) {
for (iz=izo; iz<ize; iz++) {
- muu[iy*n2*n1+ix*n1+iz] = muu[(iyo-1)*n2*n1+ix*n1+iz];
+ muxz[iy*n2*n1+ix*n1+iz] = muxz[(iyo-1)*n2*n1+ix*n1+iz];
}
}
}
@@ -1181,7 +1520,7 @@ Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
}
}
- /* muu field */
+ /* muxz field */
ixo = mod.ioTx;
ixe = mod.ieTx;
iyo = mod.ioTy;
@@ -1191,7 +1530,7 @@ Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
for (iy=iyo; iy<iye; iy++) {
for (ix=ixo; ix<ixe; ix++) {
for (iz=izo; iz<ize; iz++) {
- muu[iy*n2*n1+ix*n1+iz] = muu[iy*n2*n1+ix*n1+ize];
+ muxz[iy*n2*n1+ix*n1+iz] = muxz[iy*n2*n1+ix*n1+ize];
}
}
}
@@ -1323,7 +1662,7 @@ Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
}
}
- /* muu */
+ /* muxz */
ixo = mod.ioTx;
ixe = mod.ieTx;
iyo = mod.ioTy;
@@ -1333,7 +1672,7 @@ Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
for (iy=iyo; iy<iye; iy++) {
for (ix=ixo; ix<ixe; ix++) {
for (iz=izo; iz<ize; iz++) {
- muu[iy*n2*n1+ix*n1+iz] = muu[iy*n2*n1+ix*n1+izo-1];
+ muxz[iy*n2*n1+ix*n1+iz] = muxz[iy*n2*n1+ix*n1+izo-1];
}
}
}