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

WS15

parent a8515ea8
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
marchenko/demo/WS15/README.3
+ 7
− 3
View file @ b9b9b726
Plane waves
Marchenko emthod worls as well with plane-waves, in stead of modeling a foccus point response for the initial operaror and plane wave is modeled.
Marchenko method works as well with plane-waves, in stead of modeling a foccus point response for the initial operaror and plane wave is modeled.
Submit the initialPlane.scr to the queue to compute the intial operator. The marchenkoPlane.scr script computes the Macrehnko solution for the plane-wave response.
The script epsPlane.scr generates .eps and .png files that can be copied back to your local laptop for display.
DISCUSSION:
In epsPlane.scr there is a convolution between f1minPlane and p0lpus;
- why is this 'needed'?
- does the result have the correct amplitude?
- would this also work with focal points?
marchenko/demo/WS15/README.4
+ 6
− 0
View file @ b9b9b726
......@@ -3,3 +3,9 @@ To generate a shot record with Primaries reflections only submit the demo/twoD/p
The script epsPrimaries.scr generates .eps and .png files that can be copied back to your local laptop for display.
The directory demo/invisible contains the Reflection data to analyse:
shotsdx4_rp.su
Try to find out, by ONLY using the program marchenko_primaries, how many reflectors there are.
marchenko/demo/twoD/epsPlane.scr
+ 27
− 8
View file @ b9b9b726
#!/bin/bash
#intial plane operator
file=iniPlane_z1100_rp.su
file=iniPlane_z1800_rp.su
file_base=${file%.su}
sumax < $file mode=abs outpar=nep
clipref=`cat nep | awk '{print $1/2}'`
supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file \
label1="time (s)" label2="lateral distance (m)" \
n1tic=2 x1beg=0 x1end=2.0 d1num=0.4 \
n1tic=2 x1beg=0 x1end=2.0 d1num=0.4 title="${file_base}" \
f2=-3000 f2num=-3000 d2num=1000 clip=$clipref > ${file_base}.eps
convert -quality 90 -antialias ${file_base}.eps ${file_base}.jpg
supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < p0plus.su \
label1="time (s)" label2="lateral distance (m)" \
n1tic=2 x1beg=0 x1end=2.0 d1num=0.4 \
n1tic=2 x1beg=0 x1end=2.0 d1num=0.4 title="${file_base}_muted" \
f2=-3000 f2num=-3000 d2num=1000 clip=$clipref > ${file_base}_muted.eps
convert ${file_base}_muted.eps ${file_base}_muted.png
convert ${file_base}_muted.eps ${file_base}_muted.jpg
#Marchencko computed Greens function
file=f2Plane.su
#Zero-offset section
suwind key=offset min=0 max=0 < shots/refl_rp.su > zo.su
makewave fp=30 dt=0.004 file_out=wave4.su nt=1024 t0=0.0 scale=1
fconv file_in1=zo.su file_in2=wave4.su file_out=zow.su
file=zow.su
file_base=${file%.su}
sumax < $file mode=abs outpar=nep
clipref=`cat nep | awk '{print $1/4}'`
supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file \
label1="time (s)" label2="lateral distance (m)" \
n1tic=2 x1beg=0 x1end=2.0 d1num=0.4 title="${file_base}" \
f2=-3000 f2num=-3000 d2num=1000 clip=$clipref > ${file_base}.eps
convert ${file_base}.eps ${file_base}.jpg
#Marchencko computed result based on f1min
fconv file_in2=p0plus.su file_in1=f1minPlane.su file_out=convPlane.su verbose=3 mode=conv
fconv file_in1=convPlane.su file_in2=wave4.su file_out=convPlanew.su
file=convPlanew.su
file_base=${file%.su}
sumax < $file mode=abs outpar=nep
clipref=`cat nep | awk '{print $1/4}'`
supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file \
label1="time (s)" label2="lateral distance (m)" \
n1tic=2 x1beg=0 x1end=2.0 d1num=0.4 \
n1tic=2 x1beg=0 x1end=2.0 d1num=0.4 title="${file_base}" \
f2=-3000 f2num=-3000 d2num=1000 clip=$clipref > ${file_base}.eps
convert ${file_base}.eps ${file_base}.png
convert ${file_base}.eps ${file_base}.jpg
marchenko/demo/twoD/initialPlane.scr
+ 5
− 4
View file @ b9b9b726
......@@ -13,7 +13,8 @@ makemod sizex=6000 sizez=2000 dx=$dx dz=$dx cp0=1900 ro0=1200 \
intt=def x=-3000,3000 z=780,780 poly=0 cp=2230 ro=1700 \
intt=def x=-3000,-2200,-1500,0,1300,2100,3000 z=520,580,680,840,680,600,500 poly=2 cp=2400 ro=2800 \
makewave fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
#makewave fp=30 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
makewave w=fw fmin=0 flef=5 frig=80 fmax=100 dt=$dt file_out=wavefw.su nt=4096 t0=0.3 scale=0
dxshot=10
ishot=300
......@@ -22,12 +23,12 @@ nshots=301
export KMP_AFFINITY=disabled
export OMP_NUM_THREADS=16
file_rcv=iniPlane_z1100.su
file_rcv=iniPlane_z1800.su
fdelmodc \
file_cp=synclDown_cp.su ischeme=1 iorder=4 \
file_den=synclDown_ro.su \
file_src=wave.su \
file_src=wavefw.su \
file_rcv=$file_rcv \
src_type=1 \
src_orient=1 \
......@@ -42,7 +43,7 @@ fdelmodc \
dxrcv=10.0 \
xrcv1=-3000 xrcv2=3000 \
zrcv1=0 zrcv2=0 \
xsrc=0 zsrc=1100 plane_wave=1 nsrc=2401 \
xsrc=0 zsrc=1800 plane_wave=1 nsrc=2401 \
ntaper=200 \
left=2 right=2 top=2 bottom=2
......
marchenko/demo/twoD/marchenkoPlane.scr
+ 3
− 4
View file @ b9b9b726
......@@ -5,11 +5,10 @@ export PATH=$HOME/src/OpenSource/bin:$PATH:
export OMP_NUM_THREADS=40
#mute all events below the first arrival to get the intial focusing field
fmute file_shot=iniPlane_z1100_rp.su file_out=p0plus.su above=-1 shift=-10 verbose=1 check=1 hw=4
fmute file_shot=iniPlane_z1800_rp.su file_out=p0plus.su above=-1 shift=-10 verbose=1 check=1 hw=4
#apply the Marchenko algorithm
marchenko file_shot=shots/refl_rp.su file_tinv=p0plus.su nshots=601 verbose=1 \
tap=0 niter=15 hw=8 shift=7 smooth=3 \
file_green=pPlane.su \
file_f1plus=f1plusPlane.su file_f1min=f1minPlane.su file_f2=f2Plane.su
tap=0 niter=15 hw=8 shift=7 smooth=3 rotate=0 \
file_f1plus=f1plusPlane.su file_f1min=f1minPlane.su file_green=greenPlane.su
marchenko/marchenko.c
+ 23
− 17
View file @ b9b9b726
......@@ -105,6 +105,7 @@ char *sdoc[] = {
" file_pplus= .............. output file with p+ ",
" file_pmin= ............... output file with p- ",
" file_iter= ............... output file with -Ni(-t) for each iteration",
" rotate=1 ................. 1: t=0 at nt/2 (middle) 0: t=0 at sample 0 for f1+,-",
" verbose=0 ................ silent option; >0 displays info",
" ",
" ",
......@@ -118,7 +119,7 @@ int main (int argc, char **argv)
FILE *fp_out, *fp_f1plus, *fp_f1min;
FILE *fp_gmin, *fp_gplus, *fp_f2, *fp_pmin;
int i, j, l, ret, nshots, Nfoc, nt, nx, nts, nxs, ngath;
int size, n1, n2, ntap, tap, di, ntraces, pad;
int size, n1, n2, ntap, tap, di, ntraces, pad, rotate;
int nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn;
int reci, countmin, mode, n2out, verbose, ntfft;
int iter, niter, tracf, *muteW, *tsynW;
......@@ -170,6 +171,7 @@ int main (int argc, char **argv)
if (!getparint("tap", &tap)) tap = 0;
if (!getparint("ntap", &ntap)) ntap = 0;
if (!getparint("pad", &pad)) pad = 0;
if (!getparint("rotate", &rotate)) rotate = 1;
if(!getparint("niter", &niter)) niter = 10;
if(!getparint("hw", &hw)) hw = 15;
......@@ -699,14 +701,16 @@ sqrt(energyNi/energyN0));
}
if (file_f1plus != NULL) {
/* rotate to get t=0 in the middle */
for (i = 0; i < n2; i++) {
hdrs_out[i].f1 = -n1*0.5*dt;
memcpy(&trace[0],&f1plus[l*size+i*nts],nts*sizeof(float));
for (j = 0; j < n1/2; j++) {
f1plus[l*size+i*nts+n1/2+j] = trace[j];
}
for (j = n1/2; j < n1; j++) {
f1plus[l*size+i*nts+j-n1/2] = trace[j];
if (rotate==1) {
for (i = 0; i < n2; i++) {
hdrs_out[i].f1 = -n1*0.5*dt;
memcpy(&trace[0],&f1plus[l*size+i*nts],nts*sizeof(float));
for (j = 0; j < n1/2; j++) {
f1plus[l*size+i*nts+n1/2+j] = trace[j];
}
for (j = n1/2; j < n1; j++) {
f1plus[l*size+i*nts+j-n1/2] = trace[j];
}
}
}
ret = writeData(fp_f1plus, (float *)&f1plus[l*size], hdrs_out, n1, n2);
......@@ -714,14 +718,16 @@ sqrt(energyNi/energyN0));
}
if (file_f1min != NULL) {
/* rotate to get t=0 in the middle */
for (i = 0; i < n2; i++) {
hdrs_out[i].f1 = -n1*0.5*dt;
memcpy(&trace[0],&f1min[l*size+i*nts],nts*sizeof(float));
for (j = 0; j < n1/2; j++) {
f1min[l*size+i*nts+n1/2+j] = trace[j];
}
for (j = n1/2; j < n1; j++) {
f1min[l*size+i*nts+j-n1/2] = trace[j];
if (rotate==1) {
for (i = 0; i < n2; i++) {
hdrs_out[i].f1 = -n1*0.5*dt;
memcpy(&trace[0],&f1min[l*size+i*nts],nts*sizeof(float));
for (j = 0; j < n1/2; j++) {
f1min[l*size+i*nts+n1/2+j] = trace[j];
}
for (j = n1/2; j < n1; j++) {
f1min[l*size+i*nts+j-n1/2] = trace[j];
}
}
}
ret = writeData(fp_f1min, (float *)&f1min[l*size], hdrs_out, n1, n2);
......
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