diff --git a/marchenko/applyMute.c b/marchenko/applyMute.c
index 64f806849ee753a7da99d4c0dc7807f7bec6787e..881d1779ad3834836a924b2668517c091d2af6c8 100644
--- a/marchenko/applyMute.c
+++ b/marchenko/applyMute.c
@@ -20,7 +20,7 @@ void applyMute( float *data, int *mute, int smooth, int above, int Nfoc, int nxs
if (smooth) {
costaper = (float *)malloc(smooth*sizeof(float));
- scl = M_PI/((float)smooth);
+ scl = M_PI/((float)(smooth+1));
for (i=0; i<smooth; i++) {
costaper[i] = 0.5*(1.0+cos((i+1)*scl));
}
@@ -40,7 +40,7 @@ void applyMute( float *data, int *mute, int smooth, int above, int Nfoc, int nxs
}
}
}
- else if (above==0){
+ else if (above==0){ /* implementation of <t1-epsilon : -t1+epsilon> window */
for (i = 0; i < npos; i++) {
imute = ixpos[i];
tmute = mute[isyn*nxs+imute];
@@ -49,18 +49,36 @@ void applyMute( float *data, int *mute, int smooth, int above, int Nfoc, int nxs
memset(&data[isyn*nxs*nt+i*nt],0, sizeof(float)*nt);
continue;
}
- for (j = MAX(0,-2*ts+tmute-shift),l=0; j < MAX(0,-2*ts+tmute-shift+smooth); j++,l++) {
+ for (j = MAX(0,-2*ts+tmute-shift-smooth),l=0; j < MAX(0,-2*ts+tmute-shift); j++,l++) {
data[isyn*nxs*nt+i*nt+j] *= costaper[l];
}
- for (j = MAX(0,-2*ts+tmute-shift+smooth)+1; j < MIN(nt,nt+1-tmute+shift-smooth); j++) {
+ for (j = MAX(0,-2*ts+tmute-shift); j < MIN(nt,nt-tmute+shift); j++) {
data[isyn*nxs*nt+i*nt+j] = 0.0;
}
- for (j = MIN(nt,nt-tmute+shift-smooth),l=0; j < MIN(nt,nt-tmute+shift); j++,l++) {
+ for (j = MIN(nt,nt-tmute+shift),l=0; j < MIN(nt,nt-tmute+shift+smooth); j++,l++) {
data[isyn*nxs*nt+i*nt+j] *= costaper[smooth-l-1];
}
}
}
- else if (above==-2){
+ else if (above==4 || above==-4) { //Psi gate which is the inverse of the Theta gate (above=0)
+ for (i = 0; i < npos; i++) {
+ imute = ixpos[i];
+ tmute = mute[isyn*nxs+imute];
+ for (j = 0; j < MAX(0,tmute-shift-smooth); j++) {
+ data[isyn*nxs*nt+i*nt+j] = 0.0;
+ }
+ for (j = MAX(0,tmute-shift-smooth),l=0; j < MAX(0,tmute-shift); j++,l++) {
+ data[isyn*nxs*nt+i*nt+j] *= costaper[smooth-l-1];
+ }
+ for (j = MIN(nt,nt-tmute+shift),l=0; j < MIN(nt,nt-tmute+shift+smooth); j++,l++) {
+ data[isyn*nxs*nt+i*nt+j] *= costaper[l];
+ }
+ for (j = MIN(nt,nt-tmute+shift+smooth); j < nt; j++) {
+ data[isyn*nxs*nt+i*nt+j] = 0.0;
+ }
+ }
+ }
+ else if (above==6){ /* implementation of <t1+epsilon : -t1+epsilon> window */
for (i = 0; i < npos; i++) {
imute = ixpos[i];
tmute = mute[isyn*nxs+imute];
@@ -69,46 +87,86 @@ void applyMute( float *data, int *mute, int smooth, int above, int Nfoc, int nxs
memset(&data[isyn*nxs*nt+i*nt],0, sizeof(float)*nt);
continue;
}
- for (j = MAX(0,-2*ts+tmute+shift),l=0; j < MAX(0,-2*ts+tmute+shift+smooth); j++,l++) {
+ for (j = MAX(0,-2*ts+tmute+shift-smooth),l=0; j < MAX(0,-2*ts+tmute+shift); j++,l++) {
data[isyn*nxs*nt+i*nt+j] *= costaper[l];
}
- for (j = MAX(0,-2*ts+tmute+shift+smooth)+1; j < MIN(nt,nt+1-tmute-shift-smooth); j++) {
+ for (j = MAX(0,-2*ts+tmute+shift); j < MIN(nt,nt-tmute+shift); j++) {
data[isyn*nxs*nt+i*nt+j] = 0.0;
}
- for (j = MIN(nt,nt-tmute-shift-smooth),l=0; j < MIN(nt,nt-tmute-shift); j++,l++) {
+ for (j = MIN(nt,nt-tmute+shift),l=0; j < MIN(nt,nt-tmute+shift+smooth); j++,l++) {
data[isyn*nxs*nt+i*nt+j] *= costaper[smooth-l-1];
}
}
}
- else if (above==-1){
+ else if (above==-6){ /* time-reversed implementation of <t1+epsilon : -t1+epsilon> window */
for (i = 0; i < npos; i++) {
imute = ixpos[i];
tmute = mute[isyn*nxs+imute];
ts = tsynW[isyn*nxs+imute];
- for (j = MAX(0,ts+tmute-shift),l=0; j < MAX(0,ts+tmute-shift+smooth); j++,l++) {
+ if (tmute >= nt/2) {
+ memset(&data[isyn*nxs*nt+i*nt],0, sizeof(float)*nt);
+ continue;
+ }
+ for (j = MAX(0,-2*ts+tmute-shift-smooth),l=0; j < MAX(0,-2*ts+tmute-shift); j++,l++) {
data[isyn*nxs*nt+i*nt+j] *= costaper[l];
}
- for (j = MAX(0,ts+tmute-shift+smooth); j < nt; j++) {
+ for (j = MAX(0,-2*ts+tmute-shift); j < MIN(nt,nt-tmute-shift); j++) {
data[isyn*nxs*nt+i*nt+j] = 0.0;
}
+ for (j = MIN(nt,nt-tmute-shift),l=0; j < MIN(nt,nt-tmute-shift+smooth); j++,l++) {
+ data[isyn*nxs*nt+i*nt+j] *= costaper[smooth-l-1];
+ }
}
}
- else if (above==4) { //Psi gate which is the inverse of the Theta gate (above=0)
+ else if (above==10) { //Psi gate which is the inverse of the above=6 gate
for (i = 0; i < npos; i++) {
imute = ixpos[i];
tmute = mute[isyn*nxs+imute];
- for (j = MAX(0,tmute-shift-smooth),l=0; j < MAX(0,tmute-shift); j++,l++) {
+ for (j = 0; j < MAX(0,tmute+shift-smooth); j++) {
+ data[isyn*nxs*nt+i*nt+j] = 0.0;
+ }
+ for (j = MAX(0,tmute+shift-smooth),l=0; j < MAX(0,tmute+shift); j++,l++) {
data[isyn*nxs*nt+i*nt+j] *= costaper[smooth-l-1];
}
- for (j = 0; j < MAX(0,tmute-shift-smooth-1); j++) {
+ for (j = MIN(nt,nt-tmute+shift),l=0; j < MIN(nt,nt-tmute+shift+smooth); j++,l++) {
+ data[isyn*nxs*nt+i*nt+j] *= costaper[l];
+ }
+ for (j = MIN(nt,nt-tmute+shift+smooth); j < nt; j++) {
data[isyn*nxs*nt+i*nt+j] = 0.0;
}
- for (j = MIN(nt,nt+1-tmute+shift+smooth); j < nt; j++) {
+ }
+ }
+ else if (above==-2){
+ for (i = 0; i < npos; i++) {
+ imute = ixpos[i];
+ tmute = mute[isyn*nxs+imute];
+ ts = tsynW[isyn*nxs+imute];
+ if (tmute >= nt/2) {
+ memset(&data[isyn*nxs*nt+i*nt],0, sizeof(float)*nt);
+ continue;
+ }
+ for (j = MAX(0,-2*ts+tmute+shift),l=0; j < MAX(0,-2*ts+tmute+shift+smooth); j++,l++) {
+ data[isyn*nxs*nt+i*nt+j] *= costaper[l];
+ }
+ for (j = MAX(0,-2*ts+tmute+shift+smooth)+1; j < MIN(nt,nt+1-tmute-shift-smooth); j++) {
data[isyn*nxs*nt+i*nt+j] = 0.0;
}
- for (j = MIN(nt,nt-tmute+shift),l=0; j < MIN(nt,nt-tmute+shift+smooth); j++,l++) {
+ for (j = MIN(nt,nt-tmute-shift-smooth),l=0; j < MIN(nt,nt-tmute-shift); j++,l++) {
+ data[isyn*nxs*nt+i*nt+j] *= costaper[smooth-l-1];
+ }
+ }
+ }
+ else if (above==-1){
+ for (i = 0; i < npos; i++) {
+ imute = ixpos[i];
+ tmute = mute[isyn*nxs+imute];
+ ts = tsynW[isyn*nxs+imute];
+ for (j = MAX(0,ts+tmute-shift),l=0; j < MAX(0,ts+tmute-shift+smooth); j++,l++) {
data[isyn*nxs*nt+i*nt+j] *= costaper[l];
}
+ for (j = MAX(0,ts+tmute-shift+smooth); j < nt; j++) {
+ data[isyn*nxs*nt+i*nt+j] = 0.0;
+ }
}
}
else if (above==2){//Separates the direct part of the wavefield from the coda
diff --git a/marchenko/marchenko.c b/marchenko/marchenko.c
index c75078a5a6c638c7431a4ba41f98eede7bd9d908..0d0f1ae93bf215e21356eac23868a8680c59a2e5 100644
--- a/marchenko/marchenko.c
+++ b/marchenko/marchenko.c
@@ -492,7 +492,7 @@ int main (int argc, char **argv)
applyMute_tshift(Ni, muteW, smooth, above, Nfoc, nxs, nts, ixpos, npos, shift, iter, tsynW);
}
else {
- applyMute(Ni, muteW, smooth, above, Nfoc, nxs, nts, ixpos, npos, shift, tsynW);
+ applyMute(Ni, muteW, smooth, -above, Nfoc, nxs, nts, ixpos, npos, shift, tsynW);
}
if (iter % 2 == 0) { /* even iterations update: => f_1^-(t) */