Skip to content
Snippets Groups Projects
Commit 5961b4c7 authored by JanThorbecke's avatar JanThorbecke
Browse files

bugs in planewave calculations and adding a few extra options

parent ebca6f18
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
marchenko/demo/oneD/README_PRIMARIES
+ 33
− 32
View file @ 5961b4c7
......@@ -42,39 +42,40 @@ The full R matrix is build up from the the shot record computed with model.scr
This R is the only input of the Marchenko Primaries algorithm.
==> run './iterations.scr ./iterations.scr Figure3467' to compute the intermediate results of the first iterations of the Marchenko Primaries algorithm.
This will take 10 seconds. The generated files are:
- M0_276000.su
- iter_2760##.su (not used)
- Mi_2760##.su
- f1min_2760##.su
- pred_rr_276.su (not used)
where ## ranges from 01 to 33
To generate the postscript files for Figure 3:
==> run './epsPrimaries.scr Figure3'
This will produce the following files:
shotx0_rp.eps => Figure 2c == Figure 3a
M0_276000_flip.eps => Figure 3b
fconvN0fulltime.eps => Figure 3c
fconvN0flip.eps => Figure 3d
Mi_276001.eps => Figure 3e
To generate the postscript files for Figure 4:
==> run './epsPrimaries.scr Figure4'
To generate the postscript files for Figure 6:
==> run './epsPrimaries.scr Figure6'
To generate the postscript files for Figure 7:
==> run './epsPrimaries.scr Figure7'
==> run marchenkoIter.scr to compute the first 4 iteration of the Marchenko algorithm. This will take 1-2 minutes. The generated files are:
- p0plus.su
- pgreen_001.su
- f1plus_001.su
- f1min_001.su
- Gplus_001.su
- Gmin_001.su
- pgreen_002.su
- f1plus_002.su
- f1min_002.su
- Gplus_002.su
- Gmin_002.su
- pgreen_003.su
- f1plus_003.su
- f1min_003.su
- Gplus_003.su
- Gmin_003.su
- pgreen_004.su
- f1plus_004.su
- f1min_004.su
- Gplus_004.su
- Gmin_004.su
To generate all postscript files for Figure 4, 5 and 6
==> run './epsPrimaries.scr Figure3' to generate the postscript files of Figure 3
==> run epsMarchenkoIter.scr
shotx0_rp.eps => Figure 4 R == Figure 3c
p0plus.eps => Figure 4 G_d
iter_001.eps => Figure 4 N_0
shotx0_rp.eps => Figure 5 R == Figure 3c
f1min_001.eps => Figure 5 f^-_1,0
......
marchenko/demo/oneD/epsPrimaries.scr
+ 33
− 15
View file @ 5961b4c7
......@@ -34,6 +34,11 @@ echo "clipiter="$clipiter "clipref="$clipref "clipf1="$clipf1
clipiter=$clipref
clipf1=$clipref
# Initialisation and First iteration
if [[ "$1" == "Figure3" ]];
then
file=M0_276000.su
file_base=${file%.su}
supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
n1tic=2 d2=5 x1beg=$nstart x1end=$ns d1=1 f1=0 f1num=700 d1num=100 \
label1="time sample number" label2="lateral distance" \
......@@ -45,7 +50,6 @@ supsimage hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
label1="time sample number" label2="lateral distance" \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > ${file_base}_flip.eps
# First iteration
#convolve M0 with middle shot record of R
select=451
#original shot record from Reflection matrix
......@@ -58,9 +62,16 @@ fconv file_in1=shotsx0.su file_in2=N0.su file_out=fconvN0.su verbose=1 fmax=90
file=fconvN0.su
file_base=${file%.su}
supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
n1tic=2 d2=5 x1beg=624 x1end=1024 d1=1 f1=0 f1num=700 d1num=100 \
n1tic=2 d2=5 x1beg=0 d1=1 f1=0 f1num=0 d1num=200 \
label1="time sample number" label2="lateral distance" \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > ${file_base}fulltime.eps
file=fconvN0.su
file_base=${file%.su}
supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
n1tic=2 d2=5 x1beg=$nstart x1end=$ns d1=1 f1=0 f1num=700 d1num=100 \
label1="time sample number" label2="lateral distance" \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > ${file_base}zoom.eps
suflip < fconvN0.su flip=3 | sugain scale=-1 > fconvN0flip.su
file=fconvN0flip.su
......@@ -70,14 +81,21 @@ supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
label1="time sample number" label2="lateral distance" \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
file=fconvN0.su
iter=1
piter=$(printf %03d $iter)
file=Mi_${istart}${piter}.su
file_base=${file%.su}
supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
n1tic=2 d2=5 x1beg=0 d1=1 f1=0 f1num=0 d1num=200 \
n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 d1num=100 \
label1="time sample number" label2="lateral distance" \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > ${file_base}fulltime.eps
f2=-2250 f2num=-2000 d2num=1000 clip=$clipiter > $file_base.eps
fi
# second iteration
if [[ "$1" == "Figure4" ]];
then
suwind itmax=1023 < Mi_276001.su > N0.su
#compute R*N0
fconv file_in1=shotsx0.su file_in2=N0.su file_out=fconvN1.su verbose=1 fmax=90
......@@ -96,29 +114,27 @@ supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
n1tic=2 d2=5 x1beg=0 d1=1 f1=0 f1num=0 d1num=200 \
label1="time sample number" label2="lateral distance" \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipm1 > ${file_base}fulltime.eps
#RNi
file=iter_276000.su
file_base=${file%.su}
#clipref=$clipiter
supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
n1tic=2 d2=5 x1beg=0 d1=1 d1num=100 \
label1="time sample number" label2="lateral distance" \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
fi
#Iterations
if [[ "$1" == "Figure5" ]];
then
for (( iter=0; iter<=21; iter+=2 ))
do
piter=$(printf %03d $iter)
echo $piter
file=Mi_${istart}${piter}.su
file_base=${file%.su}
#possibly add suflip flip=3 to flip the time-axis
supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
n1tic=2 d2=5 x1beg=$nstart x1end=$ns d1=1 d1num=100 \
label1="time sample number" label2="lateral distance" \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipiter > $file_base.eps
done
fi
if [[ "$1" == "Figure5" ]];
then
for (( iter=1; iter<=21; iter+=2 ))
do
piter=$(printf %03d $iter)
......@@ -130,7 +146,9 @@ supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
label1="time sample number" label2="lateral distance" \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipiter > $file_base.eps
done
fi
exit;
#iterations
for (( iter=0; iter<=31; iter+=2 ))
do
......
marchenko/demo/oneD/iterations.scr
+ 17
− 5
View file @ 5961b4c7
......@@ -9,23 +9,35 @@
select=451
#for istart in 276
for istart in 246 256 266 276 286 296 306 316
if [[ "$1" == "Figure3467" ]];
then
for istart in 276
do
(( iend = istart + 1 ))
../../marchenko_primaries file_shot=shotsdx5_rp.su ishot=$select file_src=wave.su \
nshots=901 verbose=2 istart=$istart iend=$iend fmax=90 \
pad=512 niter=33 smooth=8 niterskip=600 shift=16 file_rr=pred_rr_${istart}.su T=0 file_iter=iter.su
done
fi
#suximage < f1min_${istart}043.su clip=1 title="${istart}" &
if [[ "$1" == "Figure8" ]];
then
for istart in 246 256 266 276 286 296 306 316
do
(( iend = istart + 1 ))
../../marchenko_primaries file_shot=shotsdx5_rp.su ishot=$select file_src=wave.su \
nshots=901 verbose=2 istart=$istart iend=$iend fmax=90 \
pad=512 niter=33 smooth=8 niterskip=600 shift=16 file_rr=pred_rr_${istart}.su T=0 file_iter=iter.su
done
fi
exit;
#T-MME
if [[ "$1" == "Figure10" ]];
then
istart=276
(( iend = istart + 1 ))
../../marchenko_primaries file_shot=shotsdx5_rp.su ishot=$select file_src=wave.su \
nshots=901 verbose=2 istart=$istart iend=$iend fmax=90 \
pad=512 niter=33 smooth=8 niterskip=600 shift=16 file_rr=pred_rr_${istart}.su T=1 file_iter=iterT.su
fi
marchenko/demo/oneD/primariesPlane.scr
+ 3
− 3
View file @ 5961b4c7
......@@ -17,11 +17,11 @@ makewave fp=20 dt=0.004 file_out=wave.su nt=1024 t0=0.0 scale=0 scfft=1
#../../marchenko_primaries file_shot=shotsdx5_rp.su plane_wave=1 file_src=wave.su src_angle=10 \
# nshots=901 verbose=2 istart=276 iend=277 fmax=90 pad=1024 \
# xorig=0 niter=20 smooth=10 niterskip=1 file_iter=iterplane.su shift=20 file_rr=plane2_10_rr.su T=0
# xorig=900 niter=2 smooth=10 niterskip=1 file_iter=iterplane.su shift=20 file_rr=plane2_10_rr.su T=0
#exit
../../marchenko_primaries file_shot=shotsdx5_rp.su plane_wave=1 file_src=wave.su src_angle=-15 \
../../marchenko_primaries file_shot=shotsdx5_rp.su plane_wave=1 file_src=wave.su src_angle=-10 \
nshots=901 verbose=2 istart=40 iend=1000 fmax=90 pad=1024 \
niter=20 smooth=10 niterskip=100 shift=20 file_rr=plane2_10_rr.su T=0
niter=20 smooth=10 niterskip=40 niterec=0 shift=20 file_rr=plane2_10_rr.su T=0
exit
......
marchenko/marchenko_primaries.c
+ 38
− 10
View file @ 5961b4c7
......@@ -70,7 +70,8 @@ char *sdoc[] = {
" plane_wave=0 ............. model plane wave",
" src_angle=0 .............. angle with horizontal of plane source array",
" src_velo=1500 ............ velocity to use in src_angle definition",
" xorig=nrecv/2 ............ center of plane-wave",
" t0=0.1 ................... time shift in plane-wave source wavelet for migration",
//" xorig=nrecv/2 ............ center of plane-wave",
" MARCHENKO ITERATIONS ",
" niter=22 ................. number of iterations to initialize and restart",
" niterec=2 ................ number of iterations in recursive part of the time-samples",
......@@ -119,8 +120,8 @@ int main (int argc, char **argv)
float fmin, fmax, tom, deltom, *tapersh, *tapersy, fxf, dxf, *xsrc, *xrcv, *zsyn, *zsrc, *xrcvsyn;
double t0, t1, t2, t3, t4, tsyn, tread, tfft, tcopy, tii;
double energyMi, *energyM0;
float d1, d2, f1, f2, fxsb, fxse, ft, fx, *xsyn, dxsrc;
float *M0, *DD, *RR, dt, dx, dxs, scl, mem;
float tt0, d1, d2, f1, f2, fxsb, fxse, ft, fx, *xsyn, dxsrc;
float *M0, *DD, *RR, *SRC, dt, dx, dxs, scl, mem;
float *rtrace, *tmpdata, *f1min, *f1plus, *RMi, *Mi, *trace;
float *Mup, *Msp;
float xmin, xmax, scale, tsq;
......@@ -168,6 +169,7 @@ int main (int argc, char **argv)
if(!getparint("plane_wave", &plane_wave)) plane_wave=0;
if(!getparfloat("src_angle", &src_angle)) src_angle = 0.0;
if (!getparfloat("src_velo",&src_velo)) src_velo=1500.;
if (!getparfloat("t0",&tt0)) tt0=0.1;
if (T>0) {
T=-1;
......@@ -223,17 +225,24 @@ int main (int argc, char **argv)
fxsb = fx;
}
assert (nxs >= nshots); /* ToDo allow other geometries */
if(!getparint("xorig", &xorig)) xorig=-(nxs-1)/2;
//if(!getparint("xorig", &xorig)) xorig=-(nxs-1)/2;
/* compute time delay for plane-wave responses */
twplane = (float *) calloc(nxs,sizeof(float));
if (plane_wave==1) {
grad2rad = 17.453292e-3;
//p = sin(src_angle*grad2rad)/src_velo;
for (i=0; i<nxs; i++) {
//twplane[i] = (xorig+i)*dxs*p;
twplane[i] = dxs*(xorig+i)*tan(src_angle*grad2rad)/src_velo;
//fprintf(stderr,"plane-wave i=%d x=%f t=%f\n", i, dxs*(xorig+i), twplane[i]);
p = sin(src_angle*grad2rad)/src_velo;
if (p < 0.0) {
for (i=0; i<nxs; i++) {
twplane[i] = fabsf((nxs-i-1)*dxs*p)+tt0;
}
}
else {
for (i=0; i<nxs; i++) {
twplane[i] = (i)*dxs*p+tt0;
//twplane[i] = dxs*(xorig+i)*tan(src_angle*grad2rad)/src_velo;
// fprintf(stderr,"plane-wave i=%d x=%f t=%f %f\n", i, dxs*(xorig+i), twplane[i], (xorig+i)*dxs*p);
}
}
}
......@@ -257,6 +266,7 @@ int main (int argc, char **argv)
Mi = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
M0 = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
DD = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
SRC = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
RR = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
Mup = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
Msp = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
......@@ -435,6 +445,20 @@ int main (int argc, char **argv)
cr1fft(ctrace, rtrace, ntfft, 1);
for (j = 0; j < nts; j++) {
DD[0*nxs*nts+l*nts+j] = -1.0*scl*rtrace[j];
}
/* compute Source wavelet for plane-wave imaging */
for (j = nw_low, m = 0; j <= nw_high; j++, m++) {
tom = j*deltom*twplane[l];
csum.r = cos(-tom);
csum.i = sin(-tom);
//cwav.r = csum.r*cwave[j].r - csum.i*cwave[j].i;
//cwav.i = csum.i*cwave[j].r + csum.r*cwave[j].i;
ctrace[j] = csum;
}
/* transfrom result back to time domain */
cr1fft(ctrace, rtrace, ntfft, 1);
for (j = 0; j < nts; j++) {
SRC[0*nxs*nts+l*nts+j] = scl*rtrace[j];
}
}
}
......@@ -563,7 +587,11 @@ int main (int argc, char **argv)
}
perc=(iend-istart)/100;if(!perc)perc=1;
if (plane_wave) writeDataIter("DDplane.su", DD, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, NINT(src_angle));
if (plane_wave) {
writeDataIter("SRCplane.su", SRC, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, NINT(src_angle));
writeDataIter("DDplane.su", DD, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, NINT(src_angle));
}
free(SRC);
/*================ start loop over number of time-samples ================*/
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment