diff --git a/fdelmodc/applySource.c b/fdelmodc/applySource.c
index f7b3a85014afdc5913b2b76051b0ad29ea79a4c5..d29016e9b0e55cbe876f66c4fd246ad439614a49 100644
--- a/fdelmodc/applySource.c
+++ b/fdelmodc/applySource.c
@@ -27,14 +27,9 @@ int applySource(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int i
int isrc, ix, iz, n1;
int id1, id2;
float src_ampl, time, scl, dt, sdx;
- float Mxx, Mzz, Mxz, rake, strike;
+ float Mxx, Mzz, Mxz;
static int first=1;
- if (!getparfloat("strike",&strike)) strike=0.5*M_PI;
- if (!getparfloat("rake",&rake)) rake=0.5*M_PI;
- if (strike!=0.5*M_PI) strike = M_PI*(strike/180.0);
- if (rake!=0.5*M_PI) rake = M_PI*(rake/180.0);
-
if (src.type==6) {
ibndz = mod.ioXz;
ibndx = mod.ioXx;
@@ -342,14 +337,10 @@ int applySource(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int i
vz[ix*n1+iz] += src_ampl*sdx;
}
}
- else if(src.type == 9) {
- Mxx = -1.0*(sin(src.dip)*cos(rake)*sin(2.0*strike)+sin(src.dip*2.0)*sin(rake)*sin(strike)*sin(strike));
- Mxz = -1.0*(cos(src.dip)*cos(rake)*cos(strike)+cos(src.dip*2.0)*sin(rake)*sin(strike));
- Mzz = sin(src.dip*2.0)*sin(rake);
-
- txx[ix*n1+iz] -= Mxx*src_ampl;
- tzz[ix*n1+iz] -= Mzz*src_ampl;
- txz[ix*n1+iz] -= Mxz*src_ampl;
+ else if(src.type == 9 || src.type == 11) {
+ txx[ix*n1+iz] -= src.Mxx*src_ampl;
+ tzz[ix*n1+iz] -= src.Mzz*src_ampl;
+ txz[ix*n1+iz] -= src.Mxz*src_ampl;
} /* src.type */
} /* ischeme */
}
diff --git a/fdelmodc/demo/fdelmodc_moment_tensor.scr b/fdelmodc/demo/fdelmodc_moment_tensor.scr
new file mode 100755
index 0000000000000000000000000000000000000000..6fa9912a7511b466c8e839bc98625cb8cefaf0ca
--- /dev/null
+++ b/fdelmodc/demo/fdelmodc_moment_tensor.scr
@@ -0,0 +1,59 @@
+#!/bin/bash
+#PBS -N fdelmodc
+#PBS -k eo
+#PBS -j eo
+#
+# Models plane wave at depth to receivers at the surface, including snapshots
+export PATH=../../utils:$PATH:
+
+makewave file_out=wavelet.su dt=0.001 nt=1024 fp=13 shift=1 w=g2 verbose=1
+
+makemod file_base=model.su \
+ cp0=1500 ro0=1000 cs0=0 sizex=2100 sizez=2100 \
+ dx=3 dz=3 orig=0,0 \
+ verbose=4
+
+export filecp=model_cp.su
+export filecs=model_cs.su
+export filero=model_ro.su
+
+export OMP_NUM_THREADS=1
+
+../fdelmodc \
+ file_cp=$filecp file_den=$filero \
+ ischeme=5 \
+ src_type=11 Mxx=0.0 Mzz=0.0 Mxz=-1.0 tmod=2.0 \
+ file_src=wavelet.su verbose=2 \
+ file_rcv=rec4S.su \
+ file_snap=snap4S0.su \
+ xrcv1=0 xrcv2=2100 dxrcv=15 \
+ zrcv1=400 zrcv2=400 \
+ rec_type_vx=1 rec_int_vx=1 \
+ dtrcv=0.004 \
+ xsrc=1050 zsrc=1050 \
+ ntaper=120 sna_type_pp=1 sna_type_ss=1 sna_type_p=1 sna_type_txx=1 sna_type_tzz=1 \
+ left=4 right=4 bottom=4 top=4 \
+ tsnap1=0.1 tsnap2=3.0 dtsnap=0.01
+
+
+exit;
+../fdelmodc \
+ file_cp=$filecp file_cs=$filecs file_den=$filero \
+ ischeme=3 \
+ file_src=wavelet.su verbose=4 \
+ file_rcv=recP.su \
+ file_snap=snapP.su \
+ xrcv1=0 xrcv2=2100 dxrcv=15 \
+ zrcv1=400 zrcv2=400 \
+ rec_type_vx=1 rec_type_pp=1 rec_type_ss=1 rec_int_vx=1 \
+ dtrcv=0.004 \
+ xsrc=1000 zsrc=1700 \
+ src_type=8 src_orient=1 tmod=2.0 \
+ ntaper=120 \
+ left=4 right=4 bottom=4 top=4 \
+ tsnap1=0.1 tsnap2=3.0 dtsnap=0.1 \
+ sna_type_ss=1 sna_type_pp=1
+
+
+exit;
+
diff --git a/fdelmodc/fdelmodc.h b/fdelmodc/fdelmodc.h
index ed07c68da00636d4e080b9b7baf0393f8218dc62..f68d00cfcb3fe8d3b0aff1ba4d2abb0c2969ae32 100644
--- a/fdelmodc/fdelmodc.h
+++ b/fdelmodc/fdelmodc.h
@@ -132,6 +132,9 @@ typedef struct _sourcePar { /* Source Array Parameters */
int orient;
int *z;
int *x;
+ float Mxx;
+ float Mzz;
+ float Mxz;
int single;
int plane;
int circle;
diff --git a/fdelmodc/getParameters.c b/fdelmodc/getParameters.c
index a316a1f53635357c73247199b83486bd47b12810..598267d5a2c6ccb1e90f86b9d2586ec6abfb7f13 100644
--- a/fdelmodc/getParameters.c
+++ b/fdelmodc/getParameters.c
@@ -52,7 +52,7 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
float xsnap1, xsnap2, zsnap1, zsnap2, xmax, zmax;
float xsrc1, xsrc2, zsrc1, zsrc2, tsrc1, tsrc2, tlength, tactive;
float src_angle, src_velo, p, grad2rad, rdelay, scaledt;
- float *xsrca, *zsrca, rrcv;
+ float *xsrca, *zsrca, rrcv, strike, rake, dip;
float rsrc, oxsrc, ozsrc, dphisrc, ncsrc;
size_t nsamp;
int i, j;
@@ -552,11 +552,21 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
if (!getparint("nshot",&shot->n)) shot->n=1;
if (!getparfloat("dxshot",&dxshot)) dxshot=dx;
if (!getparfloat("dzshot",&dzshot)) dzshot=0.0;
+ if (!getparfloat("Mxx",&src->Mxx)) src->Mxx=1.0;
+ if (!getparfloat("Mzz",&src->Mzz)) src->Mzz=1.0;
+ if (!getparfloat("Mxz",&src->Mxz)) src->Mxz=1.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;
- src->dip = M_PI*(src->dip/180.0);
+ if (!getparfloat("strike",&strike)) strike=90.0;
+ if (!getparfloat("rake",&rake)) rake=90.0;
+ strike = M_PI*(strike/180.0);
+ rake = M_PI*(rake/180.0);
+ dip = M_PI*(dip/180.0);
+
+ if (src->type==9) {
+ src->Mxx = -1.0*(sin(dip)*cos(rake)*sin(2.0*strike)+sin(dip*2.0)*sin(rake)*sin(strike)*sin(strike));
+ src->Mxz = -1.0*(cos(dip)*cos(rake)*cos(strike)+cos(dip*2.0)*sin(rake)*sin(strike));
+ src->Mzz = sin(dip*2.0)*sin(rake);
+ }
if (shot->n>1) {
idxshot=MAX(0,NINT(dxshot/dx));
@@ -960,8 +970,11 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
case 8 : fprintf(stderr,"P-potential"); break;
case 9 : fprintf(stderr,"double-couple"); break;
case 10 : fprintf(stderr,"Fz on P grid with +/-"); break;
+ case 11 : fprintf(stderr,"moment tensor"); break;
}
fprintf(stderr,"\n");
+ if (src->type==9) vmess("strike %.2f rake %.2f dip %.2f",180.0*strike/M_PI,180.0*rake/M_PI,180.0*dip/M_PI);
+ if (src->type==9 || src->type==11) vmess("Mxx %.2f Mzz %.2f Mxz %.2f",src->Mxx,src->Mzz,src->Mxz);
if (wav->random) vmess("Wavelet has a random signature with fmax=%.2f", wav->fmax);
if (src->n>1) {
vmess("*******************************************");
diff --git a/utils/HomG.c b/utils/HomG.c
index 5abf0f62b8e86c4f6d5e75b23cfb98e5a276cbe4..9857969bafa5bf0f434a8538064e03577e7dda8f 100755
--- a/utils/HomG.c
+++ b/utils/HomG.c
@@ -119,7 +119,7 @@ int main (int argc, char **argv)
float dt, dy, dx, t0, y0, x0, xmin, xmax1, sclsxgx, dxrcv, dyrcv, dzrcv;
float *conv, *conv2, *tmp1, *tmp2, cp, shift;
long nshots, ntvs, nyvs, nxvs, ntraces, ret, ix, iy, it, is, ir, ig, file_det, verbose;
- long ntr, nxr, nyr, nsr, i, l, j, k, nxvr, nyvr, nzvr, count, num;
+ long ntr, nxr, nyr, nsr, i, l, j, k, nxvr, nyvr, nzvr, count, num, isn;
float dtr, dxr, dyr, ftr, fxr, fyr, sclr, scl;
long pos1, npos, zmax, numb, dnumb, scheme, ntmax, ntshift, shift_num, zfps, zfpr, size;
long ixr, iyr, zsrc, zrcv, *xvr, *yvr, *zvr;
@@ -172,6 +172,7 @@ int main (int argc, char **argv)
num = numb;
count = 0;
+ if (num==0) count=1;
while (num != 0) {
count++;
num /= 10;
@@ -607,35 +608,38 @@ int main (int argc, char **argv)
fp_out = fopen(fout, "w+");
hdr_out = (segy *)calloc(nxvr*nyvr,sizeof(segy));
- for (ir = 0; ir < ntr; ir++) {
- for (is = 0; is < nyvr; is++) {
- for (ix = 0; ix < nxvr; ix++) {
- hdr_out[is*nxvr+ix].fldr = ir-ntr/2;
- hdr_out[is*nxvr+ix].tracl = ir*nyvr*nxvr+is*nxvr+ix+1;
- hdr_out[is*nxvr+ix].tracf = is*nxvr+ix+1;
- hdr_out[is*nxvr+ix].scalco = -1000;
- hdr_out[is*nxvr+ix].scalel = -1000;
- hdr_out[is*nxvr+ix].sdepth = hdr_shot[0].sdepth;
- hdr_out[is*nxvr+ix].selev = -hdr_shot[0].sdepth;
- hdr_out[is*nxvr+ix].trid = 1;
- hdr_out[is*nxvr+ix].ns = nzvr;
- hdr_out[is*nxvr+ix].trwf = nxvr*nyvr;
- hdr_out[is*nxvr+ix].ntr = (ir+1)*hdr_out[is*nxvr+ix].trwf;
- hdr_out[is*nxvr+ix].f1 = ((float)(zvr[0]/1000.0));
- hdr_out[is*nxvr+ix].f2 = ((float)(xvr[0]/1000.0));
- hdr_out[is*nxvr+ix].dt = dt*(1E6);
- hdr_out[is*nxvr+ix].d1 = dzrcv;
- hdr_out[is*nxvr+ix].d2 = dxrcv;
- hdr_out[is*nxvr+ix].sx = hdr_shot[0].sx;
- hdr_out[is*nxvr+ix].sy = hdr_shot[0].sy;
- hdr_out[is*nxvr+ix].gx = 1000.0*((float)(xvr[0]/1000.0)+dxrcv*ix);
- hdr_out[is*nxvr+ix].gy = 1000.0*((float)(yvr[0]/1000.0)+dyrcv*is);
- hdr_out[is*nxvr+ix].offset = (hdr_out[is*nxvr+ix].gx - hdr_out[is*nxvr+ix].sx)/1000.0;
+ for (isn = 0; isn < nshots; isn++) {
+ for (ir = 0; ir < ntr; ir++) {
+ for (is = 0; is < nyvr; is++) {
+ for (ix = 0; ix < nxvr; ix++) {
+ hdr_out[is*nxvr+ix].fldr = ir-ntr/2;
+ hdr_out[is*nxvr+ix].tracl = ir*nyvr*nxvr+is*nxvr+ix+1;
+ hdr_out[is*nxvr+ix].tracf = is*nxvr+ix+1;
+ hdr_out[is*nxvr+ix].tracr = isn;
+ hdr_out[is*nxvr+ix].scalco = -1000;
+ hdr_out[is*nxvr+ix].scalel = -1000;
+ hdr_out[is*nxvr+ix].sdepth = hdr_shot[isn*nxvs*nyvs].sdepth;
+ hdr_out[is*nxvr+ix].selev = hdr_shot[isn*nxvs*nyvs].selev;
+ hdr_out[is*nxvr+ix].trid = 1;
+ hdr_out[is*nxvr+ix].ns = nzvr;
+ hdr_out[is*nxvr+ix].trwf = nxvr*nyvr;
+ hdr_out[is*nxvr+ix].ntr = (ir+1)*hdr_out[is*nxvr+ix].trwf;
+ hdr_out[is*nxvr+ix].f1 = ((float)(zvr[0]/1000.0));
+ hdr_out[is*nxvr+ix].f2 = ((float)(xvr[0]/1000.0));
+ hdr_out[is*nxvr+ix].dt = dt*(1E6);
+ hdr_out[is*nxvr+ix].d1 = dzrcv;
+ hdr_out[is*nxvr+ix].d2 = dxrcv;
+ hdr_out[is*nxvr+ix].sx = hdr_shot[isn*nxvs*nyvs].sx;
+ hdr_out[is*nxvr+ix].sy = hdr_shot[isn*nxvs*nyvs].sy;
+ hdr_out[is*nxvr+ix].gx = 1000.0*((float)(xvr[0]/1000.0)+dxrcv*ix);
+ hdr_out[is*nxvr+ix].gy = 1000.0*((float)(yvr[0]/1000.0)+dyrcv*is);
+ hdr_out[is*nxvr+ix].offset = (hdr_out[is*nxvr+ix].gx - hdr_out[is*nxvr+ix].sx)/1000.0;
+ }
}
- }
- ret = writeData3D(fp_out, &Ghom[ir*nyvr*nxvr*nzvr], hdr_out, nzvr, nxvr*nyvr);
- if (ret < 0 ) verr("error on writing output file.");
- }
+ ret = writeData3D(fp_out, &Ghom[isn*ntr*nyvr*nxvr*nzvr+ir*nyvr*nxvr*nzvr], hdr_out, nzvr, nxvr*nyvr);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+ }
fclose(fp_out);
free(hdr_shot);
diff --git a/utils/Makefile b/utils/Makefile
index cac463dd9efd0438c87e81d69a568f0fa499869b..115653916d18f9d1294b731e133ae2f50b80511f 100644
--- a/utils/Makefile
+++ b/utils/Makefile
@@ -7,7 +7,7 @@ LIBS += -L$L -lgenfft -lzfp
#OPTC += -openmp
#OPTC += -g -O0
-ALL: makemod makewave extendModel fconv correigen green green3D basop syn2d mat2su ftr1d mutesnap padmodel truncate combine combine_induced reshape_su HomG snap2shot makeR1D pwshift
+ALL: makemod makewave extendModel fconv correigen green green3D basop syn2d mat2su ftr1d mutesnap padmodel truncate combine combine_induced reshape_su HomG snap2shot makeR1D pwshift convhomg
SRCM = \
makemod.c \
@@ -235,6 +235,16 @@ SRCPW = pwshift.c \
docpkge.c \
readSnapData3D.c
+SRCCH = convhomg.c \
+ getFileInfo3D.c \
+ writeData3D.c \
+ wallclock_time.c \
+ getpars.c \
+ verbosepkg.c \
+ atopkge.c \
+ docpkge.c \
+ readSnapData3D.c
+
OBJM = $(SRCM:%.c=%.o)
makemod: $(OBJM)
@@ -340,7 +350,12 @@ OBJPW = $(SRCPW:%.c=%.o)
pwshift: $(OBJPW)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o pwshift $(OBJPW) $(LIBS)
-install: makemod makewave extendModel fconv correigen green green3D basop syn2d mat2su ftr1d mutesnap padmodel truncate combine combine_induced reshape_su HomG snap2shot makeR1D pwshift
+OBJCH = $(SRCCH:%.c=%.o)
+
+convhomg: $(OBJCH)
+ $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o convhomg $(OBJCH) $(LIBS)
+
+install: makemod makewave extendModel fconv correigen green green3D basop syn2d mat2su ftr1d mutesnap padmodel truncate combine combine_induced reshape_su HomG snap2shot makeR1D pwshift convhomg
cp makemod $B
cp makewave $B
cp extendModel $B
@@ -362,9 +377,10 @@ install: makemod makewave extendModel fconv correigen green green3D basop syn2d
cp snap2shot $B
cp makeR1D $B
cp pwshift $B
+ cp convhomg $B
clean:
- rm -f core $(OBJM) makemod $(OBJW) makewave $(OBJE) extendModel $(OBJF) fconv $(OBJG) $(OBJC) correigen green $(OBJG3) green3D $(OBJB) basop $(OBJJ) syn2d $(OBJS) mat2su $(OBJA) ftr1d $(OBJT) mutesnap $(OBJMS) truncate $(OBJTR) padmodel $(OBJPM) combine $(OBJCO) makeR1D $(OBJMR) reshape_su $(OBJRS) combine_induced $(OBJCI) HomG $(OBJHG) snap2shot $(OBJSS) pwshift $(OBJPW)
+ rm -f core $(OBJM) makemod $(OBJW) makewave $(OBJE) extendModel $(OBJF) fconv $(OBJG) $(OBJC) correigen green $(OBJG3) green3D $(OBJB) basop $(OBJJ) syn2d $(OBJS) mat2su $(OBJA) ftr1d $(OBJT) mutesnap $(OBJMS) truncate $(OBJTR) padmodel $(OBJPM) combine $(OBJCO) makeR1D $(OBJMR) reshape_su $(OBJRS) combine_induced $(OBJCI) HomG $(OBJHG) snap2shot $(OBJSS) pwshift $(OBJPW) convhomg $(OBJCH)
realclean: clean
- rm -f $B/makemod $B/makewave $B/extendModel $B/fconv $B/correigen $B/green $B/green3D $B/basop $B/syn2d $B/mat2su $B/ftr1d $B/mutesnap $B/padmodel $B/truncate $B/combine $B/combine_induced $B/reshape_su $B/HomG $B/snap2shot $B/makeR1D $B/pwshift
+ rm -f $B/makemod $B/makewave $B/extendModel $B/fconv $B/correigen $B/green $B/green3D $B/basop $B/syn2d $B/mat2su $B/ftr1d $B/mutesnap $B/padmodel $B/truncate $B/combine $B/combine_induced $B/reshape_su $B/HomG $B/snap2shot $B/makeR1D $B/pwshift $B/convhomg
diff --git a/utils/mutesnap.c b/utils/mutesnap.c
index 037aae7dbed6892b003d5e2810ed87e979b19834..7afa393e19ddef9bd7f8c2354e8b6d742fff77b7 100644
--- a/utils/mutesnap.c
+++ b/utils/mutesnap.c
@@ -50,6 +50,8 @@ char *sdoc[] = {
" mode=0 ................... Determine first arrival by maximum (mode=0), first event above tol (mode=1) or by raytime (mode=2)",
" tol=1 .................... Tolerance for the determination of first arrival if mode=1",
" fray ..................... File containing the raytimes of the first arrivals",
+" opt=0 .................... Mute the file in the center (=0) or at the beginning (=1)",
+" verbose=1 ................ Give detailed information about the process (=1) or remain silent (=0)",
NULL};
int main (int argc, char **argv)
@@ -60,7 +62,7 @@ int main (int argc, char **argv)
float dxs, dys, dzs, scls, fzs, fxs, fys;
float dxh, dyh, dzh, sclh, fzh, fxh, fyh;
float dxrcv, dyrcv, dzrcv, dxpos, offset;
- long nts, nxs, nys, nzs, ntrs, nth, nxh, nyh, nzh, ntrh;
+ long nts, nxs, nys, nzs, ntrs, nth, nxh, nyh, nzh, ntrh, opt;
long nxyz, nxy, ret, ix, iy, iz, it, ht, indrcv, shift, rmt, mode, smooth, verbose;
segy *hdr_hom, *hdr_snap, *hdrs_mute;
@@ -78,6 +80,7 @@ int main (int argc, char **argv)
if (!getparlong("smooth", &smooth)) smooth = 5;
if (!getparlong("mode", &mode)) mode = 0;
if (!getparlong("verbose", &verbose)) verbose = 1;
+ if (!getparlong("opt", &opt)) opt = 0;
if (!getparfloat("tol", &tol)) tol = 5;
if (fhom == NULL) verr("Incorrect G_hom input");
if (mode != 2) {
@@ -163,38 +166,72 @@ int main (int argc, char **argv)
/*----------------------------------------------------------------------------*
* Apply the muting to the data
*----------------------------------------------------------------------------*/
- for (iz = 0; iz < nzh; iz++) {
- for (iy = 0; iy < nyh; iy++) {
- for (ix = 0; ix < nxh; ix++) {
- if (mode != 2) {
- for (it = 0; it < nth; it++) {
- rtrace[it] = snapdata[it*nxyz+iy*nxh*nzh+ix*nzh+iz];
+ if (opt==0) {
+ if (verbose) vmess("muting around the center");
+ for (iz = 0; iz < nzh; iz++) {
+ for (iy = 0; iy < nyh; iy++) {
+ for (ix = 0; ix < nxh; ix++) {
+ if (mode != 2) {
+ for (it = 0; it < nth; it++) {
+ rtrace[it] = snapdata[it*nxyz+iy*nxh*nzh+ix*nzh+iz];
+ }
+ }
+ if (mode == 0) {
+ indrcv = ht - topdet(rtrace,nth);
+ }
+ else if (mode == 1) {
+ indrcv = ht - farrdet(rtrace,nth,tol);
+ }
+ else if (mode == 2) {
+ indrcv = (long)roundf(timeval[iz*nxh*nyh+iy*nxh+ix]/dt);
+ }
+ rmt = MAX(MIN(nth-indrcv,indrcv)-shift-smooth,0);
+ for (it = ht-rmt+1; it < ht+1; it++) {
+ if (it-(ht-rmt+1) < smooth) {
+ homdata[it*nyh*nxh*nzh+iy*nxh*nzh+ix*nzh+iz] *= costaper[it-(ht-rmt+1)];
+ homdata[(nth-it)*nyh*nxh*nzh+iy*nxh*nzh+ix*nzh+iz] *= costaper[it-(ht-rmt+1)];
+ }
+ else{
+ homdata[it*nyh*nxh*nzh+iy*nxh*nzh+ix*nzh+iz] = 0.0;
+ homdata[(nth-it)*nyh*nxh*nzh+iy*nxh*nzh+ix*nzh+iz] = 0.0;
+ }
}
}
- if (mode == 0) {
- indrcv = ht - topdet(rtrace,nth);
- }
- else if (mode == 1) {
- indrcv = ht - farrdet(rtrace,nth,tol);
- }
- else if (mode == 2) {
- indrcv = (long)roundf(timeval[iz*nxh*nyh+iy*nxh+ix]/dt);
- }
- rmt = MAX(MIN(nth-indrcv,indrcv)-shift-smooth,0);
- for (it = ht-rmt+1; it < ht+1; it++) {
- if (it-(ht-rmt+1) < smooth) {
- homdata[it*nyh*nxh*nzh+iy*nxh*nzh+ix*nzh+iz] *= costaper[it-(ht-rmt+1)];
- homdata[(nth-it)*nyh*nxh*nzh+iy*nxh*nzh+ix*nzh+iz] *= costaper[it-(ht-rmt+1)];
+ }
+ if (verbose) vmess("Muting Homogeneous Green's function at depth %li from %li depths",iz+1,nzh);
+ }
+ }
+ else if (opt==1) {
+ if (verbose) vmess("muting at the start");
+ for (iz = 0; iz < nzh; iz++) {
+ for (iy = 0; iy < nyh; iy++) {
+ for (ix = 0; ix < nxh; ix++) {
+ if (mode != 2) {
+ for (it = 0; it < nth; it++) {
+ rtrace[it] = snapdata[it*nxyz+iy*nxh*nzh+ix*nzh+iz];
+ }
+ }
+ if (mode == 0) {
+ indrcv = topdet(rtrace,nth);
+ }
+ else if (mode == 1) {
+ indrcv = farrdet(rtrace,nth,tol);
}
- else{
+ else if (mode == 2) {
+ indrcv = (long)roundf(timeval[iz*nxh*nyh+iy*nxh+ix]/dt);
+ }
+ for (it = MAX(indrcv-shift-smooth,0); it < MAX(indrcv-shift,0); it++) {
+ homdata[it*nyh*nxh*nzh+iy*nxh*nzh+ix*nzh+iz] *= costaper[it-(indrcv-shift)-1];
+ }
+ for (it = 0; it < MAX(indrcv-shift-smooth,0); it++) {
homdata[it*nyh*nxh*nzh+iy*nxh*nzh+ix*nzh+iz] = 0.0;
- homdata[(nth-it)*nyh*nxh*nzh+iy*nxh*nzh+ix*nzh+iz] = 0.0;
}
}
}
+ if (verbose) vmess("Muting Homogeneous Green's function at depth %li from %li depths",iz+1,nzh);
}
- if (verbose) vmess("Muting Homogeneous Green's function at depth %li from %li depths",iz+1,nzh);
}
+
free(rtrace);
if (smooth) free(costaper);
if (mode == 2) {
diff --git a/utils/padmodel.c b/utils/padmodel.c
index a4723dd5e566ef6793df8c5f31ad87a4082a869d..8a2de11178fe00466d245babdab285f91aba36a9 100644
--- a/utils/padmodel.c
+++ b/utils/padmodel.c
@@ -44,9 +44,9 @@ char *sdoc[] = {
" Optional parameters: ",
" ",
" file_out=out.su .......... Filename of the output",
-" nx=0 ..................... Number of samples to pad on both the left and right",
-" ny=0 ..................... Number of samples to pad on both the front and back",
-" nz=0 ..................... Number of samples to pad on only the bottom",
+" nxpad=0 .................. Number of samples to pad on both the left and right",
+" nypad=0 .................. Number of samples to pad on both the front and back",
+" nzpad=0 .................. Number of samples to pad on both the bottom and top",
" pad=0 .................... Pad option, (=0) pad with the value at the edge, (=1) pad with a constant value",
" value=0.0 ................ Padding value if pad=1",
" verbose=1 ................ Give detailed information of process",
@@ -88,11 +88,11 @@ int main (int argc, char **argv)
nxout = nx + 2*nxpad;
nyout = ny + 2*nypad;
- nzout = nz + nzpad;
+ nzout = nz + 2*nzpad;
x0out = x0 - ((float)nxpad)*dx;
- y0out = y0 - ((float)nypad)*dx;
- z0out = z0;
+ y0out = y0 - ((float)nypad)*dy;
+ z0out = z0 - ((float)nzpad)*dz;
if (verbose) {
vmess("******************** INPUT DATA ********************");
@@ -103,7 +103,7 @@ int main (int argc, char **argv)
vmess("******************** PADDING ********************");
if (!pad) vmess("Model is padded using the edge values");
else vmess("model is padded using constant value %.3f",value);
- vmess("Number of padding samples: x: %li, y: %li, z: %li",2*nxpad,2*nypad,nzpad);
+ vmess("Number of padding samples: x: %li, y: %li, z: %li",2*nxpad,2*nypad,2*nzpad);
vmess("******************** OUTPUT DATA ********************");
vmess("Number of samples: %li, x: %li, y: %li, z: %li",nxout*nyout*nzout,nxout,nyout,nzout);
vmess("Sampling distance for x: %.3f, y: %.3f, z: %.3f",dx,dy,dz);
@@ -130,7 +130,7 @@ int main (int argc, char **argv)
for (iy = 0; iy < ny; iy++) {
for (ix = 0; ix < nx; ix++) {
for (iz = 0; iz < nz; iz++) {
- outdata[(iy+nypad)*nxout*nzout+(ix+nxpad)*nzout+iz] = indata[iy*nx*nz+ix*nz+iz];
+ outdata[(iy+nypad)*nxout*nzout+(ix+nxpad)*nzout+iz+nzpad] = indata[iy*nx*nz+ix*nz+iz];
}
}
}
@@ -177,7 +177,7 @@ int main (int argc, char **argv)
if (!pad) {
if (verbose) vmess("Padding with edge value");
- for (iz = 0; iz < nz; iz++) {
+ for (iz = nzpad; iz < nz+nzpad; iz++) {
for (iy = nypad; iy < ny+nypad; iy++) {
for (ix = nxpad-1; ix >= 0; ix--) {
outdata[iy*nxout*nzout+ix*nzout+iz] = outdata[iy*nxout*nzout+nxpad*nzout+iz];
@@ -195,6 +195,13 @@ int main (int argc, char **argv)
}
}
}
+ for (iz = 0; iz < nzpad; iz++) {
+ for (iy = 0; iy < nyout; iy++) {
+ for (ix = 0; ix < nxout; ix++) {
+ outdata[iy*nxout*nzout+ix*nzout+iz] = outdata[iy*nxout*nzout+ix*nzout+nzpad];
+ }
+ }
+ }
for (iz = nz; iz < nzout; iz++) {
for (iy = 0; iy < nyout; iy++) {
for (ix = 0; ix < nxout; ix++) {