diff --git a/Make_include b/Make_include
index c0b5d9287ba4d9bf222e68245e23c1be08f9ee14..3ad4af8460fd6e3a2bf85f99e06a1c583fb9f71b 100644
--- a/Make_include
+++ b/Make_include
@@ -13,7 +13,7 @@ ROOT=/Users/jan/src/OpenSource
#GNU
CC = gcc
-CC = gcc-mp-6
+#CC = gcc-mp-6
FC = gfortran
# Linux gcc version 4.x
OPTC = -O3 -ffast-math
diff --git a/Makefile b/Makefile
index e8fd8c5013271533077854ad10b5ed36455a589a..62c2fe254f7cebadc153b1a4ddf23d4fd9b5bbe4 100644
--- a/Makefile
+++ b/Makefile
@@ -3,7 +3,6 @@
all: mkdirs
cd FFTlib ; $(MAKE)
cd fdelmodc ; $(MAKE) install
- cd fdemmodc ; $(MAKE) install
cd utils ; $(MAKE) install
cd marchenko ; $(MAKE) install
@@ -15,16 +14,12 @@ mkdirs:
clean:
cd FFTlib ; $(MAKE) $@
cd fdelmodc ; $(MAKE) $@
- cd fdemmodc ; $(MAKE) $@
cd utils ; $(MAKE) $@
cd marchenko ; $(MAKE) $@
-# find . -name "._*"
-
realclean:
cd FFTlib ; $(MAKE) $@
cd fdelmodc ; $(MAKE) $@
- cd fdemmodc ; $(MAKE) $@
cd utils ; $(MAKE) $@
cd marchenko ; $(MAKE) $@
rm -f lib/*
diff --git a/Makefile_all b/Makefile_all
index 663d26925b58f453727b6d89a3e6f0a624862ce9..5e1ab9dbf1f11343cb240c671eb0b4ad89205e64 100644
--- a/Makefile_all
+++ b/Makefile_all
@@ -4,7 +4,9 @@
all: mkdirs
cd FFTlib ; $(MAKE)
cd fdelmodc ; $(MAKE) install
+ cd fdemmodc ; $(MAKE) install
cd utils ; $(MAKE) install
+ cd marchenko ; $(MAKE) install
cd extrap ; $(MAKE)
cd extrap3d ; $(MAKE)
@@ -16,7 +18,9 @@ mkdirs:
clean:
cd FFTlib ; $(MAKE) $@
cd fdelmodc ; $(MAKE) $@
+ cd fdemmodc ; $(MAKE) $@
cd utils ; $(MAKE) $@
+ cd marchenko ; $(MAKE) $@
cd extrap ; $(MAKE) $@
cd extrap3d ; $(MAKE) $@
@@ -24,6 +28,7 @@ realclean:
cd FFTlib ; $(MAKE) $@
cd fdelmodc ; $(MAKE) $@
cd utils ; $(MAKE) $@
+ cd marchenko ; $(MAKE) $@
cd extrap ; $(MAKE) $@
cd extrap3d ; $(MAKE) $@
rm -f lib/*
diff --git a/README b/README
index fc1107a32fefb3eefbecc2cd7995a4ce291ac248..f460387ae60b1c8136864ed93009bbc05b1e9eef 100644
--- a/README
+++ b/README
@@ -9,15 +9,19 @@ Some routines are from Seismic Unix and include the SU LEGAL_STATEMENT in the so
%%%%%%%
-REFERENCE
+REFERENCES
---------
-If this code has helped you in your research please refer to our paper in your publications:
+-1- If the Finite Difference code has helped you in your research please refer to our paper in your publications:
Finite-difference modeling experiments for seismic interferometry
Jan Thorbecke and Deyan Draganov
2011, Geophysics, Vol. 76, no 6 (November-December); p H1--H18, doi: 10.1190/GEO2010-0039.1
-The paper can be downloaded from:
+-2- If the Machenko code has helped you in your research please refer to our paper in your publications:
+
+Hopefully, a published reference to the paper can be put here.
+
+These papers can be downloaded from:
http://janth.home.xs4all.nl/Publications/Publications.html
@@ -27,12 +31,11 @@ INSTALLATION
1) To compile and link the code you first have to set the ROOT variable in the Make_include file which can be found in the directory where you have found this README.
-2) Check the compiler and CFLAGS options in the file Make_include and adapt to the system you are using. The default options are set for a the GNU C-compiler on a Linux system. A Fortran or g++ compiler is not needed to compile the code. The Makefile has been tested with GNU make.
+2) Check the compiler and CFLAGS options in the file Make_include and adapt to the system you are using. The default options are set for a the GNU C-compiler on a Linux system. A Fortran or C++ compiler are not needed to compile the code. The Makefile has been tested with GNU make.
3) If the compiler options are set in the Make_include file you can type
-==> Please make sure that SU is compiled without XDR (in $CWPROOT/Makefile.config make sure that XDRFLAG is NOT set). The SU output
-files of fdelmodc are all based on local IEEE data.
+==> Please make sure that SU is compiled without XDR (in $CWPROOT/Makefile.config make sure that XDRFLAG is NOT set). The SU output files of fdelmodc are all based on local IEEE data.
> make clean
> make
@@ -41,6 +44,7 @@ and the Makefile will make:
- FFT library
- fdelmodc
+- marchenko
- utils
The libraries will be placed in the lib/ directory and the executables in the bin/ directory.
@@ -50,7 +54,10 @@ To use the executables don't forget to include the pathname in your PATH:
export PATH='path_to_this_directory'/bin:$PATH:
setenv PATH 'path_to_this_directory'/bin:$PATH:
-If the compilation has finished without errors and produced an executable called fdelmodc you can run one of the demo programs by running
+
+Finite Difference Modeling: FDELMODC
+------------------------------------
+If the compilation has finished without errors and produced an executable called bin/fdelmodc you can run one of the demo programs by running
> ./fdelmodc_plane.scr
@@ -60,6 +67,20 @@ The demo directory contains scripts which demonstrate the different possibilitie
To reproduce the Figures shown in the GEOPHYICS manuscript "Finite-difference modeling experiments for seismic interferometry" the scripts in FiguresPaper directory can be used. Please read the README in the FiguresPaper directory for more instructions and guidelines.
+
+Marchenko method : MARCHENKO
+----------------------------
+If the compilation has finished without errors and produced an executable called bin/marchenko you can run one of the demo programs by running
+
+> ./
+
+in the directory marchenko/demo/
+
+To reproduce the Figures shown in the GEOPHYICS manuscript "Implementation of the Marchenko method" the scripts in marchenko/demo/oneD directory can be used. Please read the README in the FiguresPaper directory for more instructions and guidelines.
+
+
+MISC
+----
Other make commands which can be useful:
make clean : removes all object files, but leaves libraries and executables
@@ -69,6 +90,5 @@ The latest version of the source code and manual can be found at:
http://www.xs4all.nl/~janth/Software/Software.html
-The code is used by many different people and if there is a request for a new option in the code (for example place receivers in a circle) then I will try to implement, test and make it available.
-
+The code is used by many different people and if there is a request for a new option in the code cle) then I will try to implement, test and make it available.
diff --git a/fdelmodc/demo/fdelmodc_circ_medium.scr b/fdelmodc/demo/fdelmodc_circ_medium.scr
new file mode 100755
index 0000000000000000000000000000000000000000..ba7ebe9a7e42b4ec83c832e72bc4650be549fefc
--- /dev/null
+++ b/fdelmodc/demo/fdelmodc_circ_medium.scr
@@ -0,0 +1,93 @@
+#!/bin/bash
+#PBS -N fdelmod
+#PBS -q verylong
+#PBS -l nodes=1
+#PBS -k eo
+#PBS -j eo
+
+# receivers are placed on a circle
+
+export PATH=../../utils:$PATH:
+
+dt=0.0001
+#makewave file_out=wavelet.su dt=0.0020 nt=1024 fp=13 shift=1 w=g2 verbose=1
+makewave w=fw fmin=0 flef=6 frig=89 fmax=95 dt=$dt file_out=wavefw.su nt=16384 t0=0.3 scale=0 scfft=1
+#sufft < wavefw.su | suamp | sugain scale=0.0001 | suxgraph style=normal
+
+makemod file_base=model.su \
+ cp0=1500 ro0=1000 sizex=4000 sizez=4000 dx=2 dz=2 orig=-2000,-2000 \
+ intt=elipse var=1000,1000 x=0 z=0 cp=2000 ro=1000 verbose=1
+
+makemod file_base=hom.su \
+ cp0=1500 ro0=1000 sizex=4000 sizez=4000 dx=2 dz=2 orig=-2000,-2000 verbose=1
+
+export filecp=model_cp.su
+export filero=model_ro.su
+file_rcv=circle.su
+
+export OMP_NUM_THREADS=2
+
+set -x
+#../fdelmodc \
+ file_cp=$filecp file_den=$filero \
+ ischeme=1 \
+ file_src=wavefw.su verbose=2 \
+ file_rcv=$file_rcv \
+ rec_type_vz=0 rec_type_p=1 \
+ xrcv1=-2000 xrcv2=2000 zrcv1=-2000 zrcv2=-2000 \
+ dxrcv=10 \
+ dtrcv=0.004 \
+ xsrc=0 zsrc=-2000 nshot=1 nsrc=1 \
+ src_type=1 tmod=4.092 \
+ ntaper=100 \
+ left=2 right=2 bottom=2 top=2
+
+export filecp=hom_cp.su
+export filero=hom_ro.su
+file_rcv=hom.su
+
+../fdelmodc \
+ file_cp=$filecp file_den=$filero \
+ ischeme=1 \
+ file_src=wavefw.su verbose=2 \
+ file_rcv=$file_rcv \
+ rec_type_vz=0 rec_type_p=1 \
+ xrcv1=-2000 xrcv2=2000 zrcv1=-2000 zrcv2=-2000 \
+ dxrcv=10 \
+ dtrcv=0.004 \
+ xsrc=0 zsrc=-2000 nshot=1 nsrc=1 \
+ src_type=1 tmod=4.092 \
+ ntaper=100 \
+ left=2 right=2 bottom=2 top=2
+
+exit
+
+#for exptype in circle square doodle ;
+for exptype in circle ;
+do
+for rectype in rvx rvz rp ;
+do
+ file_rcv=${exptype}_$rectype.su
+ echo $file_rcv
+ supsimage < ${exptype}_$rectype.su hbox=4 wbox=3 titlesize=-1 labelsize=10 titlesize=-1 \
+ perc=99 label1="time [s]" f2=0 d2=2 label2="rotation in degrees" > ${exptype}_$rectype.eps
+
+done
+done
+
+supsimage < model_cp.su \
+ wbox=4 hbox=4 titlesize=-1 labelsize=10 verbose=1 \
+ d2=3 f2=0 wrgb=1.0,0,0 grgb=0,1.0,0 brgb=0,0,1.0 bps=24 \
+ label1="depth [m]" label2="lateral position [m]" > model_plane.eps
+
+supsimage < SrcRecPositions.su \
+ wbox=4 hbox=4 titlesize=-1 labelsize=10 verbose=1 \
+ d2=3 f2=0 wclip=-1 bclip=1 \
+ gabel1="depth [m]" label2="lateral position [m]" > SrcRecCircle.eps
+
+suop2 model_cp.su SrcRecPositions.su w1=1 w2=2000 op=sum | \
+ supsimage wclip=1400 bclip=2000 \
+ wbox=4 hbox=4 titlesize=-1 labelsize=10 verbose=1 \
+ d2=3 f2=0 wrgb=1.0,0,0 grgb=0,1.0,0 brgb=0,0,1.0 bps=24 \
+ label1="depth [m]" label2="lateral position [m]" > model_plane_src.eps
+
diff --git a/fdelmodc/demo/fdelmodc_fault.scr b/fdelmodc/demo/fdelmodc_fault.scr
new file mode 100755
index 0000000000000000000000000000000000000000..e2d06417aeb42c0346951b658ce69b080d896303
--- /dev/null
+++ b/fdelmodc/demo/fdelmodc_fault.scr
@@ -0,0 +1,48 @@
+#!/bin/bash
+#PBS -l nodes=1
+#PBS -N InterfModeling
+#PBS -q fourweeks
+#PBS -V
+#
+export PATH=$HOME/../thorbcke/src/OpenSource/bin:$PATH
+
+dt=0.0015
+fmax=45
+cp=2000
+
+makemod sizex=6000 sizez=2000 dx=5 dz=5 cp0=$cp ro0=1000 file_base=fault.su orig=0,0 \
+ intt=def poly=2 x=0,600,1000,2000,4000,5000 \
+ z=550,550,550,300,300,500 cp=$cp ro=2000 \
+ intt=def poly=0 x=5000,6000 z=300,200 cp=$cp ro=1500 \
+ intt=def poly=0 x=0,2500 z=900,900 cp=$cp ro=2200 \
+ intt=def poly=0 x=2000,5000 z=1000,300 cp=$cp ro=1500 \
+ intt=def poly=0 x=2000,3000,6000 z=1000,770,770 cp=$cp ro=1800 \
+ intt=def poly=0 x=2000,6000 z=1000,1000 cp=$cp ro=2200 \
+ intt=def poly=0 x=0,6000 z=1400,1400 cp=$cp ro=2400 \
+
+suximage cmap=hsv4 < fault_ro.su &
+
+makewave w=g2 fmax=45 t0=0.10 dt=$dt nt=4096 db=-40 file_out=G2.su verbose=1
+
+../fdelmodc \
+ file_cp=fault_cp.su ischeme=1 \
+ file_den=fault_ro.su \
+ file_rcv=shots.su \
+ file_src=G2.su \
+ src_type=1 \
+ dtrcv=0.003 \
+ verbose=1 \
+ tmod=3.104 \
+ nshot=101 \
+ dxshot=60 \
+ rec_delay=0.1 \
+ rec_type_vz=0 \
+ dxrcv=60.0 \
+ xsrc=0 \
+ zsrc=0 \
+ ntaper=101 \
+ left=2 right=2 top=2 bottom=2
+
+suwind key=offset min=0 max=0 < shots_rp.su | suzero itmax=30 | suximage x1end=2
+
+
diff --git a/fdelmodc/threadAffinity.c b/fdelmodc/threadAffinity.c
index 31593ceb62e58382ebf6bc178d1d26aaaa5dae2a..dcda8b6301993f52e79932e9dbd7da2d13cae501 100644
--- a/fdelmodc/threadAffinity.c
+++ b/fdelmodc/threadAffinity.c
@@ -13,6 +13,7 @@
#define CPU_SETSIZE 1024
#define SYSCTL_CORE_COUNT "machdep.cpu.core_count"
+void vmess(char *fmt, ...);
typedef struct cpu_set {
uint32_t count;
@@ -100,8 +101,7 @@ void threadAffinity(void)
#endif
(void)sched_getaffinity(0, sizeof(coremask), &coremask);
cpuset_to_cstr(&coremask, clbuf);
- printf("%s thread %d, on %s. (core affinity = %s)\n",
- prefix, thread, hnbuf, clbuf);
+ vmess("%s thread %d, on %s. (core affinity = %s)", prefix, thread, hnbuf, clbuf);
}
return;
diff --git a/marchenko/Makefile b/marchenko/Makefile
index 1cc4fd48224e7abde39aec1f18acbb4291e89de8..1a90e73a4c29a85c014e739ce5194ec7c9437b87 100644
--- a/marchenko/Makefile
+++ b/marchenko/Makefile
@@ -3,10 +3,11 @@
include ../Make_include
LIBS += -L$L -lgenfft -lm $(LIBSM)
-#OPTC += -openmp
#OPTC += -g -O0
-ALL: fmute marchenko marchenko2
+#ALL: fmute marchenko marchenko2
+
+ALL: fmute marchenko
SRCJ = fmute.c \
getFileInfo.c \
@@ -61,10 +62,11 @@ OBJH2 = $(SRCH2:%.c=%.o)
marchenko2: $(OBJH2)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko2 $(OBJH2) $(LIBS)
-install: fmute marchenko marchenko2
+install: fmute marchenko
cp fmute $B
cp marchenko $B
- cp marchenko2 $B
+
+# cp marchenko2 $B
clean:
rm -f core fmute $(OBJJ) marchenko $(OBJH) marchenko2 $(OBJH2)
diff --git a/marchenko/demo/README b/marchenko/demo/README
index 5fc50362f6ee58da1104668356b8ac28380d6013..f5a7c129e2a168d2325d162c9f886a2cc27ae43b 100644
--- a/marchenko/demo/README
+++ b/marchenko/demo/README
@@ -1,9 +1,4 @@
-Description of files:
-1) shots.scr create the shots
-2) model.scr computes the model
-3) direct_wave.scr crate the direct wave to be removed from the shots
-4) remove_direct.scr remove the direct wave from the shots and scale them
-5) first_arrival.scr computes the first arrival
-6) marchenko.scr perform the Marchenko scheme
-7) referenceShot.scr creates the reference Green's function
+The scripts to reproduce the Figures in the manuscript can be found in the directory oneD. The oneD/README explains how to run the scripts.
+
+A more complicated model can be found in the directory twoD and will takes several hours to model the reflection data.
diff --git a/marchenko/demo/oneD/.DS_Store b/marchenko/demo/oneD/.DS_Store
new file mode 100644
index 0000000000000000000000000000000000000000..5008ddfcf53c02e82d7eee2e57c38e5672ef89f6
Binary files /dev/null and b/marchenko/demo/oneD/.DS_Store differ
diff --git a/marchenko/demo/oneD/.eps b/marchenko/demo/oneD/.eps
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/marchenko/demo/oneD/README b/marchenko/demo/oneD/README
new file mode 100644
index 0000000000000000000000000000000000000000..17754221d365dd7561433df7abac8110f8b41c89
--- /dev/null
+++ b/marchenko/demo/oneD/README
@@ -0,0 +1,180 @@
+Description of files:
+1) model.scr computes the model and the 'basis' shot of R => shot5_rp.su
+2) p5all.scr create from basis shot full Reflection response matrix => shotsdx5_rp.su (3.3 GB)
+3) initialFocus.scr model G_d the intitiAla focusing function => iniFocus_rp.su
+4) referenceShot.scr creates the reference Green's function at focal point => referenceP_rp.su
+5) marchenko.scr perform the Marchenko scheme => pgreen.su, f1plus0.su, f1min0.su, f2.su
+
+extra scripts
+marchenkoIter.scr to make the figure with "Four iterations of the Marchenko method."
+
+backpropf2.scr to make Figure "Snapshots of back-propagation of f_2."
+
+eps*.scr reproduce the postscript files of the manuscript using SU postscript plotting programs.
+
+
+To reproduce the Figures in the Manuscript:
+
+--------------------------
+* Figure 2: Wavelet
+* Figure 3: Model + Initial wavefield
+
+run model.scr to generate the data .su files: this will take 3-4 minutes. The files generate are:
+ - hom_cp.su, hom_ro.su
+ - model10_cp.su, model10_ro.su
+ - shot5_fd_rp.su
+ - shot5_hom_fd_rp.su
+ - shot5_rp.su
+ - wave.su
+ - wavefw.su
+
+run initialFocus.scr to compute the direct arrival of the transmission response G_d. This will take 1-2 minutes.
+ - modelup_cp.su
+ - modelup_ro.su
+ - iniFocus_rp.su
+
+run epsModel.scr to generate the postscript files of Figure 2 and 3
+
+wavefw.eps => Figure 2a
+wavefw_freq.eps => Figure 2b
+
+model_cp_line.eps => Figure 3a
+model_ro_line.eps => Figure 3b
+shotx0_rp.eps => Figure 3c
+iniFocus_rp.eps => Figure 3d
+
+
+--------------------------
+* Figure 4: Initialisation
+* Figure 5: first update
+* Figure 6: first 4 iterations
+
+The full R matrix is build up from the the shot record computed with model.scr
+
+run p5all.scr to generate the full R matrix for a fixed spread geometry. This will take less than one minute. The file generated is
+ - shotsdx5_rp.su the file has a size of 3.3 GB
+
+This R, together with iniFocus_rp.su, is the input of the Marchenko algorithm
+
+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 Compute the reference Green's function at x=0 z=900 m in the actual model
+run referenceShot.scr This will take 1 minute.
+ - referenceP_rp.su
+
+To generate all postscript files for Figure 4, 5 and 6
+
+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
+iter_002.eps => Figure 5 -N_1
+f1plus_002.eps => Figure 5 f^+_1,0
+
+-- Figure 6 column 1
+iter_001.eps
+iter_002.eps
+iter_003.eps
+iter_004.eps
+-- Figure 6 column 2
+f1min_001.eps
+f1min_002.eps
+f1min_003.eps
+f1min_004.eps
+-- Figure 6 column 3
+p0plus_flip.eps
+f1plus_002.eps
+f1plus_003.eps
+f1plus_004.eps
+-- Figure 6 column 4
+pgreen_001.eps
+pgreen_002.eps
+pgreen_003.eps
+pgreen_004.eps
+-- Figure 6 column 5
+compare_001.eps
+compare_002.eps
+compare_003.eps
+compare_004.eps
+
+--------------------------
+* Figure 7: Comparison of Marchenko result with reference
+
+To compute the marchenko results for 8 iterations.
+
+run marchenko.scr This will take less than 1 minute. The generated files are:
+ - pgreen.su, pgreen512.su
+ - diffref.su
+ - Gplus0.su
+ - Gmin0.su
+ - f1plus0.su
+ - f1min0.su
+ - f2.su
+
+
+At the end of the run the script will display in X11 a comparison of the middle trace.
+
+To make the postscript figure
+
+run epsCompare.scr
+
+mergeGreenRef.eps => Figure 7
+
+--------------------------
+* Figure 8: snapshots of back propagating f2 in actual medium
+
+To compute the snapshots
+
+run backpropf2.scr This will take about 2 minutes. The generated output files are
+ - backpropf2_sp.su
+ - backpropf2detail_sp.su
+
+The postscript files of Figure 8 are generated with
+
+epsBackprop.scr
+
+
+-- Figure 8 column 1
+backpropf2_-0.30.eps
+backpropf2_-0.15.eps
+backpropf2_0.0.eps
+backpropf2_0.15.eps
+backpropf2_0.30.eps
+-- Figure 8 column 2
+backpropf2_-0.03.eps
+backpropf2_-0.02.eps
+backpropf2_0.00.eps
+backpropf2_0.02.eps
+backpropf2_0.03.eps
+-- Figure 8 column 3
+backpropf2sum_0.03.eps
+backpropf2sum_0.15.eps
+backpropf2sum_0.30.eps
+
+
diff --git a/marchenko/demo/oneD/backpropf2.scr b/marchenko/demo/oneD/backpropf2.scr
new file mode 100755
index 0000000000000000000000000000000000000000..e1e614c8e8b7569a8d61c1d1f48ec512358a39a8
--- /dev/null
+++ b/marchenko/demo/oneD/backpropf2.scr
@@ -0,0 +1,90 @@
+#!/bin/bash
+
+export PATH=$HOME/src/OpenSource/bin/:$PATH:
+
+cd /Users/jan/LaTeX/Articles/Marchenko/scripts/oneD
+which ftr1d
+
+#dx=1.25
+#dt=0.00025
+
+dx=2.5
+dt=0.0005
+
+file_cp=model10_cp.su
+file_ro=model10_ro.su
+
+export OMP_NUM_THREADS=4
+
+# t=0 focal time is at 2.0445 seconds back=propagating
+#shift f2.su such that t=0 is positioned in the middle of the time axis
+# the extra shift of 0.000250 is needed because of the staggered time implementation of the Finite Difference program.
+ns=1024
+dtrcv=`surange < f2.su | grep dt | awk '{print $2/1000000.0}'`
+suwind key=gx min=-2250000 max=2250000 itmax=1023 < f2.su > nep.su
+shift=$(echo "scale=6; ($dtrcv*($ns/2.0-1)+0.5*$dt-0.000250)" | bc -l)
+echo $shift
+basop choice=shift shift=$shift file_in=nep.su verbose=1 > pplus.su
+
+# the f2.su is sampled with 4ms the FD program need 0.5ms
+# time axis is interpolated by making use of FFT's: sinc interpolation
+ftr1d file_in=pplus.su file_out=freq.su
+sushw < freq.su key=nhs,dt a=8192,500 > fr.su
+ftr1d file_in=fr.su n1=8194 file_out=pplusdt.su verbose=1
+
+#backpropagate f2.su and collect snapshots
+fdelmodc \
+ file_cp=$file_cp ischeme=1 iorder=4 \
+ file_den=$file_ro \
+ file_src=pplusdt.su \
+ file_rcv=backprop_f2_z900.su \
+ grid_dir=0 \
+ src_type=1 \
+ src_injectionrate=1 \
+ src_orient=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.004 \
+ rec_delay=0.0 \
+ verbose=2 \
+ tmod=3.10 \
+ dxrcv=5.0 \
+ xrcv1=-2250 xrcv2=2250 \
+ zrcv1=900 zrcv2=900 \
+ zsrc=0 xsrc=0 \
+ npml=101 \
+ file_snap=backpropf2.su tsnap1=1.0445 dtsnap=0.010 tsnap2=3.0445 dxsnap=5 dzsnap=5 zsnap1=0 zsnap2=1250 xsnap1=-1000 xsnap2=1000 \
+ sna_type_vz=0 \
+ sna_type_p=1 \
+ left=2 right=2 top=2 bottom=2
+
+
+#for detailed, middle column,
+#backpropagate f2.su and collect snapshots
+fdelmodc \
+ file_cp=$file_cp ischeme=1 iorder=4 \
+ file_den=$file_ro \
+ file_src=pplusdt.su \
+ file_rcv=backprop_f2_z900.su \
+ grid_dir=0 \
+ src_type=1 \
+ src_injectionrate=1 \
+ src_orient=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.004 \
+ rec_delay=0.0 \
+ verbose=2 \
+ tmod=3.10 \
+ dxrcv=5.0 \
+ xrcv1=-2250 xrcv2=2250 \
+ zrcv1=900 zrcv2=900 \
+ zsrc=0 xsrc=0 \
+ npml=101 \
+ file_snap=backpropf2detail.su tsnap1=2.0245 dtsnap=0.0005 tsnap2=2.0645 dxsnap=2.5 dzsnap=2.5 zsnap1=0 zsnap2=1250 xsnap1=-1000 xsnap2=1000 \
+ sna_type_vz=0 \
+ sna_type_p=1 \
+ left=2 right=2 top=2 bottom=2
+
diff --git a/marchenko/demo/oneD/clean b/marchenko/demo/oneD/clean
new file mode 100755
index 0000000000000000000000000000000000000000..3890128152ba3f4b11471dfdb5ddd1399840bc08
--- /dev/null
+++ b/marchenko/demo/oneD/clean
@@ -0,0 +1,4 @@
+#!/bin/bash
+
+rm *.su *.bin *.eps nep line* *.asci
+
diff --git a/marchenko/demo/oneD/eps.scr b/marchenko/demo/oneD/eps.scr
new file mode 100755
index 0000000000000000000000000000000000000000..37fe990f319a60e196dce785d83c02d23e6edd5c
--- /dev/null
+++ b/marchenko/demo/oneD/eps.scr
@@ -0,0 +1,186 @@
+#!/bin/bash
+
+export PATH=$HOME/src/OpenSource/bin/:$PATH:
+
+#Postscript file of model and shot record at xsrc=0
+cat << EOF1 > line1
+400 -2500
+400 2500
+EOF1
+
+cat << EOF2 > line2
+700 -2500
+700 2500
+EOF2
+
+cat << EOF3 > line3
+1100 -2500
+1100 2500
+EOF3
+
+#model
+supsimage hbox=4 wbox=6 labelsize=12 < model10_cp.su \
+ x1beg=0 x1end=1400.0 d1num=200 legend=1 threecolor=0 \
+ curve=line1,line2,line3 npair=2,2,2 curvecolor=black,black,black \
+ n1tic=5 x2beg=-2250 f2num=-2000 d2num=1000 x2end=2250 > model_cp_line.eps
+
+supsimage hbox=4 wbox=6 labelsize=12 < model10_ro.su \
+ x1beg=0 x1end=1400.0 d1num=200 legend=1 threecolor=0 \
+ curve=line1,line2,line3 npair=2,2,2 curvecolor=black,black,black \
+ n1tic=5 x2beg=-2250 f2num=-2000 d2num=1000 x2end=2250 > model_ro_line.eps
+
+#wavelet
+dt=0.0005
+supsgraph < wavefw.su \
+ labelsize=12 d1=$dt style=normal \
+ label1="time (s)" label2="amplitude" \
+ d1num=0.15 wbox=6 hbox=3 x1end=0.9 > wavefw.eps
+
+sufft < wavefw.su | suamp | sugain scale=$dt | supsgraph \
+ labelsize=12 style=normal \
+ label1="time (s)" label2="amplitude" \
+ d1num=10 wbox=6 hbox=3 x1end=125 x2end=1.1 > wavefw_freq.eps
+
+
+#shot record
+file=shot5_rp.su
+sumax < $file mode=abs outpar=nep
+clipref=`cat nep | awk '{print $1/3}'`
+suwind key=gx min=-2250000 max=2250000 < $file | \
+ supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 \
+ n1tic=2 d2=5 f1=0.0 x1beg=0 x1end=2.004 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > shotx0_rp.eps
+
+#Initial focusing operator
+file=iniFocus_rp.su
+file_base=${file%.su}
+sumax < $file mode=abs outpar=nep
+clipref=`cat nep | awk '{print $1/3}'`
+suwind key=gx min=-2250000 max=2250000 < $file | \
+ supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 \
+ n1tic=2 d2=5 f1=0.0 x1beg=0 x1end=2.004 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
+
+file=p0plus.su
+file_base=${file%.su}
+ns=1024
+dtrcv=`surange < p0plus.su | grep dt | awk '{print $2/1000000.0}'`
+suwind key=gx min=-2250000 max=2250000 itmax=1023 < $file > nep.su
+shift=$(echo "scale=4; ($dtrcv*($ns/2.0-1))" | bc -l)
+basop choice=shift shift=$shift file_in=nep.su | \
+ suflip flip=3 | \
+ supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 \
+ n1tic=2 d2=5 x1beg=-1.504 x1end=1.5 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > ${file_base}_flip.eps
+rm nep.su
+
+file=p0plus.su
+file_base=${file%.su}
+suwind key=gx min=-2250000 max=2250000 < $file | \
+ suflip flip=3 | \
+ supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 \
+ n1tic=2 d2=5 f1=-2.044 f1num=-2.000 x1beg=-2.004 x1end=0.0 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
+
+#convolution of G_d with middle shot from R - not used -
+#suwind key=gx min=-2250000 max=2250000 < shot5_rp.su > shot0.su
+#fconv file_in1=iniFocus_rp.su file_in2=shot0.su file_out=GdRconv.su
+
+#mute to get pslinepos.asci files used in plotting only
+fmute file_shot=iniFocus_rp.su file_out=nep.su above=0 shift=8 verbose=1 check=1 hw=4
+
+#set clip for iteration updates
+file=iter_001.su
+sumax < $file mode=abs outpar=nep
+clipiter=`cat nep | awk '{print $1/6}'`
+
+#iterations
+for (( iter=1; iter<=4; iter+=1 ))
+do
+piter=$(printf %03d $iter)
+echo $piter
+
+file=iter_$piter.su
+#ns=`surange < iter_001.su | grep ns | awk '{print $2}'`
+#dtrcv=`surange < iter_001.su | grep dt | awk '{print $2/1000000.0}'`
+#shift=$(echo "scale=4; ($dtrcv*($ns/2.0-1))" | bc -l)
+#basop choice=shift shift=$shift file_in=$file | \
+file_base=${file%.su}
+#sumax < $file mode=abs outpar=nep
+#clipref=`cat nep | awk '{print $1/3}'`
+clipref=$clipiter
+supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file \
+ n1tic=2 d2=5 x1beg=-1.504 x1end=1.5 d1num=0.4 \
+ curve=pslinepos.asci,pslineneg.asci npair=901,901 curvewidth=2,2 curvecolor=black,black curvedash=3,3 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
+
+file=f1min_$piter.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\
+ n1tic=2 d2=5 x1beg=-1.504 x1end=1.5 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
+
+file=f1plus_$piter.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\
+ n1tic=2 d2=5 x1beg=-1.504 x1end=1.5 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
+
+file=pgreen_$piter.su
+file_base=${file%.su}
+sumax < $file mode=abs outpar=nep
+clipref=`cat nep | awk '{print $1/3}'`
+suwind key=gx min=-2250000 max=2250000 < $file | \
+ supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 \
+ n1tic=2 d2=5 f1=0.0 x1beg=0 x1end=2.004 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
+
+#compare Green's funtions on Marhcenko and reference result
+suwind key=gx min=0 max=0 itmax=511 < pgreen_$piter.su | sumax mode=abs outpar=nepmg
+suwind key=gx min=0 max=0 itmax=511 < referenceP_rp.su | sumax mode=abs outpar=neprf
+mg=`cat nepmg | awk '{print $1}'`
+rf=`cat neprf | awk '{print $1}'`
+value=${value/[eE][+][0]/*10^}
+mg=${mg/[eE][+][0]/*10^}
+rf=${rf/[eE][+][0]/*10^}
+rm nep*
+scale=$(echo "scale=3; ($rf)/($mg)" | bc -l)
+#echo $scale
+
+(suwind key=gx min=0 max=0 itmax=511 < pgreen_$piter.su | sugain scale=$scale; \
+ suwind key=gx min=0 max=0 < referenceP_rp.su) | \
+ supsgraph hbox=6 wbox=1 labelsize=10 linegray=0.3,0.0 style=seismic \
+ lineon=3.5,1.0 lineoff=1.3,0.0 > compare_$piter.eps
+
+done
+
+file=diffref.su
+file_base=${file%.su}
+sumax < referenceP_rp.su mode=abs outpar=nep
+clipref=`cat nep | awk '{print $1}'`
+suwind key=gx min=-2250000 max=2250000 < $file | \
+ supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 \
+ n1tic=2 d2=5 f1=0.0 x1beg=0 x1end=2.004 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
+suwind key=gx min=-2250000 max=2250000 < referenceP_rp.su | \
+ supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 \
+ n1tic=2 d2=5 f1=0.0 x1beg=0 x1end=2.004 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > referenceP_rp.eps
+
+suwind < pgreen512.su j=50 s=1 | \
+ supswigp n2=19 fill=0 \
+ hbox=4 wbox=8 labelsize=10 linewidth=1.0 \
+ n1tic=2 d2=250 f1=0.0 x1beg=0 x1end=2.004 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=500 > green.eps
+suwind < referenceP_rp.su j=50 s=1 | \
+ supswigp n2=19 fill=0 tracecolor=red \
+ hbox=4 wbox=8 labelsize=10 linewidth=1.0 \
+ n1tic=2 d2=250 f1=0.0 x1beg=0 x1end=2.004 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=500 > ref.eps
+
+psmerge in=green.eps in=ref.eps > mergeGreenRef.eps
+
diff --git a/marchenko/demo/oneD/epsBackprop.scr b/marchenko/demo/oneD/epsBackprop.scr
new file mode 100755
index 0000000000000000000000000000000000000000..a4b4f17fa8a32c90bf368c8dd6d867f2e0485582
--- /dev/null
+++ b/marchenko/demo/oneD/epsBackprop.scr
@@ -0,0 +1,54 @@
+#!/bin/bash
+
+export PATH=$HOME/src/OpenSource/bin/:$PATH:
+
+dx=2.5
+dt=0.0005
+
+file_snap="backpropf2"
+dtsnap=0.01
+nsnap=101
+for fldr in 71 86 98 99 101 103 104 116 131;
+do
+ times=$(echo "scale=2; $dtsnap*(${fldr}-$nsnap)" | bc -l)
+ atime=`printf "%4.2f" $times`
+ suwind key=fldr min=$fldr max=$fldr < ${file_snap}_sp.su | \
+ supsimage hbox=4 wbox=6 labelsize=10 \
+ x1beg=0 x1end=1250.0 clip=1e9 \
+ curve=line1,line2,line3 npair=2,2,2 curvecolor=black,black,black curvedash=3,3,3 \
+ n1tic=4 f2=-1000 d2=$dx x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > ${file_snap}_$atime.eps
+done
+
+#for the detailed, middle column, replace this line with the file_snap line above.
+file_snap="backpropf2detail"
+dtsnap=0.0005
+nsnap=41
+for fldr in 1 21 41 61 81;
+do
+ times=$(echo "scale=2; $dtsnap*(${fldr}-$nsnap)" | bc -l)
+ atime=`printf "%4.2f" $times`
+ suwind key=fldr min=$fldr max=$fldr < ${file_snap}_sp.su | \
+ supsimage hbox=4 wbox=6 labelsize=10 \
+ x1beg=0 x1end=1250.0 clip=1e9 \
+ curve=line1,line2,line3 npair=2,2,2 curvecolor=black,black,black curvedash=3,3,3 \
+ n1tic=4 f2=-1000 d2=$dx x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > ${file_snap}_$atime.eps
+done
+
+#select files for snapshot between -0.7 => 0 <= +0.07 (fldr 31-101-171)
+#add pos and negative times to get response of homogenoeus Green's function
+
+file_snap="backpropf2"
+for fldr in 71 86 98 101;
+do
+ times=$(echo "scale=2; -0.01*(${fldr}-101)" | bc -l)
+ atime=`printf "%4.2f" $times`
+ suwind key=fldr min=$fldr max=$fldr < ${file_snap}_sp.su > neg.su
+ (( fldr = 101+(101-$fldr) ))
+ suwind key=fldr min=$fldr max=$fldr < ${file_snap}_sp.su > pos.su
+ susum neg.su pos.su | \
+ supsimage hbox=4 wbox=6 labelsize=10 \
+ x1beg=0 x1end=1250.0 clip=1e9 \
+ curve=line1,line2,line3 npair=2,2,2 curvecolor=black,black,black curvedash=3,3,3 \
+ n1tic=4 x2beg=-1000 d2=$dx f2num=-1000 d2num=500 x2end=1000 > ${file_snap}sum_$atime.eps
+done
+
diff --git a/marchenko/demo/oneD/epsCompare.scr b/marchenko/demo/oneD/epsCompare.scr
new file mode 100755
index 0000000000000000000000000000000000000000..facd0e03845f5409ac27dfd03a3fa76c56e629d1
--- /dev/null
+++ b/marchenko/demo/oneD/epsCompare.scr
@@ -0,0 +1,32 @@
+#!/bin/bash
+
+export PATH=$HOME/src/OpenSource/bin/:$PATH:
+
+#mke figures for reference and Marchenko result an merge into one file
+
+file=diffref.su
+file_base=${file%.su}
+sumax < referenceP_rp.su mode=abs outpar=nep
+clipref=`cat nep | awk '{print $1}'`
+suwind key=gx min=-2250000 max=2250000 < $file | \
+ supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 \
+ n1tic=2 d2=5 f1=0.0 x1beg=0 x1end=2.004 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
+suwind key=gx min=-2250000 max=2250000 < referenceP_rp.su | \
+ supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 \
+ n1tic=2 d2=5 f1=0.0 x1beg=0 x1end=2.004 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > referenceP_rp.eps
+
+suwind < pgreen512.su j=50 s=1 | \
+ supswigp n2=19 fill=0 \
+ hbox=4 wbox=8 labelsize=10 linewidth=1.0 \
+ n1tic=2 d2=250 f1=0.0 x1beg=0 x1end=2.004 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=500 > green.eps
+suwind < referenceP_rp.su j=50 s=1 | \
+ supswigp n2=19 fill=0 tracecolor=red \
+ hbox=4 wbox=8 labelsize=10 linewidth=1.0 \
+ n1tic=2 d2=250 f1=0.0 x1beg=0 x1end=2.004 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=500 > ref.eps
+
+psmerge in=green.eps in=ref.eps > mergeGreenRef.eps
+
diff --git a/marchenko/demo/oneD/epsMarchenkoIter.scr b/marchenko/demo/oneD/epsMarchenkoIter.scr
new file mode 100755
index 0000000000000000000000000000000000000000..35197df60210aa69e9e45f696faa45ba2909806c
--- /dev/null
+++ b/marchenko/demo/oneD/epsMarchenkoIter.scr
@@ -0,0 +1,100 @@
+#!/bin/bash
+
+export PATH=$HOME/src/OpenSource/bin/:$PATH:
+
+#Direct field of transmission repsponse
+file=p0plus.su
+file_base=${file%.su}
+ns=1024
+dtrcv=`surange < p0plus.su | grep dt | awk '{print $2/1000000.0}'`
+suwind key=gx min=-2250000 max=2250000 itmax=1023 < $file > nep.su
+shift=$(echo "scale=4; ($dtrcv*($ns/2.0-1))" | bc -l)
+basop choice=shift shift=$shift file_in=nep.su | \
+ suflip flip=3 | \
+ supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 \
+ n1tic=2 d2=5 x1beg=-1.504 x1end=1.5 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > ${file_base}_flip.eps
+rm nep.su
+
+file=p0plus.su
+file_base=${file%.su}
+suwind key=gx min=-2250000 max=2250000 < $file | \
+ suflip flip=3 | \
+ supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 \
+ n1tic=2 d2=5 f1=-2.044 f1num=-2.000 x1beg=-2.004 x1end=0.0 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
+
+#convolution of G_d with middle shot from R - not used -
+#suwind key=gx min=-2250000 max=2250000 < shot5_rp.su > shot0.su
+#fconv file_in1=iniFocus_rp.su file_in2=shot0.su file_out=GdRconv.su
+
+#mute to get pslinepos.asci files used in plotting only
+fmute file_shot=iniFocus_rp.su file_out=nep.su above=0 shift=8 verbose=1 check=1 hw=4
+
+#set same clip factor for iteration updates
+file=iter_001.su
+sumax < $file mode=abs outpar=nep
+clipiter=`cat nep | awk '{print $1/6}'`
+
+#iterations
+for (( iter=1; iter<=4; iter+=1 ))
+do
+piter=$(printf %03d $iter)
+echo $piter
+
+file=iter_$piter.su
+#ns=`surange < iter_001.su | grep ns | awk '{print $2}'`
+#dtrcv=`surange < iter_001.su | grep dt | awk '{print $2/1000000.0}'`
+#shift=$(echo "scale=4; ($dtrcv*($ns/2.0-1))" | bc -l)
+#basop choice=shift shift=$shift file_in=$file | \
+file_base=${file%.su}
+clipref=$clipiter
+supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file \
+ n1tic=2 d2=5 x1beg=-1.504 x1end=1.5 d1num=0.4 \
+ curve=pslinepos.asci,pslineneg.asci npair=901,901 curvewidth=2,2 curvecolor=black,black curvedash=3,3 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
+
+file=f1min_$piter.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\
+ n1tic=2 d2=5 x1beg=-1.504 x1end=1.5 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
+
+file=f1plus_$piter.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\
+ n1tic=2 d2=5 x1beg=-1.504 x1end=1.5 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
+
+file=pgreen_$piter.su
+file_base=${file%.su}
+sumax < $file mode=abs outpar=nep
+clipref=`cat nep | awk '{print $1/3}'`
+suwind key=gx min=-2250000 max=2250000 < $file | \
+ supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 \
+ n1tic=2 d2=5 f1=0.0 x1beg=0 x1end=2.004 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
+
+#compare Green's funtions on Marhcenko and reference result
+suwind key=gx min=0 max=0 itmax=511 < pgreen_$piter.su | sumax mode=abs outpar=nepmg
+suwind key=gx min=0 max=0 itmax=511 < referenceP_rp.su | sumax mode=abs outpar=neprf
+mg=`cat nepmg | awk '{print $1}'`
+rf=`cat neprf | awk '{print $1}'`
+value=${value/[eE][+][0]/*10^}
+mg=${mg/[eE][+][0]/*10^}
+rf=${rf/[eE][+][0]/*10^}
+rm nep*
+scale=$(echo "scale=3; ($rf)/($mg)" | bc -l)
+#echo $scale
+
+(suwind key=gx min=0 max=0 itmax=511 < pgreen_$piter.su | sugain scale=$scale; \
+ suwind key=gx min=0 max=0 < referenceP_rp.su) | \
+ supsgraph hbox=6 wbox=1 labelsize=10 linegray=0.3,0.0 style=seismic \
+ lineon=3.5,1.0 lineoff=1.3,0.0 > compare_$piter.eps
+
+done
+
diff --git a/marchenko/demo/oneD/epsModel.scr b/marchenko/demo/oneD/epsModel.scr
new file mode 100755
index 0000000000000000000000000000000000000000..6a31b126a03c9db8ec78547b1fe2abff5d9a419d
--- /dev/null
+++ b/marchenko/demo/oneD/epsModel.scr
@@ -0,0 +1,64 @@
+#!/bin/bash
+
+export PATH=$HOME/src/OpenSource/bin/:$PATH:
+
+#Postscript file of model and shot record at xsrc=0
+cat << EOF1 > line1
+400 -2500
+400 2500
+EOF1
+
+cat << EOF2 > line2
+700 -2500
+700 2500
+EOF2
+
+cat << EOF3 > line3
+1100 -2500
+1100 2500
+EOF3
+
+#model
+supsimage hbox=4 wbox=6 labelsize=12 < model10_cp.su \
+ x1beg=0 x1end=1400.0 d1num=200 legend=1 threecolor=0 \
+ curve=line1,line2,line3 npair=2,2,2 curvecolor=black,black,black \
+ n1tic=5 x2beg=-2250 f2num=-2000 d2num=1000 x2end=2250 > model_cp_line.eps
+
+supsimage hbox=4 wbox=6 labelsize=12 < model10_ro.su \
+ x1beg=0 x1end=1400.0 d1num=200 legend=1 threecolor=0 \
+ curve=line1,line2,line3 npair=2,2,2 curvecolor=black,black,black \
+ n1tic=5 x2beg=-2250 f2num=-2000 d2num=1000 x2end=2250 > model_ro_line.eps
+
+#wavelet
+dt=0.0005
+supsgraph < wavefw.su \
+ labelsize=12 d1=$dt style=normal \
+ label1="time (s)" label2="amplitude" \
+ d1num=0.15 wbox=6 hbox=3 x1end=0.9 > wavefw.eps
+
+sufft < wavefw.su | suamp | sugain scale=$dt | supsgraph \
+ labelsize=12 style=normal \
+ label1="time (s)" label2="amplitude" \
+ d1num=10 wbox=6 hbox=3 x1end=125 x2end=1.1 > wavefw_freq.eps
+
+
+#shot record
+file=shot5_rp.su
+sumax < $file mode=abs outpar=nep
+clipref=`cat nep | awk '{print $1/3}'`
+suwind key=gx min=-2250000 max=2250000 < $file | \
+ supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 \
+ n1tic=2 d2=5 f1=0.0 x1beg=0 x1end=2.004 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > shotx0_rp.eps
+
+#Initial focusing operator
+file=iniFocus_rp.su
+file_base=${file%.su}
+sumax < $file mode=abs outpar=nep
+clipref=`cat nep | awk '{print $1/3}'`
+suwind key=gx min=-2250000 max=2250000 < $file | \
+ supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 \
+ n1tic=2 d2=5 f1=0.0 x1beg=0 x1end=2.004 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
+
+rm nep
diff --git a/marchenko/demo/oneD/initialFocus.scr b/marchenko/demo/oneD/initialFocus.scr
new file mode 100755
index 0000000000000000000000000000000000000000..1bbacbf480f4fc599fa49b13c42543a6671179db
--- /dev/null
+++ b/marchenko/demo/oneD/initialFocus.scr
@@ -0,0 +1,39 @@
+#!/bin/bash
+
+export PATH=$HOME/src/OpenSource/bin:$PATH:
+
+dx=2.5
+dt=0.0005
+
+#the model upto 900 m depth, deeper reflections are not needed to model the direct transmission response
+makemod sizex=10000 sizez=1400 dx=$dx dz=$dx cp0=1800 ro0=1000 \
+ orig=-5000,0 file_base=modelup.su verbose=2 \
+ intt=def x=-5000,5000 z=400,400 poly=0 cp=2300 ro=3000 \
+ intt=def x=-5000,5000 z=700,700 poly=0 cp=2000 ro=1100
+
+#makewave fp=25 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
+
+export OMP_NUM_THREADS=1
+
+fdelmodc \
+ file_cp=modelup_cp.su ischeme=1 iorder=4 \
+ file_den=modelup_ro.su \
+ file_src=wave.su \
+ file_rcv=iniFocus.su \
+ src_type=1 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.004 \
+ rec_delay=0.1 \
+ verbose=2 \
+ tmod=2.144 \
+ dxrcv=5 \
+ xrcv1=-2250 xrcv2=2250 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=0 zsrc=900 \
+ ntaper=101 \
+ left=2 right=2 top=2 bottom=2
+
diff --git a/marchenko/demo/oneD/marchenko.scr b/marchenko/demo/oneD/marchenko.scr
new file mode 100755
index 0000000000000000000000000000000000000000..16502d5cfb9fb64cf5154af5cfd51ca9c27d64f8
--- /dev/null
+++ b/marchenko/demo/oneD/marchenko.scr
@@ -0,0 +1,32 @@
+#!/bin/bash -x
+
+export PATH=$HOME/src/OpenSource/bin:$PATH:
+export OMP_NUM_THREADS=1
+
+#mute all events below the first arrival to get the intial focusing field
+fmute file_shot=iniFocus_rp.su file_out=p0plus.su above=-1 shift=-8 verbose=1 check=0 hw=4
+
+#apply the Marchenko algorithm
+marchenko file_shot=shotsdx5_rp.su file_tinv=p0plus.su nshots=901 verbose=1 \
+ tap=0 ntap=41 niter=8 hw=12 shift=8 smooth=5 \
+ file_green=pgreen.su file_gplus=Gplus0.su file_gmin=Gmin0.su \
+ file_f1plus=f1plus0.su file_f1min=f1min0.su file_f2=f2.su
+
+#compare Green's funtions on Marhcenko and reference result
+suwind key=gx min=0 max=0 itmax=511 < pgreen.su | sumax mode=abs outpar=nepmg
+suwind key=gx min=0 max=0 itmax=511 < referenceP_rp.su | sumax mode=abs outpar=neprf
+mg=`cat nepmg | awk '{print $1}'`
+rf=`cat neprf | awk '{print $1}'`
+value=${value/[eE][+][0]/*10^}
+mg=${mg/[eE][+][0]/*10^}
+rf=${rf/[eE][+][0]/*10^}
+rm nep*
+scale=$(echo "scale=3; ($rf)/($mg)" | bc -l)
+echo $scale
+
+(suwind key=gx min=0 max=0 itmax=511 < pgreen.su | sugain scale=$scale; \
+ suwind key=gx min=0 max=0 < referenceP_rp.su) | suxgraph
+
+suwind itmax=511 < pgreen.su > pgreen512.su
+suop2 pgreen512.su referenceP_rp.su op=diff w2=1 w1=$scale > diffref.su
+
diff --git a/marchenko/demo/oneD/marchenkoIter.scr b/marchenko/demo/oneD/marchenkoIter.scr
new file mode 100755
index 0000000000000000000000000000000000000000..401f97f7c2108e92e0ff5ca813d9fdfd2b4d183a
--- /dev/null
+++ b/marchenko/demo/oneD/marchenkoIter.scr
@@ -0,0 +1,21 @@
+#!/bin/bash -x
+
+export PATH=$HOME/src/OpenSource/bin:$PATH:
+export OMP_NUM_THREADS=1
+
+#mute all events below the first arrival to get the intial focusing field
+fmute file_shot=iniFocus_rp.su file_out=p0plus.su above=-1 shift=-8 verbose=1 check=0 hw=4
+
+for (( iter=1; iter<=4; iter+=1 ))
+do
+echo "doing iteration $iter"
+piter=$(printf %03d $iter)
+
+#apply the Marchenko algorithm
+marchenko file_shot=shotsdx5_rp.su file_tinv=p0plus.su nshots=901 verbose=1 \
+ tap=0 ntap=41 niter=$iter hw=12 shift=8 smooth=5 \
+ file_green=pgreen_$piter.su file_iter=iter.su \
+ file_f1plus=f1plus_$piter.su file_f1min=f1min_$piter.su
+
+done
+
diff --git a/marchenko/demo/oneD/model.scr b/marchenko/demo/oneD/model.scr
new file mode 100755
index 0000000000000000000000000000000000000000..58c8e1beb99af636393cc203489901a1f40a0c1a
--- /dev/null
+++ b/marchenko/demo/oneD/model.scr
@@ -0,0 +1,77 @@
+#!/bin/bash
+
+#adjust this PATH to where the code is installed
+export PATH=$HOME/src/OpenSource/bin:$PATH:
+
+dx=2.5
+dt=0.0005
+
+#define gridded model for FD computations
+makemod sizex=10000 sizez=1400 dx=$dx dz=$dx cp0=1800 ro0=1000 \
+ orig=-5000,0 file_base=model10.su verbose=2 \
+ intt=def x=-5000,5000 z=400,400 poly=0 cp=2300 ro=3000 \
+ intt=def x=-5000,5000 z=700,700 poly=0 cp=2000 ro=1100 \
+ intt=def x=-5000,5000 z=1100,1100 poly=0 cp=2500 ro=4000
+
+#define wavelet for modeling R
+makewave w=fw fmin=0 flef=5 frig=80 fmax=100 dt=$dt file_out=wavefw.su nt=4096 t0=0.3 scale=0 scfft=1
+
+#define wavelet for reference and intial focusing field.
+makewave fp=25 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
+
+export OMP_NUM_THREADS=4
+
+#Model shot record in middle of model
+fdelmodc \
+ file_cp=model10_cp.su ischeme=1 iorder=4 \
+ file_den=model10_ro.su \
+ file_src=wavefw.su \
+ file_rcv=shot5_fd.su \
+ src_type=7 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.004 \
+ rec_delay=0.3 \
+ verbose=2 \
+ tmod=4.392 \
+ dxrcv=5.0 \
+ xrcv1=-4500 xrcv2=4500 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=0 zsrc=0 \
+ ntaper=101 \
+ left=2 right=2 top=2 bottom=2
+
+#define homogenoeus model to compute direct wave only
+makemod sizex=10000 sizez=1200 dx=$dx dz=$dx cp0=1800 ro0=1000 \
+ orig=-5000,0 file_base=hom.su verbose=2
+
+#Model direct wave only in middle of model
+fdelmodc \
+ file_cp=hom_cp.su ischeme=1 iorder=4 \
+ file_den=hom_ro.su \
+ file_src=wavefw.su \
+ file_rcv=shot5_hom_fd.su \
+ src_type=7 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.004 \
+ rec_delay=0.3 \
+ verbose=2 \
+ tmod=4.392 \
+ dxrcv=5.0 \
+ xrcv1=-4500 xrcv2=4500 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=0 zsrc=0 \
+ ntaper=101 \
+ left=2 right=2 top=2 bottom=2
+
+#subtract direct wave from full model shot record: this defines R
+sudiff shot5_fd_rp.su shot5_hom_fd_rp.su > shot5_rp.su
+
+
diff --git a/marchenko/demo/oneD/p5all.scr b/marchenko/demo/oneD/p5all.scr
new file mode 100755
index 0000000000000000000000000000000000000000..a8f2115734fa167fa664290d1afcff785e95df02
--- /dev/null
+++ b/marchenko/demo/oneD/p5all.scr
@@ -0,0 +1,30 @@
+#!/bin/bash -x
+
+export PATH=$HOME/src/OpenSource/bin:$PATH:
+
+# Generate the full R matrix for a fixed spread geometry.
+
+dxshot=5000 # with scalco factor of 1000
+ishot=0
+nshots=901
+
+echo $1
+
+rm shotsdx5_rp.su
+
+while (( ishot < nshots ))
+do
+
+ (( xsrc = -2250000 + ${ishot}*${dxshot} ))
+ (( tr1 = 901 - ${ishot} ))
+ (( tr2 = ${tr1} + 900 ))
+ echo xsrc=$xsrc tr1=$tr1 tr2=$tr2
+
+ (( ishot = $ishot + 1))
+
+ suwind < shot5_rp.su key=tracl min=$tr1 max=$tr2 | \
+ sushw key=sx,gx,fldr,trwf \
+ a=$xsrc,-2250000,$ishot,901 b=0,5000,0,0 j=0,901,0,0 >> shotsdx5_rp.su
+
+done
+
diff --git a/marchenko/demo/oneD/referenceShot.scr b/marchenko/demo/oneD/referenceShot.scr
new file mode 100755
index 0000000000000000000000000000000000000000..b7a2b771341b3115d71bfebe2ec06e308846cbc6
--- /dev/null
+++ b/marchenko/demo/oneD/referenceShot.scr
@@ -0,0 +1,34 @@
+#!/bin/bash
+
+export PATH=$HOME/src/OpenSource/bin:$PATH:
+
+#Compute the reference Green's fucntion at x=0 z=900 m in the actual model
+dx=2.5
+dt=0.0005
+
+makewave fp=25 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
+
+export OMP_NUM_THREADS=2
+
+fdelmodc \
+ file_cp=model10_cp.su ischeme=1 iorder=4 \
+ file_den=model10_ro.su \
+ file_src=wave.su \
+ file_rcv=referenceP.su \
+ src_type=1 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.004 \
+ rec_delay=0.1 \
+ verbose=2 \
+ tmod=2.144 \
+ dxrcv=5.0 \
+ xrcv1=-2250 xrcv2=2250 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=0 zsrc=900 \
+ ntaper=101 \
+ left=2 right=2 top=2 bottom=2
+
diff --git a/marchenko/demo/twoD/README b/marchenko/demo/twoD/README
new file mode 100644
index 0000000000000000000000000000000000000000..5fc50362f6ee58da1104668356b8ac28380d6013
--- /dev/null
+++ b/marchenko/demo/twoD/README
@@ -0,0 +1,9 @@
+Description of files:
+1) shots.scr create the shots
+2) model.scr computes the model
+3) direct_wave.scr crate the direct wave to be removed from the shots
+4) remove_direct.scr remove the direct wave from the shots and scale them
+5) first_arrival.scr computes the first arrival
+6) marchenko.scr perform the Marchenko scheme
+7) referenceShot.scr creates the reference Green's function
+
diff --git a/marchenko/demo/twoD/clean b/marchenko/demo/twoD/clean
new file mode 100755
index 0000000000000000000000000000000000000000..085f0557f80ea4ee096c73abfa9b42a078e823bb
--- /dev/null
+++ b/marchenko/demo/twoD/clean
@@ -0,0 +1,4 @@
+#!/bin/bash
+
+rm *.su *.bin *.txt *.eps nep
+
diff --git a/marchenko/demo/direct.scr b/marchenko/demo/twoD/direct.scr
similarity index 100%
rename from marchenko/demo/direct.scr
rename to marchenko/demo/twoD/direct.scr
diff --git a/marchenko/demo/first_arrival.scr b/marchenko/demo/twoD/first_arrival.scr
similarity index 100%
rename from marchenko/demo/first_arrival.scr
rename to marchenko/demo/twoD/first_arrival.scr
diff --git a/marchenko/demo/marchenko.scr b/marchenko/demo/twoD/marchenko.scr
similarity index 100%
rename from marchenko/demo/marchenko.scr
rename to marchenko/demo/twoD/marchenko.scr
diff --git a/marchenko/demo/model.scr b/marchenko/demo/twoD/model.scr
similarity index 100%
rename from marchenko/demo/model.scr
rename to marchenko/demo/twoD/model.scr
diff --git a/marchenko/demo/referenceShot.scr b/marchenko/demo/twoD/referenceShot.scr
similarity index 100%
rename from marchenko/demo/referenceShot.scr
rename to marchenko/demo/twoD/referenceShot.scr
diff --git a/marchenko/demo/remove_direct.scr b/marchenko/demo/twoD/remove_direct.scr
similarity index 100%
rename from marchenko/demo/remove_direct.scr
rename to marchenko/demo/twoD/remove_direct.scr
diff --git a/marchenko/demo/shots.scr b/marchenko/demo/twoD/shots.scr
similarity index 100%
rename from marchenko/demo/shots.scr
rename to marchenko/demo/twoD/shots.scr
diff --git a/marchenko/marchenko.c b/marchenko/marchenko.c
index 0106ea32a89dcd2d2ac361afd7ab3d952bc6179f..535085b4d04a06c7cfa25226851a67f8720ce207 100644
--- a/marchenko/marchenko.c
+++ b/marchenko/marchenko.c
@@ -38,7 +38,7 @@ int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
int disp_fileinfo(char *file, int n1, int n2, float f1, float f2, float d1, float d2, segy *hdrs);
double wallclock_time(void);
-void synthesis(complex *Refl, complex *Fop, float *syndata, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int reci, int nshots, int verbose);
+void synthesis(complex *Refl, complex *Fop, float *iRN, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int reci, int nshots, int verbose);
void synthesisPosistions(int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int reci, int nshots, int *ixpossyn, int *npossyn, int verbose);
@@ -56,12 +56,9 @@ char *sdoc[] = {
" ",
" Optional parameters: ",
" ",
-" SYNTHESIS ",
-" ixa=0 .................... number of traces after focus point",
-" ixb=ixa .................. number of traces before focus point",
+" INTEGRATION ",
" tap=0 .................... lateral taper focusing(1), shot(2) or both(3)",
" ntap=0 ................... number of taper points at boundaries",
-" reci=0 ................... 1; add focusing in emission 2; emission only",
" fmin=0 ................... minimum frequency",
" fmax=70 .................. maximum frequency",
" MARCHENKO ITERATIONS ",
@@ -84,9 +81,8 @@ char *sdoc[] = {
" file_iter= ............... output file with N for each iteration",
" verbose=0 ................ silent option; >0 displays info",
" ",
-" Note that if ixa=0 and ixb=0 all shots are used.",
" ",
-" author : Jan Thorbecke : 2013 (j.w.thorbecke@tudelft.nl)",
+" author : Jan Thorbecke : 2016 (j.w.thorbecke@tudelft.nl)",
" ",
NULL};
/**************** end self doc ***********************************/
@@ -105,7 +101,7 @@ int main (int argc, char **argv)
double t0, t1, t2, t3, tsyn, tread, tfft;
float *shotdata, d1, d2, f1, f2, fts, fxs, ft, fx, *xsyn, dxsrc;
float *green, *pplus, *f2p, *pmin, *tinv, *mute, dt, dx, dts, dxs, scl, mem;
- float *f1plus, *f1min, *Nk, *Nk_1, *trace, *Gmin, *Gplus;
+ float *f1plus, *f1min, *iRN, *Ni, *trace, *Gmin, *Gplus;
float max, scel, xmin, xmax, weight;
complex *Refl, *Fop, *ctrace;
char *file_tinv, *file_shot, *file_green, *file_iter;
@@ -153,7 +149,7 @@ int main (int argc, char **argv)
if (reci && ntap) vwarn("tapering influences the reciprocal result");
-/*================ Reading info about shot and focusing operator sizes ================*/
+/*================ Reading info about shot and initial operator sizes ================*/
ngath = 0; /* setting ngath=0 scans all traces; n2 contains maximum traces/gather */
ret = getFileInfo(file_tinv, &n1, &n2, &ngath, &d1, &d2, &f1, &f2, &xmin, &xmax, &scl, &ntraces);
@@ -171,7 +167,6 @@ int main (int argc, char **argv)
if (!getparfloat("dt", &dt)) dt = d1;
ntfft = optncr(MAX(nt, nts));
-// nts = ntfft;
nf = ntfft/2+1;
df = 1.0/(ntfft*dt);
nfreq = ntfft/2+1;
@@ -196,8 +191,8 @@ int main (int argc, char **argv)
Gplus = (float *)calloc(Nsyn*nxs*ntfft,sizeof(float));
f1plus = (float *)calloc(Nsyn*nxs*ntfft,sizeof(float));
f1min = (float *)calloc(Nsyn*nxs*ntfft,sizeof(float));
- Nk = (float *)calloc(Nsyn*nxs*ntfft,sizeof(float));
- Nk_1 = (float *)calloc(Nsyn*nxs*ntfft,sizeof(float));
+ iRN = (float *)calloc(Nsyn*nxs*ntfft,sizeof(float));
+ Ni = (float *)calloc(Nsyn*nxs*ntfft,sizeof(float));
ctrace = (complex *)malloc(ntfft*sizeof(complex));
trace = (float *)malloc(ntfft*sizeof(float));
mute = (float *)calloc(Nsyn*nxs,sizeof(float));
@@ -217,9 +212,10 @@ int main (int argc, char **argv)
mode=-1; /* apply complex conjugate to read in data */
readTinvData(file_tinv, xrcvsyn, xsyn, zsyn, xnxsyn, Fop, nw, nw_low, Nsyn, nxs, ntfft,
mode, mute, tinv, hw, verbose);
-/* reading data added zero's to the number of time samples to be the same as ntfft */
+ /* reading data added zero's to the number of time samples to be the same as ntfft */
nts = ntfft;
+ /* define tapers to taper edges of acquisition */
if (tap == 1 || tap == 3) {
for (j = 0; j < ntap; j++)
tapersy[j] = (cos(PI*(j-ntap)/ntap)+1)/2.0;
@@ -243,6 +239,7 @@ int main (int argc, char **argv)
}
}
+ /* check consistency of header values */
if (xrcvsyn[0] != 0 || xrcvsyn[1] != 0 ) fxs = xrcvsyn[0];
fxs2 = fxs + (float)(nxs-1)*dxs;
dxf = (xrcvsyn[nxs-1] - xrcvsyn[0])/(float)(nxs-1);
@@ -283,6 +280,7 @@ int main (int argc, char **argv)
}
free(tapersh);
+ /* check consistency of header values */
fxf = xsrc[0];
if (nx > 1) dxf = (xrcv[0] - xrcv[nx-1])/(float)(nx-1);
else dxf = d2;
@@ -389,9 +387,9 @@ int main (int argc, char **argv)
t2 = wallclock_time();
-/*================ construction of Nk(-t) = - \int R(x,t) Fop(t) ================*/
+/*================ construction of Ni(-t) = - \int R(x,t) Fop(t) ================*/
- synthesis(Refl, Fop, Nk, nx, nt, nxs, nts, dt, xsyn, Nsyn,
+ synthesis(Refl, Fop, iRN, nx, nt, nxs, nts, dt, xsyn, Nsyn,
xrcv, xsrc, fxs2, fxs, dxs, dxsrc, dx, ixa, ixb, ntfft, nw, nw_low, nw_high,
reci, nshots, verbose);
@@ -399,7 +397,7 @@ int main (int argc, char **argv)
memset(&Fop[0].r, 0, Nsyn*nxs*nw*2*sizeof(float));
if (file_iter != NULL) {
- writeDataIter(file_iter, Nk, hdrs_out, ntfft, nxs, d2, f2, n2out, Nsyn, xsyn, zsyn, iter);
+ writeDataIter(file_iter, iRN, hdrs_out, ntfft, nxs, d2, f2, n2out, Nsyn, xsyn, zsyn, iter);
}
/* initialization */
@@ -408,32 +406,32 @@ int main (int argc, char **argv)
for (l = 0; l < Nsyn; l++) {
for (i = 0; i < npossyn; i++) {
j = 0;
- Nk_1[l*nxs*nts+i*nts+j] = -Nk[l*nxs*nts+i*nts+j];
+ Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
- Nk_1[l*nxs*nts+i*nts+j] = -Nk[l*nxs*nts+i*nts+nts-j];
+ Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+i*nts+nts-j];
}
}
}
- /* p0^-(x,t) = Nk = (R * T_d^inv)(t) */
+ /* p0^-(x,t) = iRN = (R * T_d^inv)(t) */
for (l = 0; l < Nsyn; l++) {
for (i = 0; i < npossyn; i++) {
j = 0;
- pmin[l*nxs*nts+i*nts+j] = Nk[l*nxs*nts+i*nts+j];
+ pmin[l*nxs*nts+i*nts+j] = iRN[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
- pmin[l*nxs*nts+i*nts+j] = Nk[l*nxs*nts+i*nts+j];
+ pmin[l*nxs*nts+i*nts+j] = iRN[l*nxs*nts+i*nts+j];
}
}
}
- applyMute(Nk_1, mute, smooth, above, Nsyn, nxs, nts, xsrc, ixpossyn, npossyn, shift);
+ applyMute(Ni, mute, smooth, above, Nsyn, nxs, nts, xsrc, ixpossyn, npossyn, shift);
- /* even iterations: => - f_1^-(-t) = windowed(Nk) */
+ /* even iterations: => - f_1^-(-t) = windowed(iRN) */
for (l = 0; l < Nsyn; l++) {
for (i = 0; i < npossyn; i++) {
j = 0;
- f1min[l*nxs*nts+i*nts+j] -= Nk_1[l*nxs*nts+i*nts+j];
+ f1min[l*nxs*nts+i*nts+j] -= Ni[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
- f1min[l*nxs*nts+i*nts+j] -= Nk_1[l*nxs*nts+i*nts+nts-j];
+ f1min[l*nxs*nts+i*nts+j] -= Ni[l*nxs*nts+i*nts+nts-j];
}
}
}
@@ -444,9 +442,9 @@ int main (int argc, char **argv)
//ix = NINT((xsrc[i]-fxs)/dxs);
ix = ixpossyn[i];
//fprintf(stderr,"i=%d xsrc=%f ix=%d ixpossyn=%d\n", i, xsrc[i], ix, ixpossyn[i]);
- f2p[l*nxs*nts+i*nts+j] = tinv[l*nxs*nts+ix*nts+j] + Nk_1[l*nxs*nts+i*nts+j];
+ f2p[l*nxs*nts+i*nts+j] = tinv[l*nxs*nts+ix*nts+j] + Ni[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
- f2p[l*nxs*nts+i*nts+j] = tinv[l*nxs*nts+ix*nts+j] + Nk_1[l*nxs*nts+i*nts+j];
+ f2p[l*nxs*nts+i*nts+j] = tinv[l*nxs*nts+ix*nts+j] + Ni[l*nxs*nts+i*nts+j];
}
}
}
@@ -463,26 +461,26 @@ int main (int argc, char **argv)
}
}
else if (iter==1) {
- /* Nk_1(x,t) = -\int R(x,t) M_0(x,-t) dxdt*/
+ /* Ni(x,t) = -\int R(x,t) M_0(x,-t) dxdt*/
for (l = 0; l < Nsyn; l++) {
for (i = 0; i < npossyn; i++) {
j = 0;
- Nk_1[l*nxs*nts+i*nts+j] = -Nk[l*nxs*nts+i*nts+j];
+ Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
- Nk_1[l*nxs*nts+i*nts+j] = -Nk[l*nxs*nts+i*nts+nts-j];
+ Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+i*nts+nts-j];
}
}
}
for (l = 0; l < Nsyn; l++) {
for (i = 0; i < npossyn; i++) {
j = 0;
- pmin[l*nxs*nts+i*nts+j] -= Nk_1[l*nxs*nts+i*nts+j];
+ pmin[l*nxs*nts+i*nts+j] -= Ni[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
- pmin[l*nxs*nts+i*nts+j] -= Nk_1[l*nxs*nts+i*nts+nts-j];
+ pmin[l*nxs*nts+i*nts+j] -= Ni[l*nxs*nts+i*nts+nts-j];
}
}
}
- applyMute(Nk_1, mute, smooth, above, Nsyn, nxs, nts, xsrc, ixpossyn, npossyn, shift);
+ applyMute(Ni, mute, smooth, above, Nsyn, nxs, nts, xsrc, ixpossyn, npossyn, shift);
/* Pressure based scheme */
for (l = 0; l < Nsyn; l++) {
for (i = 0; i < npossyn; i++) {
@@ -496,9 +494,9 @@ int main (int argc, char **argv)
for (l = 0; l < Nsyn; l++) {
for (i = 0; i < npossyn; i++) {
j = 0;
- f2p[l*nxs*nts+i*nts+j] += Nk_1[l*nxs*nts+i*nts+j];
+ f2p[l*nxs*nts+i*nts+j] += Ni[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
- f2p[l*nxs*nts+i*nts+j] += Nk_1[l*nxs*nts+i*nts+j];
+ f2p[l*nxs*nts+i*nts+j] += Ni[l*nxs*nts+i*nts+j];
}
}
}
@@ -507,9 +505,9 @@ int main (int argc, char **argv)
for (i = 0; i < npossyn; i++) {
j = 0;
ix = ixpossyn[i];
- f1plus[l*nxs*nts+i*nts+j] = tinv[l*nxs*nts+ix*nts+j] + Nk_1[l*nxs*nts+i*nts+j];
+ f1plus[l*nxs*nts+i*nts+j] = tinv[l*nxs*nts+ix*nts+j] + Ni[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
- f1plus[l*nxs*nts+i*nts+j] = tinv[l*nxs*nts+ix*nts+j] + Nk_1[l*nxs*nts+i*nts+j];
+ f1plus[l*nxs*nts+i*nts+j] = tinv[l*nxs*nts+ix*nts+j] + Ni[l*nxs*nts+i*nts+j];
}
}
}
@@ -520,22 +518,22 @@ int main (int argc, char **argv)
for (l = 0; l < Nsyn; l++) {
for (i = 0; i < npossyn; i++) {
j = 0;
- Nk_1[l*nxs*nts+i*nts+j] = -Nk[l*nxs*nts+i*nts+j];
+ Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
- Nk_1[l*nxs*nts+i*nts+j] = -Nk[l*nxs*nts+i*nts+nts-j];
+ Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+i*nts+nts-j];
}
}
}
for (l = 0; l < Nsyn; l++) {
for (i = 0; i < npossyn; i++) {
j = 0;
- pmin[l*nxs*nts+i*nts+j] -= Nk_1[l*nxs*nts+i*nts+j];
+ pmin[l*nxs*nts+i*nts+j] -= Ni[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
- pmin[l*nxs*nts+i*nts+j] -= Nk_1[l*nxs*nts+i*nts+nts-j];
+ pmin[l*nxs*nts+i*nts+j] -= Ni[l*nxs*nts+i*nts+nts-j];
}
}
}
- applyMute(Nk_1, mute, smooth, above, Nsyn, nxs, nts, xsrc, ixpossyn, npossyn, shift);
+ applyMute(Ni, mute, smooth, above, Nsyn, nxs, nts, xsrc, ixpossyn, npossyn, shift);
/* compute full Green's function G = p^+(-t) + p^-(t) */
if (iter == niter-1) {
@@ -554,9 +552,9 @@ int main (int argc, char **argv)
for (l = 0; l < Nsyn; l++) {
for (i = 0; i < npossyn; i++) {
j = 0;
- f2p[l*nxs*nts+i*nts+j] += Nk_1[l*nxs*nts+i*nts+j];
+ f2p[l*nxs*nts+i*nts+j] += Ni[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
- f2p[l*nxs*nts+i*nts+j] += Nk_1[l*nxs*nts+i*nts+j];
+ f2p[l*nxs*nts+i*nts+j] += Ni[l*nxs*nts+i*nts+j];
}
}
}
@@ -566,9 +564,9 @@ int main (int argc, char **argv)
for (l = 0; l < Nsyn; l++) {
for (i = 0; i < npossyn; i++) {
j = 0;
- f1min[l*nxs*nts+i*nts+j] -= Nk_1[l*nxs*nts+i*nts+j];
+ f1min[l*nxs*nts+i*nts+j] -= Ni[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
- f1min[l*nxs*nts+i*nts+j] -= Nk_1[l*nxs*nts+i*nts+nts-j];
+ f1min[l*nxs*nts+i*nts+j] -= Ni[l*nxs*nts+i*nts+nts-j];
}
}
}
@@ -577,9 +575,9 @@ int main (int argc, char **argv)
for (l = 0; l < Nsyn; l++) {
for (i = 0; i < npossyn; i++) {
j = 0;
- f1plus[l*nxs*nts+i*nts+j] += Nk_1[l*nxs*nts+i*nts+j];
+ f1plus[l*nxs*nts+i*nts+j] += Ni[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
- f1plus[l*nxs*nts+i*nts+j] += Nk_1[l*nxs*nts+i*nts+j];
+ f1plus[l*nxs*nts+i*nts+j] += Ni[l*nxs*nts+i*nts+j];
}
}
}
@@ -609,7 +607,7 @@ int main (int argc, char **argv)
}
}
- synthesis(Refl, Fop, Nk, nx, nt, nxs, nts, dt, xsyn, Nsyn,
+ synthesis(Refl, Fop, iRN, nx, nt, nxs, nts, dt, xsyn, Nsyn,
xrcv, xsrc, fxs2, fxs, dxs, dxsrc, dx, ixa, ixb, ntfft, nw, nw_low, nw_high,
reci, nshots, verbose);
@@ -617,9 +615,9 @@ int main (int argc, char **argv)
for (l = 0; l < Nsyn; l++) {
for (i = 0; i < npossyn; i++) {
j=0;
- Gmin[l*nxs*nts+i*nts+j] = Nk[l*nxs*nts+i*nts+j] - f1min[l*nxs*nts+i*nts+j];
+ Gmin[l*nxs*nts+i*nts+j] = iRN[l*nxs*nts+i*nts+j] - f1min[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
- Gmin[l*nxs*nts+i*nts+j] = Nk[l*nxs*nts+i*nts+j] - f1min[l*nxs*nts+i*nts+j];
+ Gmin[l*nxs*nts+i*nts+j] = iRN[l*nxs*nts+i*nts+j] - f1min[l*nxs*nts+i*nts+j];
}
}
}
@@ -641,7 +639,7 @@ int main (int argc, char **argv)
}
}
- synthesis(Refl, Fop, Nk, nx, nt, nxs, nts, dt, xsyn, Nsyn,
+ synthesis(Refl, Fop, iRN, nx, nt, nxs, nts, dt, xsyn, Nsyn,
xrcv, xsrc, fxs2, fxs, dxs, dxsrc, dx, ixa, ixb, ntfft, nw, nw_low, nw_high,
reci, nshots, verbose);
@@ -649,18 +647,18 @@ int main (int argc, char **argv)
for (l = 0; l < Nsyn; l++) {
for (i = 0; i < npossyn; i++) {
j=0;
- Gplus[l*nxs*nts+i*nts+j] = -Nk[l*nxs*nts+i*nts+j] + f1plus[l*nxs*nts+i*nts+j];
+ Gplus[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+i*nts+j] + f1plus[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
- Gplus[l*nxs*nts+i*nts+j] = -Nk[l*nxs*nts+i*nts+j] + f1plus[l*nxs*nts+i*nts+nts-j];
+ Gplus[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+i*nts+j] + f1plus[l*nxs*nts+i*nts+nts-j];
}
}
}
} /* end if for last iteration */
- /* transform muted Nk_1 to frequency domain */
+ /* transform muted Ni to frequency domain */
for (l = 0; l < Nsyn; l++) {
for (i = 0; i < npossyn; i++) {
- rc1fft(&Nk_1[l*nxs*nts+i*nts],ctrace,ntfft,-1);
+ rc1fft(&Ni[l*nxs*nts+i*nts],ctrace,ntfft,-1);
ix = ixpossyn[i];
for (iw=0; iw<nw; iw++) {
Fop[l*nxs*nw+iw*nxs+ix].r = ctrace[nw_low+iw].r;
@@ -812,11 +810,15 @@ int main (int argc, char **argv)
free(hdrs_out);
free(tapersy);
+ free(Ni);
exit(0);
}
-void synthesis(complex *Refl, complex *Fop, float *syndata, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int reci, int nshots, int verbose)
+
+/*================ Convolution and Integration ================*/
+
+void synthesis(complex *Refl, complex *Fop, float *iRN, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int reci, int nshots, int verbose)
{
int nfreq, size, iox, inx;
float scl;
@@ -836,7 +838,7 @@ void synthesis(complex *Refl, complex *Fop, float *syndata, int nx, int nt, int
t0 = wallclock_time();
/* reset output data to zero */
- memset(&syndata[0], 0, Nsyn*nxs*nts*sizeof(float));
+ memset(&iRN[0], 0, Nsyn*nxs*nts*sizeof(float));
for (k=0; k<nshots; k++) {
@@ -870,7 +872,7 @@ void synthesis(complex *Refl, complex *Fop, float *syndata, int nx, int nt, int
#endif
#pragma omp parallel default(none) \
- shared(syndata, dx, npe, nw, verbose) \
+ shared(iRN, dx, npe, nw, verbose) \
shared(Refl, Nsyn, reci, xrcv, xsrc, xsyn, fxs, nxs, dxs) \
shared(nx, ixa, ixb, dxsrc, iox, inx, k, nfreq, nw_low, nw_high) \
shared(Fop, size, nts, ntfft, scl, ixsrc, stderr) \
@@ -881,44 +883,10 @@ void synthesis(complex *Refl, complex *Fop, float *syndata, int nx, int nt, int
#pragma omp for
for (l = 0; l < Nsyn; l++) {
-/*
- if (ixa || ixb) {
- if (reci == 0) {
- x0 = xsyn[l]-ixb*dxsrc;
- x1 = xsyn[l]+ixa*dxsrc;
- if ((xsrc[k] < x0) || (xsrc[k] > x1)) continue;
- ix = NINT((xsrc[k]-x0)/dxsrc);
- dosrc = 1;
- }
- else if (reci == 1) {
- x0 = xsyn[l]-ixb*dxs;
- x1 = xsyn[l]+ixa*dxs;
- if (((xsrc[k] < x0) || (xsrc[k] > x1)) &&
- (xrcv[k*nx+0] < x0) && (xrcv[k*nx+nx-1] < x0)) continue;
- if (((xsrc[k] < x0) || (xsrc[k] > x1)) &&
- (xrcv[k*nx+0] > x1) && (xrcv[k*nx+nx-1] > x1)) continue;
- if ((xsrc[k] < x0) || (xsrc[k] > x1)) dosrc = 0;
- else dosrc = 1;
- ix = NINT((xsrc[k]-x0)/dxs);
- }
- else if (reci == 2) {
- if (NINT(dxsrc/dx)*dx != NINT(dxsrc)) dx = dxs;
- x0 = xsyn[l]-ixb*dx;
- x1 = xsyn[l]+ixa*dx;
- if ((xrcv[k*nx+0] < x0) && (xrcv[k*nx+nx-1] < x0)) continue;
- if ((xrcv[k*nx+0] > x1) && (xrcv[k*nx+nx-1] > x1)) continue;
- }
- }
- else {
-*/
ix = k;
x0 = fxs;
x1 = fxs+dxs*nxs;
dosrc = 1;
-// }
-// if (reci == 1 && dosrc) ix = NINT((xsrc[k]-x0)/dxs);
-
-// if (reci < 2 && dosrc) {
for (j = 0; j < nfreq; j++) sum[j].r = sum[j].i = 0.0;
for (j = nw_low, m = 0; j <= nw_high; j++, m++) {
for (i = iox; i < inx; i++) {
@@ -934,37 +902,9 @@ void synthesis(complex *Refl, complex *Fop, float *syndata, int nx, int nt, int
{
cr1fft(sum, rdata, ntfft, 1);
}
-// fprintf(stderr,"synthesis[%d] = %d ix=%d\n", k, ixsrc, ix);
/* dx = receiver distance */
for (j = 0; j < nts; j++)
- syndata[l*size+ix*nts+j] += rdata[j]*scl*dx;
-// }
-
-/*
- if (reci == 1 || reci == 2) {
- for (j = 0; j < nfreq; j++) sum[j].r = sum[j].i = 0.0;
- for (i = iox; i < inx; i++) {
- if ((xrcv[k*nx+i] < x0) || (xrcv[k*nx+i] > x1)) continue;
- if (reci == 1) ix = NINT((xrcv[k*nx+i]-x0)/dxs);
- else ix = NINT((xrcv[k*nx+i]-x0)/dx);
-
- for (j = nw_low, m = 0; j < nw_high; j++, m++) {
- tmp = Fop[l*nw*nxs+m*nxs+ixsrc];
- sum[j].r = Refl[k*nw*nx+m*nx+i].r*tmp.r -
- Refl[k*nw*nx+m*nx+i].i*tmp.i;
- sum[j].i = Refl[k*nw*nx+m*nx+i].i*tmp.r +
- Refl[k*nw*nx+m*nx+i].r*tmp.i;
- }
-#pragma omp critical
-{
- cr1fft(sum, rdata, ntfft, 1);
-}
- // dxsrc = source distance
- for (j = 0; j < nts; j++)
- syndata[l*size+ix*nts+j] += rdata[j]*scl*dxsrc;
- }
- }
-*/
+ iRN[l*size+ix*nts+j] += rdata[j]*scl*dx;
} /* end of parallel Nsyn loop */
@@ -979,7 +919,7 @@ void synthesis(complex *Refl, complex *Fop, float *syndata, int nx, int nt, int
}
} /* end of parallel region */
- if (verbose>3) vmess("*** Shot gather %d processed ***", k);
+ if (verbose>3) vmess("*** Shot gather %d processed ***", k);
} /* end of nshots (k) loop */
diff --git a/utils/Makefile b/utils/Makefile
index e898e9667499b00b41cc3ace5133cde2095dbd28..ba62e4963987177461d43e464264be23fb9b5370 100644
--- a/utils/Makefile
+++ b/utils/Makefile
@@ -6,7 +6,7 @@ LIBS += -L$L -lgenfft -lm $(LIBSM)
#OPTC += -openmp
#OPTC += -g -O0
-ALL: makemod makewave extendModel fconv correigen green basop syn2d mat2su
+ALL: makemod makewave extendModel fconv correigen green basop syn2d mat2su ftr1d
SRCM = \
makemod.c \
@@ -111,6 +111,15 @@ SRCA = mat2su.c \
docpkge.c \
getpars.c
+SRCT = ftr1d.c \
+ getFileInfo.c \
+ readData.c \
+ writeData.c \
+ wallclock_time.c \
+ verbosepkg.c \
+ atopkge.c \
+ docpkge.c \
+ getpars.c
OBJM = $(SRCM:%.c=%.o)
@@ -157,7 +166,12 @@ OBJA = $(SRCA:%.c=%.o)
mat2su: $(OBJA)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o mat2su $(OBJA) $(LIBS)
-install: makemod makewave extendModel fconv correigen green basop syn2d mat2su
+OBJT = $(SRCT:%.c=%.o)
+
+ftr1d: $(OBJT)
+ $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o ftr1d $(OBJT) $(LIBS)
+
+install: makemod makewave extendModel fconv correigen green basop syn2d mat2su ftr1d
cp makemod $B
cp makewave $B
cp extendModel $B
@@ -167,12 +181,13 @@ install: makemod makewave extendModel fconv correigen green basop syn2d mat2su
cp basop $B
cp syn2d $B
cp mat2su $B
+ cp ftr1d $B
clean:
- rm -f core $(OBJM) makemod $(OBJW) makewave $(OBJE) extendModel $(OBJF) fconv $(OBJG) $(OBJC) correigen green $(OBJB) basop $(OBJJ) syn2d $(OBJS) $(OBJA) mat2su $(OBJH)
+ rm -f core $(OBJM) makemod $(OBJW) makewave $(OBJE) extendModel $(OBJF) fconv $(OBJG) $(OBJC) correigen green $(OBJB) basop $(OBJJ) syn2d $(OBJS) $(OBJA) mat2su $(OBJH) ftr1d $(OBJT)
realclean: clean
- rm -f $B/makemod $B/makewave $B/extendModel $B/fconv $B/correigen $B/green $B/basop $B/syn2d $B/mat2su
+ rm -f $B/makemod $B/makewave $B/extendModel $B/fconv $B/correigen $B/green $B/basop $B/syn2d $B/mat2su $B/ftr1d
diff --git a/utils/ftr1d.c b/utils/ftr1d.c
new file mode 100644
index 0000000000000000000000000000000000000000..4406c889152d8b5005e9d600a2daba8513e9d752
--- /dev/null
+++ b/utils/ftr1d.c
@@ -0,0 +1,389 @@
+#include "par.h"
+#include "segy.h"
+#include <genfft.h>
+#include <time.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+ float r,i;
+} complex;
+#endif/* complex */
+
+int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *d2, float *f1, float *f2, float *xmin, float *xmax, float *sclsxgx, int *nxm);
+int readData(FILE *fp, float *data, segy *hdrs, int n1);
+int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
+int disp_fileinfo(char *file, int n1, int n2, float f1, float f2, float d1, float d2, segy *hdrs);
+double wallclock_time(void);
+
+/*********************** self documentation **********************/
+char *sdoc[] = {
+" ",
+" ftr1d - 1 Dimensional FFT",
+" ",
+" ftr1d file_in= [optional parameters]",
+" ",
+" Required parameters:",
+" ",
+" file_in= .................. input file (DELPHI format)",
+" ",
+" Optional parameters:",
+" ",
+" file_out= ................ Output file with the result (DELPHI format)",
+" n1=from input ............ number of samples to be Fourier transformed",
+" opt=default .............. kind of Fourier Transform (see below)",
+" scale=0 .................. scale for the FFT (0 gives defaults)",
+" sign=1/-1 ................ sign in the kernel of the FFT (see below)",
+" nxmax=512 ................ maximum number of traces in input file",
+" ntmax=1024 ............... maximum number of samples/trace in input file",
+" key=fldr ................. SU search key",
+" verbose=0 ................ silent option; >0 display info",
+" ",
+" Options for opt and defaults for scale and sign:",
+" - 1 = real -> complex scale = 1 sign=-1 (time-axis)",
+" - 2 = complex -> real scale = 1/n1 sign=1 (freq-axis)",
+" - 3 = complex -> complex scale = 1 sign=1 (x-axis)",
+" - 4 = complex -> complex scale = 1/n1 sign=-1 (kx-axis)",
+" ",
+" Which Fourier transform is carried out by default depends on the ",
+" axis of the array to be transformed. ",
+" If necessary zeros are added to the data. Note that if n1 is chosen bigger",
+" than the actual number of samples zeros are added to the data.",
+" ",
+" author : Jan Thorbecke : 09-08-1995 (J.W.Thorbecke@TUDelft.nl)",
+" ",
+NULL};
+/**************** end self doc ***********************************/
+
+void main(int argc, char *argv[])
+{
+ FILE *fp_in, *fp_out;
+ short trid, trid_out;
+ int optn, choice, sign, ntmax, nxmax, verbose, first, ngath, ntraces;
+ int nfreq, error, n1, n2, ret, n1_org, nf_org, i, j;
+ size_t size;
+ float *rdata, *datin, *datout, d1, d2, f1, f2, scale;
+ float scl, xmin, xmax;
+ double t0, t1, t2;
+ complex *cdata;
+ segy *hdrs;
+ char *file_in, *file_out;
+
+/* Reading input parameters */
+
+ t0 = wallclock_time();
+ initargs(argc, argv);
+ requestdoc(1);
+
+ if(!getparint("verbose", &verbose)) verbose=0;
+ if(!getparstring("file_in", &file_in)) {
+ if (verbose) vwarn("parameter file_in not found, assume pipe");
+ file_in = NULL;
+ }
+ if(!getparstring("file_out", &file_out)){
+ if (verbose) vwarn("parameter file_out not found, assume pipe");
+ file_out = NULL;
+ }
+ if(!getparint("n1", &optn)) optn = -1;
+ if(!getparint("opt", &choice)) choice = -1;
+ if(!getparfloat("scale", &scale)) scale = 0;
+ if(!getparint("sign", &sign)) sign = 0;
+ if(!getparint("ntmax", &ntmax)) ntmax = 1024;
+ if(!getparint("nxmax", &nxmax)) nxmax = 512;
+ n1 = 0;
+
+/* Opening input file for reading */
+ if (file_in != NULL) fp_in = fopen(file_in, "r");
+ else fp_in=stdin;
+ if (fp_in == NULL) verr("error on opening input file_in=%s", file_in);
+
+/* get dimensions */
+ ngath = 1;
+ error = getFileInfo(file_in, &n1, &n2, &ngath, &d1, &d2, &f1, &f2, &xmin, &xmax, &scl, &ntraces);
+ if (error == 0) {
+ if (!getparint("ntmax", &ntmax)) ntmax = n1;
+ if (!getparint("nxmax", &nxmax)) nxmax = n2;
+ if (verbose>=2 && (ntmax!=n1 || nxmax!=n2))
+ vmess("dimensions overruled: %d x %d",ntmax,nxmax);
+ }
+ else {
+ if (verbose>=2) vmess("dimensions used: %d x %d",ntmax,nxmax);
+ }
+ size = 2*ntmax * nxmax;
+// trid = hdrs[0].trid;
+// if (trid > FCMPLX) size *= 2;
+ datin = (float *) malloc(size*sizeof(float));
+ hdrs = (segy *) calloc(nxmax,sizeof(segy));
+ if (datin == NULL || hdrs==NULL ) {
+ verr("memory allocation error for input data");
+ }
+
+/* Opening output file for writing */
+ if (file_out==NULL) fp_out = stdout;
+ else {
+ fp_out = fopen(file_out, "w+");
+ if (fp_out==NULL) verr("error on creating output file");
+ }
+
+/* start reading data from input file */
+
+ error = 0;
+ first = 1;
+ while (error >= 0) {
+
+ n2 = readData(fp_in, datin, hdrs, n1);
+ if (n2 == 0) {
+ fclose(fp_in);
+ if (verbose) {
+ if (first) verr("error in reading first gather of file %s", file_in);
+ else {
+ vmess("end of data reached");
+ fclose(fp_out);
+ }
+ }
+ break;
+ }
+ n1 = hdrs[0].ns;
+ f1 = hdrs[0].f1;
+ f2 = hdrs[0].f2;
+ d1 = hdrs[0].d1;
+ d2 = hdrs[0].d2;
+ if (verbose) {
+ disp_fileinfo(file_in, n1, n2, f1, f2, d1, d2, hdrs);
+ }
+ trid = hdrs[0].trid;
+ if (!getparint("n1", &optn)) optn = -1;
+ if (optn < 0 ) optn = n1;
+ n1_org = n1;
+
+/* Determine from axis information which transform is desired */
+
+ if (choice == -1) {
+ if (trid == TREAL) choice = 1;
+ else if (trid == FUNPACKNYQ) choice = 2;
+ else if (trid == TCMPLX) choice = 3;
+ else if (trid == KOMEGA) choice = 4;
+ else vwarn("No sensible trid available use opt to define the kind of FFT");
+ if (verbose) vmess("Option %d for FFT is selected", choice);
+ }
+
+ t1 = wallclock_time();
+ if (first) if (verbose) vmess("CPU-time input = %.3f", t1-t0);
+
+/* The FFT */
+
+ if (choice == 1 ) {
+ optn = optncr(optn);
+ if(verbose) vmess("Transforming %d real points", optn);
+ nfreq = optn/2+1;
+ if (first) {
+ rdata = (float *) malloc(optn*n2*sizeof(float));
+ if(rdata == NULL) verr("memory allocation error rdata");
+ datout = (float *) malloc(2*nfreq*n2*sizeof(float));
+ if(datout == NULL) verr("memory allocation error datout");
+ }
+
+ if (sign == 0) sign = -1;
+ if (scale == 0) scale = 1.0;
+ for(i = 0; i < n2; i++) {
+ for(j = 0; j < n1; j++)
+ rdata[i*optn+j] = datin[j+n1*i]*scale;
+ for(j = n1; j < optn; j++)
+ rdata[i*optn+j] = 0.0;
+ }
+
+ rcmfft(&rdata[0], (complex *)&datout[0], optn, n2, optn, nfreq, sign);
+ trid_out=FUNPACKNYQ;
+ n1 = 2*nfreq;
+ d1 = 1.0/(optn*d1);
+ f1 = 0.0;
+
+ for(i = 0; i < n2; i++) {
+ hdrs[i].trid = trid_out;
+ hdrs[i].ns = n1;
+ hdrs[i].trwf = n2;
+ hdrs[i].nhs = n1_org;
+ hdrs[i].d1 = d1;
+ hdrs[i].f1 = f1;
+ }
+ }
+ else if (choice == 2) {
+ optn = optncr(optn-2);
+ nfreq = optn/2+1;
+ if(verbose) vmess("Transforming to %d real points", optn);
+ if (first) {
+ cdata = (complex *) malloc(nfreq*n2*sizeof(complex));
+ datout = (float *) malloc(optn*n2*sizeof(float));
+ }
+
+// TODO change n1 to nfreq of original samples
+ nf_org = n1_org/2;
+ if (scale == 0) scale = 1.0/optn;
+ for(i = 0; i < n2; i++) {
+ for(j = 0; j < nf_org; j++) {
+ cdata[i*nfreq+j].r = datin[2*j+2*nf_org*i]*scale;
+ cdata[i*nfreq+j].i = datin[2*j+2*nf_org*i+1]*scale;
+ }
+ for(j = nf_org; j < nfreq; j++) {
+ cdata[i*nfreq+j].r = 0.0;
+ cdata[i*nfreq+j].i = 0.0;
+ }
+ }
+
+ if (sign == 0 ) sign = 1;
+
+ crmfft(&cdata[0], &datout[0], optn, n2, nfreq, optn, sign);
+ trid_out = TREAL;
+ n1 = optn;
+ d1 = 1.0/(optn*d1);
+ f1 = 0.0;
+
+ /* copy data to original trace length */
+ if (hdrs[0].nhs>0 && hdrs[0].nhs<n1) {
+ n1_org=hdrs[0].nhs;
+ if (verbose)
+ vmess("reduce output trace length to original size %d",n1_org);
+ for (i = 0; i < n2; i++)
+ for (j = 0; j < n1_org; j++)
+ datout[i*n1_org+j]=datout[i*n1+j];
+ n1 = n1_org;
+ }
+
+ /* update headers for output file */
+ for(i = 0; i < n2; i++) {
+ hdrs[i].trid = trid_out;
+ hdrs[i].ns = n1;
+ hdrs[i].trwf = n2;
+ hdrs[i].d1 = d1;
+ hdrs[i].f1 = f1;
+ }
+ }
+ else if (choice == 3 ) {
+ optn = optncc(optn);
+ if(verbose) vmess("Transforming %d complex points", optn);
+ if (first) {
+ cdata = (complex *) malloc(optn*n2*sizeof(complex));
+ }
+
+ if (scale == 0) scale = 1.0;
+ if (trid == TCMPLX) {
+ for(i = 0; i < n2; i++) {
+ for(j = 0; j < n1; j++) {
+ cdata[i*optn+j].r = datin[2*j+2*n1*i]*scale;
+ cdata[i*optn+j].i = datin[2*j+2*n1*i+1]*scale;
+ }
+ for(j = n1; j < optn; j++) {
+ cdata[i*optn+j].r = 0.0;
+ cdata[i*optn+j].i = 0.0;
+ }
+ }
+ }
+ else {
+ vwarn("Real data is first made complex and then transformed");
+ for(i = 0; i < n2; i++) {
+ for(j = 0; j < n1; j++) {
+ cdata[i*optn+j].r = datin[j+n1*i]*scale;
+ cdata[i*optn+j].i = 0.0;
+ }
+ for(j = n1; j < optn; j++) {
+ cdata[i*optn+j].r = 0.0;
+ cdata[i*optn+j].i = 0.0;
+ }
+ }
+ }
+
+ if (sign == 0) sign = 1;
+ ccmfft(&cdata[0], optn, n2, optn, sign);
+ datout = (float *)&cdata[0].r;
+
+ trid_out = KOMEGA;
+ n1 = optn;
+ d1 = 2.0*PI/(optn*d1);
+ f1 = 0.0;
+
+ for(i = 0; i < n2; i++) {
+ hdrs[i].trid = trid_out;
+ hdrs[i].ns = n1;
+ hdrs[i].trwf = n2;
+ hdrs[i].d1 = d1;
+ hdrs[i].f1 = f1;
+ }
+ }
+ else if (choice == 4) {
+ optn = optncc(optn);
+ if(verbose) vmess("Transforming %d complex points", optn);
+ if (first) {
+ cdata = (complex *) malloc(optn*n2*sizeof(complex));
+ datout = (float *) malloc(2*optn*n2*sizeof(float));
+ }
+
+ if (scale == 0) scale = 1.0/optn;
+ if (trid == KOMEGA) {
+ for(i = 0; i < n2; i++) {
+ for(j = 0; j < n1; j++) {
+ cdata[i*optn+j].r = datin[2*j+2*n1*i]*scale;
+ cdata[i*optn+j].i = datin[2*j+2*n1*i+1]*scale;
+ }
+ for(j = n1; j < optn; j++) {
+ cdata[i*optn+j].r = 0.0;
+ cdata[i*optn+j].i = 0.0;
+ }
+ }
+ }
+ else {
+ vwarn("Real data is first made complex and then transformed");
+ for(i = 0; i < n2; i++) {
+ for(j = 0; j < n1; j++) {
+ cdata[i*optn+j].r = datin[j+n1*i]*scale;
+ cdata[i*optn+j].i = 0.0;
+ }
+ for(j = n1; j < optn; j++) {
+ cdata[i*optn+j].r = 0.0;
+ cdata[i*optn+j].i = 0.0;
+ }
+ }
+ }
+
+ if (sign == 0) sign = -1;
+ ccmfft(&cdata[0], optn, n2, optn, sign);
+ datout = (float *)&cdata[0].r;
+
+ trid_out = TCMPLX;
+ n1 = optn;
+ d1 = 2.0*PI/(optn*d1);
+ f1 = 0.0;
+
+ for(i = 0; i < n2; i++) {
+ hdrs[i].trid = trid_out;
+ hdrs[i].ns = n1;
+ hdrs[i].trwf = n2;
+ hdrs[i].d1 = d1;
+ hdrs[i].f1 = f1;
+ }
+ }
+ t2 = wallclock_time();
+ if (verbose)vmess("CPU-time FFT = %.3f", t2-t1);
+
+/* Write data to output */
+
+ ret = writeData(fp_out, datout, hdrs, n1, n2);
+ if (ret < 0 ) verr("error on writing output file.");
+
+ first = 0;
+ if (verbose) disp_fileinfo(file_out, n1, n2, f1, f2, d1, d2, hdrs);
+ }
+
+ free(hdrs);
+ free(datin);
+ if (datout !=NULL) free(datout);
+ if (rdata != NULL) free(rdata);
+ if (cdata != NULL) free(cdata);
+
+ t2 = wallclock_time();
+ if (verbose)vmess("Total CPU-time for ftr1d = %.3f", t2-t0);
+
+ exit(0);
+}
diff --git a/utils/getFileInfo.c b/utils/getFileInfo.c
index 9fb7710930214b601e6c714c0096e7595a764e0d..490ba4ad7e382117ca7612b6f28b039cc4f4b7f1 100644
--- a/utils/getFileInfo.c
+++ b/utils/getFileInfo.c
@@ -48,7 +48,7 @@ int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *
bytes = ftello( fp );
*n1 = hdr.ns;
- if (hdr.trid == 1 || hdr.dt != 0) {
+ if ( (hdr.trid == 1) && (hdr.dt != 0) ) {
*d1 = ((float) hdr.dt)*1.e-6;
*f1 = ((float) hdr.delrt)/1000.;
}