diff --git a/.gitignore b/.gitignore
index 64ba6df2394e1590a987ec1ecae0ab82345fc1b3..a0caee250f5bbd368528285a1171eea3b070906c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -27,3 +27,5 @@ marchenko_applications/fmute
marchenko_applications/marchenko
marchenko_full/fmute
marchenko_full/marchenko
+corrvir/corrvir
+fdemmodc/fdemmodc
diff --git a/corrvir/correlate.c b/corrvir/correlate.c
index 1a4ae8de7879ad6520fb56e116fb23346f14e2c7..e491ec1c8e879e240c19653ff38704ac29668cc6 100644
--- a/corrvir/correlate.c
+++ b/corrvir/correlate.c
@@ -64,7 +64,6 @@ int coherence(complex *cmaster, complex *cslaves, int nfreq, int ncor, float rep
den[j+icc] = am1*am2;
}
}
- leps = reps*maxden;
for (istation=0; istation<ncor; istation++) {
icc = istation*nfreq;
diff --git a/doc/fdelmodcManual.pdf b/doc/fdelmodcManual.pdf
index 3cebcbe08381ab0575edb73d8192bfa034eea3e5..e0c2229cc93226f611ad7c7d7a2c88c0c6bea7d5 100644
Binary files a/doc/fdelmodcManual.pdf and b/doc/fdelmodcManual.pdf differ
diff --git a/fdelmodc/demo/fdelmodc_plane_txt.scr b/fdelmodc/demo/fdelmodc_plane_txt.scr
new file mode 100755
index 0000000000000000000000000000000000000000..d2d687b835dd330e843903dd53c3a163ea2ffd68
--- /dev/null
+++ b/fdelmodc/demo/fdelmodc_plane_txt.scr
@@ -0,0 +1,98 @@
+#!/bin/bash
+#PBS -N fdelmodc
+#PBS -k eo
+#PBS -j eo
+#
+# Models plane wave at depth to receivers at the surface, including snapshots
+export PATH=../../utils:$PATH:
+
+makewave file_out=wavelet.su dt=0.001 nt=1024 fp=13 shift=1 w=g2 verbose=1
+
+makemod file_base=model.su \
+ cp0=1500 ro0=1000 sizex=2100 sizez=2100 \
+ dx=3 dz=3 orig=0,0 \
+ intt=def poly=0 cp=1650 ro=2000 \
+ x=0,2100 z=500,500 gradcp=0.5 grad=100 \
+ intt=def poly=1 cp=1800 ro=2500 \
+ x=0,800,1200,2100 z=900,1400,1400,1200 gradcp=0 grad=0 \
+ verbose=4
+
+export filecp=model_cp.su
+export filecs=model_cs.su
+export filero=model_ro.su
+
+export OMP_NUM_THREADS=1
+
+rm Src.txt
+for i in `seq 1 700 `; do (( x = i*3 )); (( z = i*3 )); echo $x $z >> Src.txt ; done
+
+time ../fdelmodc \
+ file_cp=$filecp file_den=$filero \
+ ischeme=1 \
+ file_src=wavelet.su verbose=4 \
+ file_rcv=rec.su \
+ file_snap=snap.su \
+ xrcv1=0 xrcv2=2100 dxrcv=15 \
+ zrcv1=400 zrcv2=400 \
+ rec_type_vx=1 rec_type_pp=1 rec_type_ss=1 rec_int_vx=1 \
+ dtrcv=0.004 \
+ src_txt=Src.txt \
+ src_type=1 tmod=3.0 src_velo=1800 src_angle=5 \
+ ntaper=21 src_window=11 \
+ left=4 right=4 bottom=4 top=4 \
+ tsnap1=0.1 tsnap2=3.0 dtsnap=0.1 \
+ sna_type_ss=1 sna_type_pp=1
+
+exit;
+
+# to show a movie of the snapshots
+#suxmovie < snap_svz.su perc=99 loop=1
+
+# to reproduce the images in the manual use:
+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]" > SrcRecPositions.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
+
+supsimage < rec_rvz.su \
+ wbox=3 hbox=4 titlesize=-1 labelsize=10 clip=1e-10 verbose=1 \
+ label1="time [s]" label2="lateral position [m]" > rec_plane_rvz.eps
+
+supsimage < rec_rpp.su \
+ wbox=3 hbox=4 titlesize=-1 labelsize=10 clip=1e-11 verbose=1 \
+ label1="time [s]" label2="lateral position [m]" > rec_plane_rpp.eps
+
+supsimage < rec_rss.su \
+ wbox=3 hbox=4 titlesize=-1 labelsize=10 clip=1e-11 verbose=1 \
+ label1="time [s]" label2="lateral position [m]" > rec_plane_rss.eps
+
+#snapshots
+for file in snap_svz snap_spp snap_sss; do
+
+suwind < $file.su key=fldr min=$ifldr max=$ifldr > nep1.su
+while (( ifldr < 12 ))
+do
+(( ifldr += 4 ))
+echo $ifldr
+suwind < $file.su key=fldr min=$ifldr max=$ifldr > nep2.su
+susum nep2.su nep1.su > snaps.su
+mv snaps.su nep1.su
+done
+
+supsimage < nep1.su \
+ wbox=4 hbox=4 titlesize=-1 labelsize=10 verbose=1 \
+ d2=3 f2=0 perc=99 \
+ label1="depth [m]" label2="lateral position [m]" > ${file}_snap.eps
+done
+
diff --git a/fdelmodc/demo/matlab/FD_matlab_interface.m b/fdelmodc/demo/matlab/FD_matlab_interface.m
new file mode 100644
index 0000000000000000000000000000000000000000..da2ff4e44e454de8bacadfaa831d2ef226d61507
--- /dev/null
+++ b/fdelmodc/demo/matlab/FD_matlab_interface.m
@@ -0,0 +1,97 @@
+function [ P, Vz] = FD_mod( xS, xR, tmod, dtrcv, dx, cgrid, rhogrid, orig)
+%Summary of this function goes here
+% Detailed explanation goes here
+
+% save Velocity and density grid
+dims = size(cgrid);
+fileID =fopen('mod_cp.bin','w+','l');
+fwrite(fileID,cgrid,'float32');
+fclose(fileID);
+fileID =fopen('mod_ro.bin','w+','l');
+fwrite(fileID,rhogrid,'float32');
+fclose(fileID);
+
+% Compute sizes for makemod
+sizez=(dims(1)-1)*dx;
+sizex=(dims(2)-1)*dx;
+origz=orig(1);
+origx=orig(2);
+zsrc=xS(2);
+xsrc=xS(1);
+
+%write receiver arrays to file
+dlmwrite('rcv.txt',xR, ' ');
+
+%compute dt for modeling dt < 0.606*dx/Cmax
+Cmax=max(cgrid(:));
+dxmod=dx;
+dtmax=0.606*dxmod/Cmax;
+Sdt=ceil(dtrcv/dtmax);
+dtmod=dtrcv/(Sdt);
+ntwave=16384;
+fmax=0.8/(2*dtrcv); % fmax is 80% of Nyquist frequency
+frig=0.75/(2*dtrcv);
+flef=0.05/(2*dtrcv);
+%tmod=(ntrcv-Sdt+1)*dtrcv;
+
+fileID = fopen('run.scr','w+');
+fprintf(fileID,'#!/bin/bash\n');
+fprintf(fileID,'export PATH=$HOME/src/OpenSource/bin:/opt/CWP/bin/:.:$PATH\n');
+fprintf(fileID,'which fdelmodc\n');
+%fprintf(fileID,'set -x\n');
+fprintf(fileID,'dt=%e\n',dtmod);
+fprintf(fileID,'suaddhead < mod_ro.bin ntrpr=%d ns=%d | \\\n',dims(2), dims(1));
+fprintf(fileID,'sushw key=f1,f2,d1,d2,gx,scalco a=%f,%f,%f,%f,%d,-1000 b=0,0,0,0,%d,0 > mod_ro.su\n',origz, origx, dx, dx, int32(origx*1000), int32(dx*1000));
+fprintf(fileID,'suaddhead < mod_cp.bin ntrpr=%d ns=%d | \\\n',dims(2), dims(1));
+fprintf(fileID,'sushw key=f1,f2,d1,d2,gx,scalco a=%f,%f,%f,%f,%d,-1000 b=0,0,0,0,%d,0 > mod_cp.su\n',origz, origx, dx, dx, int32(origx*1000), int32(dx*1000));
+fprintf(fileID,'makewave w=fw fmin=0 flef=%f frig=%f fmax=%f dt=$dt file_out=wavefw.su nt=%d shift=1 scale=0 scfft=1 verbose=1 >& nep\n', flef, frig, fmax, ntwave);
+fprintf(fileID,'t0=`grep shift nep | awk ''{print $6}''`\n');
+fprintf(fileID,'echo rec_delay for shift in wavelet: t0=$t0\n');
+fprintf(fileID,'tmod=$(echo "scale=4; %f+${t0}" | bc -l)\n',tmod);
+fprintf(fileID,'export filecp=mod_cp.su\n');
+fprintf(fileID,'export filero=mod_ro.su\n');
+fprintf(fileID,'export OMP_NUM_THREADS=4\n');
+fprintf(fileID,'fdelmodc \\\n');
+fprintf(fileID,'file_cp=$filecp file_den=$filero \\\n');
+fprintf(fileID,'ischeme=1 \\\n');
+fprintf(fileID,'file_src=wavefw.su verbose=1 \\\n');
+fprintf(fileID,'dt=$dt \\\n');
+fprintf(fileID,'file_rcv=recv.su \\\n');
+fprintf(fileID,'rec_type_vz=1 rec_type_p=1 rec_int_vz=2 \\\n');
+fprintf(fileID,'rcv_txt=rcv.txt \\\n');
+fprintf(fileID,'dtrcv=%e \\\n', dtrcv);
+fprintf(fileID,'xsrc=%f zsrc=%f \\\n', xsrc, zsrc);
+fprintf(fileID,'src_type=1 tmod=$tmod rec_delay=$t0 \\\n');
+fprintf(fileID,'ntaper=100 \\\n');
+fprintf(fileID,'left=2 right=2 bottom=2 top=2 \n\n');
+fprintf(fileID,'\n');
+fprintf(fileID,'sustrip < recv_rp.su > recv_rp.bin\n');
+fprintf(fileID,'sustrip < recv_rvz.su > recv_rvz.bin\n');
+fprintf(fileID,'surange < recv_rp.su | grep ns | awk ''{print $2}'' > samples\n');
+fprintf(fileID,'surange < recv_rp.su | grep traces | awk ''{print $1}'' > traces\n');
+fclose(fileID);
+!chmod +x run.scr
+system('./run.scr');
+
+path = getenv('PATH');
+path = [path ':$HOME/src/OpenSource/bin:/opt/CWP/bin/:.:'];
+setenv('PATH', path);
+
+% get number of samples and traces
+ns=dlmread('samples');
+ntr=dlmread('traces');
+
+% Pressure field
+file='recv_rp.bin';
+fid=fopen(file,'r');
+P=fread(fid,[ns,ntr],'float32');
+fclose(fid);
+
+% Particle velocity field
+file='recv_rvz.bin';
+fid=fopen(file,'r');
+Vz=fread(fid,[ns,ntr],'float32');
+fclose(fid);
+
+end
+
diff --git a/fdelmodc/demo/matlab/test_matlab_interface.m b/fdelmodc/demo/matlab/test_matlab_interface.m
new file mode 100644
index 0000000000000000000000000000000000000000..7046233db61ebecbfc62ce21a53c48a73fd9ca79
--- /dev/null
+++ b/fdelmodc/demo/matlab/test_matlab_interface.m
@@ -0,0 +1,47 @@
+clear all; clc; close all; clear workspace
+
+% set number of dimensions
+global nDIM; nDIM = 2; % set dimension of space
+
+ % set up spatial grid
+ N1fd = 1280; N2fd = 1280; dxfd =2.13313822/5;
+ x1fd = -(N1fd+1)*dxfd/2 + (1:N1fd)*dxfd;
+ x2fd = -(N2fd+1)*dxfd/2 + (1:N2fd)*dxfd;
+ orig = [-(N1fd+1)*dxfd/2,-(N2fd+1)*dxfd/2];
+ [X1fd,X2fd] = ndgrid(x1fd,x2fd);
+
+ % load model
+% filn=sprintf('sos_fidmod.bin'); fid = fopen(filn,'r'); sos_fd = fread(fid,[N1fd*N2fd],'float32'); fclose('all');
+% filn=sprintf('rho_fidmod.bin'); fid = fopen(filn,'r'); rho_fd = fread(fid,[N1fd*N2fd],'float32'); fclose('all');
+% sos_fd = reshape(sos_fd,[N1fd,N2fd]);
+% rho_fd = reshape(rho_fd,[N1fd,N2fd]);
+ sos_fd = ones(N1fd,N2fd)*1500; % speed of sound
+ rho_fd = ones(N1fd,N2fd)*1500; % density
+
+
+ % time parameters
+ Ntdf = 1024; dtdf = 10^(-3);
+
+ % set up acquisition grid
+ r_rec = 200; % radius of circle
+ Nr = 250; % number of receivers
+ rcvr_phi(1:Nr) = (1:Nr) * 2*pi/Nr; % angles
+ xR = zeros(2,Nr);
+ xR(1,1:Nr) = r_rec * cos(rcvr_phi);
+ xR(2,1:Nr) = r_rec * sin(rcvr_phi);
+ xS = xR(:,1) % choose source at first position
+
+ % plot the impedance and acquisition geometry
+ figure; imagesc(x1fd,x2fd,sos_fd.*rho_fd);
+ hold on; scatter(xR(1,:),xR(2,:),'*b');
+ hold on; scatter(xS(1,:),xS(2,:),'*r');
+ xlabel('x (m)');
+ ylabel('y (m)');
+
+ % Call finite difference code
+ [P, Vz]=FD_mod( xS.', xR.', 0.5, dtdf, dxfd, sos_fd, rho_fd, orig);
+
+ % make a plot
+ figure;imagesc(P);
+ xlabel('angle (degrees)');
+ ylabel('time (t)');
diff --git a/marchenko/demo/ScientificReports/README b/marchenko/demo/ScientificReports/README
new file mode 100644
index 0000000000000000000000000000000000000000..cbecad73b225273369f72c33d01774720013e790
--- /dev/null
+++ b/marchenko/demo/ScientificReports/README
@@ -0,0 +1,61 @@
+The scripts in this directory create Figure 1 of the paper in Scientific Reports
+
+Forward model the data
+1a) model.scr computes the gridded model
+ => tutodel_dx0.5_ro.su, tutodel_dx0.5_cp.su
+1b) shots_slurm.scr creates jobs to model the shots and submit jobs to slurm workload manager
+ => shotsnew/shots_x*_rp.su : ranging from -3000 to 3000 with dxshot=10
+ => to model one shots takes ~25 minutes, 601 shots are modeled
+1c) check.scr after the jobs on shots_*.scr are finished checks if all shots are there
+2) direct.scr creates the direct arrival to be removed from the shots
+ => direct_rp.su
+3) remove_direct.scr remove the direct wave from the shots
+ => shotsnew/refl_rp.su
+4) initialFocus.scr model G_d the intitial focusing function
+ => iniFocusz800x0_rp.su
+
+Apply the marchenko method to compute f2
+5) marchenko.scr perform the Marchenko scheme
+ => f2.su, f1plus.su, f1min.su, Gmin.su, Gplus.su, pgreen.su
+
+Backpropagation of the recorded data from one side and all-sides: snaphots of the wavefield are recorded
+6) back_injrate_planes.scr
+ => snapinj_planevzvxsum_sp.su : backpropagated from all 4 sides of the model
+ => snapinj_surf_sp.su : backpropagated from surface only
+
+Backpropagation of the Marchenko computed result
+7) backpropf2.scr
+ => snapinj_f2_sp.su : backpropagated f2
+
+Generate the postscript files of Figure 1
+8) epsBack.scr
+ => results of backpropagation from all 4 sides: first column of Figure 1
+snapinj_planevzvxsum_-0.60.eps
+snapinj_planevzvxsum_-0.30.eps
+snapinj_planevzvxsum_-0.03.eps
+snapinj_planevzvxsum_0.00.eps
+snapinj_planevzvxsum_0.03.eps
+snapinj_planevzvxsum_0.30.eps
+snapinj_planevzvxsum_0.60.eps
+ => results of backpropagation from surface only: second column of Figure 1
+snapinj_surf_-0.60.eps
+snapinj_surf_-0.30.eps
+snapinj_surf_-0.03.eps
+snapinj_surf_0.00.eps
+snapinj_surf_0.03.eps
+snapinj_surf_0.30.eps
+snapinj_surf_0.60.eps
+ => results of backpropagation of f2 from surface only: third column of Figure 1
+snapinj_f2_-0.60.eps
+snapinj_f2_-0.30.eps
+snapinj_f2_-0.03.eps
+snapinj_f2_0.00.eps
+snapinj_f2_0.03.eps
+snapinj_f2_0.30.eps
+snapinj_f2_0.60.eps
+ => results of symmetrized backpropagation of f2 from surface only: fourth column of Figure 1
+snapinj_f2sum_0.00.eps
+snapinj_f2sum_0.03.eps
+snapinj_f2sum_0.30.eps
+snapinj_f2sum_0.60.eps
+
diff --git a/marchenko/demo/ScientificReports/back_injrate_planes.scr b/marchenko/demo/ScientificReports/back_injrate_planes.scr
new file mode 100755
index 0000000000000000000000000000000000000000..e5e0ebdb86d52a9334a3f9f36e886aa3caa9e106
--- /dev/null
+++ b/marchenko/demo/ScientificReports/back_injrate_planes.scr
@@ -0,0 +1,141 @@
+#!/bin/bash
+
+export PATH=:$HOME/src/OpenSource/bin:$HOME/bin:$PATH:
+
+export OMP_NUM_THREADS=8
+dt=0.00050
+dx=2.5
+
+dt=0.00010
+dx=0.5
+
+makewave fp=30 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
+./model.scr $dx
+ntr=$(echo "scale=0; 1+6000/${dx}" | bc -l)
+echo $ntr
+
+ix1a=1
+ix1b=$(echo "scale=0; ${ix1a}+6000/${dx}" | bc -l)
+ix2a=$(echo "scale=0; ${ix1b}+1" | bc -l)
+ix2b=$(echo "scale=0; ${ix2a}+6000/${dx}" | bc -l)
+ix3a=$(echo "scale=0; ${ix2b}+1" | bc -l)
+ix3b=$(echo "scale=0; ${ix3a}+1200/${dx}" | bc -l)
+ix4a=$(echo "scale=0; ${ix3b}+1" | bc -l)
+ix4b=$(echo "scale=0; ${ix4a}+1200/${dx}" | bc -l)
+
+file_cp=tutodel_dx${dx}_cp.su
+file_ro=tutodel_dx${dx}_ro.su
+
+#model data to be propagated back into medium
+fdelmodc \
+ file_cp=$file_cp ischeme=1 iorder=4 \
+ file_den=$file_ro \
+ file_src=wave.su \
+ file_rcv=inj_rate_plane_dx${dx}.su \
+ src_type=1 \
+ src_injectionrate=1 \
+ rec_type_vz=1 \
+ rec_type_vx=1 \
+ rec_type_p=0 \
+ rec_int_vz=2 \
+ rec_int_vx=2 \
+ dtrcv=$dt \
+ rec_delay=0.1 \
+ verbose=1 \
+ tmod=4.4000 \
+ xrcv1=-3000,-3000,-3000,3000 xrcv2=3000,3000,-3000,3000 zrcv1=0,1200,0,0 zrcv2=0,1200,1200,1200 \
+ dxrcv=$dx,$dx,0,0 dzrcv=0,0,$dx,$dx \
+ xsrc=0 zsrc=800 \
+ ntaper=400 \
+ left=2 right=2 top=2 bottom=2
+
+
+##################################
+# back propagate from all sides
+# scale with -1 for outward pointing vector
+scale=-1.0
+suwind key=tracl min=$ix1a max=$ix1b < inj_rate_plane_dx${dx}_rvz.su | sugain scale=$scale > inj_rate_plane_dx${dx}vz.su
+suwind key=tracl min=$ix2a max=$ix2b < inj_rate_plane_dx${dx}_rvz.su >> inj_rate_plane_dx${dx}vz.su
+
+suwind < inj_rate_plane_dx${dx}_rvx.su key=tracl min=$ix3a max=$ix3b | sugain scale=$scale > inj_rate_plane_dx${dx}vx.su
+suwind < inj_rate_plane_dx${dx}_rvx.su key=tracl min=$ix4a max=$ix4b >> inj_rate_plane_dx${dx}vx.su
+
+# at 4.3000 seconds (tmod - rec_delay) the focus is at t=0
+fdelmodc \
+ file_cp=$file_cp ischeme=1 iorder=4 \
+ file_den=$file_ro \
+ file_src=inj_rate_plane_dx${dx}vz.su \
+ file_snap=snapinj_planevz.su \
+ grid_dir=1 \
+ src_type=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_p=0 \
+ dtrcv=$dt \
+ rec_delay=0.0 \
+ verbose=3 \
+ tmod=5.004 \
+ tsnap1=3.000 tsnap2=5.00 dtsnap=0.01 \
+ xsnap1=-1000 xsnap2=1000 dxsnap=2.5 dzsnap=2.5 \
+ sna_type_vz=0 \
+ sna_type_p=1 \
+ ntaper=400 \
+ left=2 right=2 top=2 bottom=2
+
+# tsnap1=4.200 tsnap2=4.50 dtsnap=0.004 \
+
+fdelmodc \
+ file_cp=$file_cp ischeme=1 iorder=4 \
+ file_den=$file_ro \
+ file_src=inj_rate_plane_dx${dx}vx.su \
+ file_snap=snapinj_planevx.su \
+ grid_dir=1 \
+ src_type=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_p=0 \
+ dtrcv=$dt \
+ rec_delay=0.0 \
+ verbose=1 \
+ tmod=5.004 \
+ tsnap1=3.000 tsnap2=5.00 dtsnap=0.01 \
+ xsnap1=-1000 xsnap2=1000 dxsnap=2.5 dzsnap=2.5 \
+ sna_type_vz=0 \
+ sna_type_p=1 \
+ ntaper=400 \
+ left=2 right=2 top=2 bottom=2
+
+# tsnap1=4.200 tsnap2=4.50 dtsnap=0.004 \
+
+suop2 snapinj_planevz_sp.su snapinj_planevx_sp.su op=sum w1=1 w2=1 > snapinj_planevzvxsum_sp.su
+
+##################################
+# back propagate from surface only
+scale=-1.0
+suwind key=tracl min=1 max=$ntr < inj_rate_plane_dx${dx}_rvz.su | sutaper tr1=100 tr2=100 ntr=$ntr | sugain scale=$scale > inj_rate_surf_dx${dx}_rvz.su
+
+# at 4.3000 seconds the focus is at t=0
+fdelmodc \
+ file_cp=$file_cp ischeme=1 iorder=4 \
+ file_den=$file_ro \
+ file_src=inj_rate_surf_dx${dx}_rvz.su \
+ file_snap=snapinj_surf.su \
+ grid_dir=1 \
+ src_type=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_p=0 \
+ dtrcv=$dt \
+ rec_delay=0.0 \
+ verbose=1 \
+ tmod=5.004 \
+ tsnap1=3.000 tsnap2=5.00 dtsnap=0.01 \
+ xsnap1=-1000 xsnap2=1000 dxsnap=2.5 dzsnap=2.5 \
+ sna_type_vz=0 \
+ sna_type_p=1 \
+ ntaper=400 \
+ left=2 right=2 top=2 bottom=2
+
diff --git a/marchenko/demo/ScientificReports/backpropf2.scr b/marchenko/demo/ScientificReports/backpropf2.scr
new file mode 100755
index 0000000000000000000000000000000000000000..a19fca25cc61b0157ecdaf86b0cc9c52d3ce726e
--- /dev/null
+++ b/marchenko/demo/ScientificReports/backpropf2.scr
@@ -0,0 +1,57 @@
+#!/bin/bash
+
+export PATH=:$HOME/src/OpenSource/bin:$HOME/bin:$PATH:
+
+export OMP_NUM_THREADS=8
+dx=2.5
+dt=0.00050
+dx=0.5
+dt=0.000125
+
+ns=`surange < f2.su | grep ns | awk '{print $2}'`
+dtrcv=`surange < f2.su | grep dt | awk '{print $2/1000000.0}'`
+#suwind key=gx min=-2250000 max=2250000 itmax=1023 < f2.su | sushw key=f2 a=-2250 > nep.su
+#shift=$(echo "scale=6; ($dtrcv*($ns/2.0-1)+0.5*$dt-0.000250)" | bc -l)
+shift=$(echo "scale=6; ($dtrcv*($ns/2.0-1))" | bc -l)
+echo $shift
+basop choice=shift shift=$shift file_in=f2.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
+# this is now done in fdelmodc
+
+#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
+
+file_cp=tutodel_dx${dx}_cp.su
+file_ro=tutodel_dx${dx}_ro.su
+
+# back propagate with f2 from marchenko.scr
+
+set -x
+fdelmodc \
+ file_cp=$file_cp ischeme=1 iorder=4 \
+ file_den=$file_ro \
+ file_src=pplus.su \
+ file_snap=snapinj_f2.su \
+ grid_dir=0 \
+ src_type=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=0 \
+ rec_int_p=0 \
+ verbose=4 \
+ dt=$dt \
+ tmod=3.0000 \
+ tsnap1=0.744125 tsnap2=2.744125 dtsnap=0.01 \
+ xsnap1=-1000 xsnap2=1000 dxsnap=2.5 dzsnap=2.5 \
+ sna_type_vz=0 \
+ sna_type_p=1 \
+ fmax=120 \
+ ntaper=400 \
+ left=2 right=2 top=2 bottom=2
+
+#2.044125
+# tmod=2.1000 \
+# tsnap1=2.0400 tsnap2=2.0500 dtsnap=0.000125 \
diff --git a/marchenko/demo/ScientificReports/check.scr b/marchenko/demo/ScientificReports/check.scr
new file mode 100755
index 0000000000000000000000000000000000000000..b369721f40779a714ff6a020a404bacfc846e21c
--- /dev/null
+++ b/marchenko/demo/ScientificReports/check.scr
@@ -0,0 +1,29 @@
+#!/bin/bash
+#PBS -N fdelmod
+#PBS -q long
+#PBS -l nodes=1
+#PBS -k eo
+#PBS -j eo
+
+dxshot=10
+ishot=0
+nshots=601
+zsrc=0
+
+while (( ishot < nshots ))
+do
+
+ (( xsrc = -3000 + ${ishot}*${dxshot} ))
+
+ file_rcv=shotsnew/shots_${xsrc}_rp.su
+
+ if [ ! -e "$file_rcv" ]
+ then
+ echo $xsrc is missing
+ sbatch jobs/slurm_$ishot.job
+ fi
+
+ (( ishot = $ishot + 1))
+
+done
+
diff --git a/marchenko/demo/ScientificReports/clean b/marchenko/demo/ScientificReports/clean
new file mode 100755
index 0000000000000000000000000000000000000000..26987b43f2604488b03655488d6daaad171eb7ad
--- /dev/null
+++ b/marchenko/demo/ScientificReports/clean
@@ -0,0 +1,4 @@
+#!/bin/bash
+
+rm *.su *.bin *.eps nep line* *.asci curve* *.mod
+
diff --git a/marchenko/demo/ScientificReports/direct.scr b/marchenko/demo/ScientificReports/direct.scr
new file mode 100755
index 0000000000000000000000000000000000000000..fa06b4215d8a8908a7bcbe37ff16177f265d146a
--- /dev/null
+++ b/marchenko/demo/ScientificReports/direct.scr
@@ -0,0 +1,40 @@
+#!/bin/bash
+
+cd /vardim/home/thorbcke/data/Kees/Marchenko/Tutorial
+
+dx=2.5
+dt=0.00050
+fast="fast"
+
+dx=0.5
+dt=0.0001
+fast=""
+
+makemod sizex=12000 sizez=250 dx=$dx dz=$dx cp0=1500 ro0=1000 \
+ orig=-6000,0 file_base=noContrast.su
+
+export OMP_NUM_THREADS=8
+makewave w=fw fmin=0 flef=5 frig=80 fmax=100 dt=$dt file_out=wavefw${fast}.su nt=8192 t0=0.4 scale=0 scfft=1
+
+fdelmodc \
+ file_cp=noContrast_cp.su ischeme=1 iorder=4 \
+ file_den=noContrast_ro.su \
+ file_src=wavefw${fast}.su \
+ file_rcv=direct${fast}.su \
+ src_type=7 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ rec_delay=0.4 \
+ dtrcv=0.004 \
+ verbose=2 \
+ tmod=4.492 \
+ dxrcv=10.0 \
+ xrcv1=-6000 xrcv2=6000 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=0 zsrc=0 \
+ ntaper=200 \
+ left=2 right=2 top=2 bottom=2
+
diff --git a/marchenko/demo/ScientificReports/epsBack.scr b/marchenko/demo/ScientificReports/epsBack.scr
new file mode 100755
index 0000000000000000000000000000000000000000..fad63523dfa74f535cbaf7dd71847af4d4b3349c
--- /dev/null
+++ b/marchenko/demo/ScientificReports/epsBack.scr
@@ -0,0 +1,131 @@
+#!/bin/bash
+
+export PATH=:$HOME/src/OpenSource/bin:$HOME/bin:$PATH:
+
+dx=0.5
+file_snap_all=snapinj_planevzvxsum_sp.su
+file_snap_surf=snapinj_surf_sp.su
+file_snap_f2=snapinj_f2_sp.su
+
+
+sumax < ${file_snap_all} mode=abs outpar=nep
+clip_all=`cat nep | awk '{print $1/25}'`
+clip_t0all=`cat nep | awk '{print $1/10}'`
+sumax < ${file_snap_surf} mode=abs outpar=nep
+clip_surf=`cat nep | awk '{print $1/10}'`
+clip_t0surf=`cat nep | awk '{print $1/4}'`
+sumax < ${file_snap_f2} mode=abs outpar=nep
+clip_f2=`cat nep | awk '{print $1/10}'`
+clip_t0f2=`cat nep | awk '{print $1/4}'`
+echo $clip_all $clip_surf $clip_f2
+
+#131 => t=0
+for file_snap in $file_snap_all $file_snap_surf $file_snap_f2
+do
+ if [ $file_snap == $file_snap_all ]; then clip=$clip_all; fi;
+ if [ $file_snap == $file_snap_surf ]; then clip=$clip_surf; fi;
+ if [ $file_snap == $file_snap_f2 ]; then clip=$clip_f2; fi;
+ echo $clip
+ file_base=${file_snap%_sp.su}
+ suwind key=fldr min=61 max=201 < $file_snap | sugain scale=-1 | sustrip > $file_base.bin
+ for fldr in 71 101 128 131 134 161 191;
+ do
+ if [ $fldr == 131 ];
+ then
+ clipt=$(echo "scale=2; ${clip}*4.0" | bc -l);
+ else
+ clipt=$clip;
+ fi;
+ echo clip=$clip clipt=$clipt
+ times=$(echo "scale=2; 0.01*(${fldr}-131)" | bc -l)
+ atime=`printf "%4.2f" $times`
+ suwind key=fldr min=$fldr max=$fldr < ${file_snap} | \
+ supsimage hbox=2.5 wbox=4 labelsize=-1 \
+ x1beg=0 x1end=1250.0 clip=$clipt \
+ curve=curve1,curve2,curve3 npair=25,25,25 curvecolor=black,black,black curvedash=3,3,3 \
+ n1tic=5 x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > ${file_base}_$atime.eps
+ done
+done
+
+file_base=${file_snap_f2%_sp.su}
+rm ${file_base}sum.su
+for shot in `seq 61 1 131`;
+do
+ suwind key=fldr min=$shot max=$shot < ${file_snap_f2} > neg.su
+ (( a = 131+(131-$shot) ))
+ echo $shot $a
+ suwind key=fldr min=$a max=$a < ${file_snap_f2} > pos.su
+ susum neg.su pos.su | sushw key=fldr a=$shot >> ${file_base}sum.su
+done
+
+for shot in `seq 132 1 201`;
+do
+ suwind key=fldr min=$shot max=$shot < ${file_snap_f2}> pos.su
+ (( a = 131+(131-$shot) ))
+ echo $shot $a
+ suwind key=fldr min=$a max=$a < ${file_snap_f2} > neg.su
+ susum neg.su pos.su | sushw key=fldr a=$fldr >> ${file_base}sum.su
+done
+
+sugain scale=-1 < ${file_base}sum.su | sustrip > ${file_base}sum.bin
+
+exit
+
+#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_base=${file_snap%_sp.su}
+for fldr in 71 101 128 131;
+do
+ times=$(echo "scale=2; 0.01*(131-${fldr})" | bc -l)
+ atime=`printf "%4.2f" $times`
+ suwind key=fldr min=$fldr max=$fldr < ${file_snap_f2} > neg.su
+ (( fldr = 131+(131-$fldr) ))
+ suwind key=fldr min=$fldr max=$fldr < ${file_snap_f2}> pos.su
+ susum neg.su pos.su | \
+ supsimage hbox=2.5 wbox=4 labelsize=-1 \
+ x1beg=0 x1end=1250.0 clip=$clip_f2 \
+ curve=curve1,curve2,curve3 npair=25,25,25 curvecolor=black,black,black curvedash=3,3,3 \
+ n1tic=5 x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > ${file_base}sum_$atime.eps
+done
+
+#supsgraph < wave.su hbox=6 wbox=2 labelsize=10 f1=0 x1end=0.4 > wave.eps
+#basop file_in=wave.su choice=15 | supsgraph hbox=6 wbox=2 labelsize=10 f1=0 x1end=0.4 > wavesjom.eps
+
+file=inj_rate_surf_dx${dx}_rvz.su
+file_base=${file%.su}
+sumax < ${file} mode=abs outpar=nep
+clip=`cat nep | awk '{print $1/5}'`
+
+ns=`surange < $file | grep ns | awk '{print $2}'`
+dtrcv=`surange < $file | grep dt | awk '{print $2/1000000.0}'`
+shift=$(echo "scale=6; ($dtrcv*($ns/2.0-1))" | bc -l)
+skip=$(echo "scale=0; (10/${dx})" | bc -l)
+echo $shift $skip
+basop choice=shift shift=$shift file_in=$file verbose=1 | suwind s=1 j=$skip | sushw key=d2,tracl a=10,1 b=0,1 > nep.su
+supsimage hbox=2.5 wbox=4 labelsize=-1 < nep.su \
+ x1beg=-0.6 x1end=1.5 clip=$clip grid1=dot \
+ n1tic=5 x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > ${file_base}.eps
+
+suwind s=1 j=10 < nep.su | \
+ supswigp hbox=2.5 wbox=4 labelsize=-1 d2=100 \
+ x1beg=-1.5 x1end=1.5 xcur=1.5 grid1=dot n1tic=1 d1num=2.145 \
+ f2=-3000 x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > ${file_base}_wiggle.eps
+
+file=pplus.su
+file_base=${file%.su}
+sumax < ${file} mode=abs outpar=nep
+clip=`cat nep | awk '{print $1/5}'`
+
+suflip flip=3 < $file | \
+ supsimage hbox=2.5 wbox=4 labelsize=-1 x2beg=-1000 x2end=1000 \
+ x1beg=-0.6 x1end=1.5 clip=$clip grid1=dot \
+ n1tic=5 x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > ${file_base}.eps
+
+suflip flip=3 < $file | suwind s=1 j=10 | \
+ supswigp hbox=2.5 wbox=4 labelsize=-1 d2=100 \
+ x1beg=-1.5 x1end=1.5 xcur=1.5 grid1=dot n1tic=1 d1num=2.044 \
+ f2=-3000 x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > ${file_base}_wiggle.eps
+
+
+
diff --git a/marchenko/demo/ScientificReports/initialFocus.scr b/marchenko/demo/ScientificReports/initialFocus.scr
new file mode 100755
index 0000000000000000000000000000000000000000..4353514c42aafbaaa8a4799b4a31c68a43d3be36
--- /dev/null
+++ b/marchenko/demo/ScientificReports/initialFocus.scr
@@ -0,0 +1,86 @@
+#!/bin/bash
+#PBS -N fdelmod
+#PBS -q verylong
+#PBS -l nodes=1
+#PBS -k eo
+#PBS -j eo
+
+export PATH=$HOME/bin64:$HOME/src/OpenSource/utils:$PATH:
+which makewave
+which makemod
+
+dx=2.5
+dt=0.0005
+dx=0.5
+dt=0.0001
+
+makemod sizex=6000 sizez=810 dx=$dx dz=$dx cp0=1500 ro0=1000 \
+ orig=-3000,0 file_base=tutodeldown.su verbose=2 \
+ intt=def x=-3000,-2000,-1000,0,1000,2000,3000 z=240,130,250,300,350,380,320 poly=2 cp=1950 ro=4500 grad=0 \
+ intt=def x=-3000,-2200,-1500,0,1300,2100,3000 z=620,640,590,600,740,700,600 poly=2 cp=2000 ro=1400 grad=0
+
+file_cp=tutodeldown_cp.su
+file_ro=tutodeldown_cp.su
+
+makewave fp=30 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
+
+dxshot=10
+ishot=300
+nshots=301
+
+export OMP_NUM_THREADS=1
+mkdir -p shots${fast}
+mkdir -p jobs
+
+while (( ishot < nshots ))
+do
+
+ (( xsrc = -3000 + ${ishot}*${dxshot} ))
+ echo xsrc=$xsrc
+ file_rcv=iniFocusz800x${xsrc}.su
+
+ cat << EOF > jobs/job_$ishot.slurm
+#!/bin/bash
+#
+#SBATCH -J marchenko46
+#SBATCH --cpus-per-task=8
+#SBATCH --ntasks=1
+#SBATCH --time=0:45:00
+#SBATCH -o iniFocus-%A.out
+
+cd \$SLURM_SUBMIT_DIR
+
+export OMP_NUM_THREADS=8
+
+ fdelmodc \
+ file_cp=$file_cp ischeme=1 iorder=4 \
+ file_den=$file_ro \
+ file_src=wave.su \
+ file_rcv=$file_rcv \
+ src_type=1 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ rec_delay=0.1 \
+ dtrcv=0.004 \
+ verbose=2 \
+ tmod=4.192 \
+ dxrcv=10.0 \
+ xrcv1=-3000 xrcv2=3000 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=$xsrc zsrc=800 \
+ ntaper=200 \
+ left=2 right=2 top=2 bottom=2
+EOF
+
+ chmod +x jobs/job_$ishot.slurm
+ jobs/job_$ishot.slurm
+ #sbatch jobs/job_$ishot.slurm
+
+ (( ishot = $ishot + 1))
+done
+
+
+
diff --git a/marchenko/demo/ScientificReports/marchenko.scr b/marchenko/demo/ScientificReports/marchenko.scr
new file mode 100755
index 0000000000000000000000000000000000000000..427f2c096b3da9633d1f0a4bee3c86007f2d17a4
--- /dev/null
+++ b/marchenko/demo/ScientificReports/marchenko.scr
@@ -0,0 +1,22 @@
+#!/bin/bash -x
+#PBS -N fdelmod
+#PBS -q verylong
+#PBS -l nodes=1
+#PBS -k eo
+#PBS -j eo
+
+export PATH=$HOME/src/OpenSource/bin:$PATH:
+export OMP_NUM_THREADS=8
+
+#for backprop.scr use fast model (with stair-case interfaces)
+#for backprop.scr use pml model (with 0.5 fine-grids)
+
+smooth=3
+
+fmute file_shot=iniFocusz800x0_rp.su file_out=p0plus.su above=-1 shift=-10 verbose=1 check=1 hw=4
+
+marchenko file_shot=shotsnew/refl_rp.su file_tinv=p0plus.su nshots=601 \
+ file_green=pgreen.su verbose=1 tap=0 ntap=10 niter=15 hw=4 shift=10 smooth=$smooth \
+ file_gplus=Gplus.su file_gmin=Gmin.su file_f1plus=f1plus.su file_f1min=f1min.su \
+ file_f2=f2.su fmax=90 file_f1plus=f1plus.su
+
diff --git a/marchenko/demo/ScientificReports/model.scr b/marchenko/demo/ScientificReports/model.scr
new file mode 100755
index 0000000000000000000000000000000000000000..7a37739be2a78183a3f2582b7f7409d7d8001327
--- /dev/null
+++ b/marchenko/demo/ScientificReports/model.scr
@@ -0,0 +1,98 @@
+#!/bin/bash -x
+
+#cd /vardim/home/thorbcke/data/Kees/Marchenko/Nature_SAGA
+
+if [ -n "$1" ]
+then
+ dx=$1
+else
+ dx=2.5
+fi
+
+makemod sizex=6000 sizez=1250 dx=$dx dz=$dx cp0=1500 ro0=1000 \
+ orig=-3000,0 file_base=tutodel_dx$dx.su rayfile=1 skip=100 verbose=2 \
+ intt=def x=-3000,-2000,-1000,0,1000,2000,3000 z=240,130,250,300,350,380,320 poly=2 cp=1950 ro=4500 \
+ intt=def x=-3000,-2200,-1500,0,1300,2100,3000 z=620,640,590,600,740,700,600 poly=2 cp=2000 ro=1400 \
+ intt=def x=-3000,-1800,0,2200,3000 z=920,1000,900,1000,1010 poly=2 cp=2300 ro=1600
+
+suwind key=gx min=0 max=0 < tutodel_dx${dx}_cp.su > tracemodel.su
+
+supsgraph < tracemodel.su hbox=6 wbox=3 labelsize=10 style=seismic \
+ x1beg=0 d1=$dx f1=0 d2num=100 label1="depth (m)" label2="velocity (m/s)"> tracemodel_x_0.eps
+
+supsimage hbox=4 wbox=6 labelsize=10 < tutodel_dx2.5_cp.su \
+ x1beg=0 x1end=1200.0 legend=1 threecolor=1 \
+ d1s=1.0 d2s=1.0 \
+ wrgb=1.0,.5,0 grgb=0,.7,1.0 brgb=0,1.0,0 \
+ bps=24 bclip=2300 wclip=1800 \
+ curve=curve1,curve2,curve3 npair=25,25,25 curvecolor=black,black,black curvedash=3,3,3 \
+ n1tic=5 x2beg=-3000 f2num=-3000 d2num=1000 x2end=3000 > tutodel_cp.eps
+
+#itmax=$(echo "scale=0; (1250/${dx}+1)" | bc -l)
+suwind key=gx min=-1000000 max=1000000 < tutodel_dx${dx}_cp.su | sugain scale=-1 dt=1 | sustrip > tutodelsnap_cp.bin
+
+exit;
+#model reference for SI results
+makewave fp=30 dt=$dt file_out=wave${dx}.su nt=4096 t0=0.1 scale=1
+
+dx=0.5
+file_cp=tutodelbelow_dx${dx}_cp.su
+file_ro=tutodelbelow_dx${dx}_ro.su
+
+fdelmodc \
+ file_cp=$file_cp ischeme=1 iorder=4 \
+ file_den=$file_ro \
+ file_src=wave${dx}.su \
+ file_rcv=reference_${dx}Z800.su \
+ src_type=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_p=0 \
+ dtrcv=0.004 \
+ dxrcv=2.5 \
+ rec_delay=0.1 \
+ verbose=1 \
+ tmod=1.104 \
+ xrcv1=-3000 xrcv2=3000 zrcv1=800 zrcv2=800 \
+ xsrc=0 zsrc=800 \
+ ntaper=400 \
+ left=4 right=4 top=4 bottom=4
+
+sumax < reference_${dx}Z800_rp.su mode=abs outpar=nep
+clip=`cat nep | awk '{print $1/30}'`
+sugain epow=0.8 < reference_${dx}Z800_rp.su | \
+ supsimage hbox=4 wbox=8 labelsize=10 linewidth=0.0 \
+ n1tic=2 d2=2.5 f1=0.0 x1beg=0.0 x1end=1.004 \
+ f2=-3000 f2num=-3000 d2num=1000 clip=$clip > reference_${dx}Z800_rp.eps
+
+file_cp=tutodel_dx${dx}_cp.su
+file_ro=tutodel_dx${dx}_ro.su
+
+fdelmodc \
+ file_cp=$file_cp ischeme=1 iorder=4 \
+ file_den=$file_ro \
+ file_src=wave${dx}.su \
+ file_rcv=green_${dx}Z800.su \
+ src_type=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_p=0 \
+ dtrcv=0.004 \
+ dxrcv=2.5 \
+ rec_delay=0.1 \
+ verbose=1 \
+ tmod=1.104 \
+ xrcv1=-3000 xrcv2=3000 zrcv1=800 zrcv2=800 \
+ xsrc=0 zsrc=800 \
+ ntaper=400 \
+ left=4 right=4 top=4 bottom=4
+
+sumax < green_${dx}Z800_rp.su mode=abs outpar=nep
+clip=`cat nep | awk '{print $1/30}'`
+sugain epow=0.8 < green_${dx}Z800_rp.su | \
+ supsimage hbox=4 wbox=8 labelsize=10 linewidth=0.0 \
+ n1tic=2 d2=2.5 f1=0.0 x1beg=0.0 x1end=1.004 \
+ f2=-3000 f2num=-3000 d2num=1000 clip=$clip > green_${dx}Z800_rp.eps
+
diff --git a/marchenko/demo/ScientificReports/remove_direct.scr b/marchenko/demo/ScientificReports/remove_direct.scr
new file mode 100755
index 0000000000000000000000000000000000000000..71a99c127a193cd2d9c3cc79ef4e72b7ee95f368
--- /dev/null
+++ b/marchenko/demo/ScientificReports/remove_direct.scr
@@ -0,0 +1,42 @@
+#!/bin/bash
+#PBS -N fdelmod
+#PBS -q verylong
+#PBS -l nodes=1
+#PBS -k eo
+#PBS -j eo
+
+export PATH=$HOME/bin:$HOME/src/OpenSource/bin:$PATH:
+
+dxshot=10
+ishot=0
+nshots=601
+
+file_R=shotsnew/refl_rp.su
+rm $file_R
+
+while (( ishot < nshots ))
+do
+
+ (( xsrc = -3000 + ${ishot}*${dxshot} ))
+ (( sx = ${xsrc}*1000 ))
+ (( iishot = ${ishot}*${dxshot}/10 ))
+ (( tr1 = 601 - ${iishot} ))
+ (( tr2 = ${tr1} + 600 ))
+ echo xsrc=$xsrc tr1=$tr1 tr2=$tr2
+
+# direct wave
+ suwind < direct_rp.su key=tracl min=$tr1 max=$tr2 > direct.su
+
+# 2D shot
+ file_rcv=shotsnew/shots_${xsrc}_rp.su
+ #suwind key=tracl min=1 max=601 < $file_rcv > shotz0.su
+ sudiff $file_rcv direct.su > refl.su
+
+ (( ishot = $ishot + 1))
+
+ sushw < refl.su key=fldr a=$ishot | sugain scale=0.5 >> $file_R
+
+done
+
+rm direct.su refl.su
+
diff --git a/marchenko/demo/ScientificReports/shotslurm.scr b/marchenko/demo/ScientificReports/shotslurm.scr
new file mode 100755
index 0000000000000000000000000000000000000000..fb5feee5c5172108e37479485121bc23c6a9c2ba
--- /dev/null
+++ b/marchenko/demo/ScientificReports/shotslurm.scr
@@ -0,0 +1,67 @@
+#!/bin/bash
+
+export PATH=:$HOME/src/OpenSource/bin:$PATH:
+
+dt=0.0001
+makewave w=fw fmin=0 flef=5 frig=80 fmax=100 dt=$dt file_out=wavefw.su nt=8192 t0=0.4 scale=0 scfft=1
+
+./model.scr 0.5
+
+mkdir -p shotsnew
+mkdir -p jobs
+
+dxshot=10
+ishot=0
+nshots=601
+zsrc=0
+
+while (( ishot < nshots ))
+do
+
+ (( xsrc = -3000 + ${ishot}*${dxshot} ))
+
+ echo ishot=$ishot xsrc=$xsrc zsrc=$zsrc
+
+ cat << EOF > jobs/slurm_$ishot.job
+#!/bin/bash
+#SBATCH -J shot_${xsrc}
+#SBATCH --cpus-per-task=6
+#SBATCH --ntasks=1
+#SBATCH --time=03:00:00
+#SBATCH -o shot_${xsrc}.out
+
+cd \$SLURM_SUBMIT_DIR
+
+ export PATH=:\$HOME/src/OpenSource/bin:\$PATH:
+ export OMP_NUM_THREADS=6
+ file_rcv=shotsnew/shots_${xsrc}.su
+
+ fdelmodc \
+ file_cp=tutodel_cp.su ischeme=1 iorder=4 \
+ file_den=tutodel_ro.su \
+ file_src=wavefw.su \
+ file_rcv=\$file_rcv \
+ src_type=7 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ rec_delay=0.4 \
+ dtrcv=0.004 \
+ verbose=2 \
+ tmod=4.492 \
+ dxrcv=10.0 \
+ xrcv1=-3000 xrcv2=3000 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=$xsrc zsrc=$zsrc \
+ ntaper=200 \
+ left=2 right=2 top=2 bottom=2
+EOF
+
+ sbatch jobs/slurm_$ishot.job
+
+ (( ishot = $ishot + 1))
+
+done
+
diff --git a/marchenko/demo/test2/Joeri/mar.scr b/marchenko/demo/test2/Joeri/mar.scr
deleted file mode 100755
index b03a1bcf07436ac0ca51e466b59fc7d876bbf086..0000000000000000000000000000000000000000
--- a/marchenko/demo/test2/Joeri/mar.scr
+++ /dev/null
@@ -1,37 +0,0 @@
-#!/bin/bash -x
-#makewave fp=30 dt=0.004 file_out=wave2.su nt=2048 t0=0.0 scale=1
-#PBS -N fdelmod
-#PBS -q verylong
-#PBS -l nodes=1
-#PBS -k eo
-#PBS -j eo
-
-#export PATH=$HOME/bin64:$HOME/src/OpenSource/utils:$PATH:
-export PATH=$HOME/Marchenko/utils:$PATH:
-which makewave
-which makemod
-which fmute
-which syn2d
-export OMP_NUM_THREADS=8
-
-#cd /vardim/home/bracken/Marchenko/Redatum/Altmod
-
-#raytime file_vel=model10_cp.su xsrc1=0 zsrc1=1550 file_out=rayA.su verbose=1 dxrcv=10
-#q1=0.5
-
-#~/Marchenko/utils/marchenkocost file_shot=../shotsdx_rp.su file_tinv=above.su file_wav=wave2.su file_green=greenA1.su verbose=2 \
-#key=sx tap=1 ntap=0 reci=0 niter=20 hw=5 shift=6 smooth=5 w=1 file_gplus=GplusA1.su file_gmin=GminA1.su file_f1plus=f1plusA1.su \
-#file_f1min=f1minA1.su file_Gm0=Gm0A1.su file_direct=tinvA1.su file_Udk=UdkA1.su startb=1.0 endb=1.0 nb=0
-
-#~/Marchenko/utils/marchenkocost file_shot=../shotsdx_rp.su file_tinv=above.su file_wav=wave2.su file_green=greenA2.su verbose=2 \
-#key=sx tap=1 ntap=0 reci=0 niter=20 hw=5 shift=6 smooth=5 w=1 file_gplus=GplusA2.su file_gmin=GminA2.su file_f1plus=f1plusA2.su \
-#file_f1min=f1minA2.su file_Gm0=Gm0A2.su file_direct=tinvA2.su file_Udk=UdkA2.su startb=2.0 endb=2.0 nb=0
-
-#~/Marchenko/utils/marchenkocost file_shot=../shotsdx_rp.su file_tinv=above.su file_wav=wave2.su file_green=greenA05.su verbose=2 \
-#key=sx tap=1 ntap=0 reci=0 niter=20 hw=5 shift=6 smooth=5 w=1 file_gplus=GplusA05.su file_gmin=GminA05.su file_f1plus=f1plusA05.su \
-#file_f1min=f1minA05.su file_Gm0=Gm0A05.su file_direct=tinvA05.su file_Udk=UdkA05.su startb=0.5 endb=0.5 nb=0
-
-#~bracken/Marchenko/utils/marchenkocost file_shot=shots.su file_tinv=nep.su file_wav=wave2.su file_green=green.su verbose=2 \
-/vardim/home/thorbcke/src/OpenSource/bin/marchenko file_shot=shots.su file_tinv=nep.su file_green=green.su verbose=2 \
-key=sx tap=1 ntap=0 reci=0 niter=30 hw=5 shift=9 smooth=8 w=1 scaling=0 file_gmin=gmin.su \
-file_f1plus=f1plus.su file_gplus=gplus.su
diff --git a/marchenko/demo/test2/README b/marchenko/demo/test2/README
deleted file mode 100644
index e39ec201319654383a97fff0812df5a9dfe7da22..0000000000000000000000000000000000000000
--- a/marchenko/demo/test2/README
+++ /dev/null
@@ -1,198 +0,0 @@
-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 intitial 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.
-+) clean : remove all produced files and start with a clean directory
-
-
-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
-Note if you model the initial Focusing operator also with a w=fw wavelet the length of the wavelet becomes very long. The
-mute-windows applied in Marchenko will then also mute a big part of this very long fw wavelet and will not converge anymore.
-
-
-==> 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 this 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 and generates the file;
- - 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
-
-
-Note that the script epsIterwithLabels.scr produces the same figures, but with axis-labels.
-
---------------------------
-* 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 1 minute. The generated output file is
- - backpropf2_sp.su
-
-The postscript files of Figure 8 are generated with
-
-==> run epsBackprop.scr
-
--- Figure 8 column 1
-backpropf2_-0.30.eps
-backpropf2_-0.15.eps
-backpropf2_-0.03.eps
-backpropf2_-0.02.eps
-backpropf2_0.00.eps
--- Figure 8 column 2
-backpropf2_0.30.eps
-backpropf2_0.15.eps
-backpropf2_0.03.eps
-backpropf2_0.02.eps
-backpropf2_0.00.eps
--- Figure 8 column 3
-backpropf2sum_0.30.eps
-backpropf2sum_0.15.eps
-backpropf2sum_0.03.eps
-backpropf2sum_0.02.eps
-backpropf2_0.00.eps
-
-
-The figures in the appendix, to explain the different options in the programs, are reproduced by
-
-==> run figAppendi.scr
-
--- Figure A-1
-noise_above0.eps
-noise_above1.eps
-noise_above-1.eps
-noise_above2.eps
-noise_above4.eps
-
--- Figure A-2
-iniFocus_shifts.eps
-
diff --git a/marchenko/demo/test2/backpropf2.scr b/marchenko/demo/test2/backpropf2.scr
deleted file mode 100755
index 75535781bece263a61cf4256c3b1852d59a729dc..0000000000000000000000000000000000000000
--- a/marchenko/demo/test2/backpropf2.scr
+++ /dev/null
@@ -1,56 +0,0 @@
-#!/bin/bash
-
-export PATH=$HOME/src/OpenSource/bin/:$PATH:
-
-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
-
-
diff --git a/marchenko/demo/test2/clean b/marchenko/demo/test2/clean
deleted file mode 100755
index 3890128152ba3f4b11471dfdb5ddd1399840bc08..0000000000000000000000000000000000000000
--- a/marchenko/demo/test2/clean
+++ /dev/null
@@ -1,4 +0,0 @@
-#!/bin/bash
-
-rm *.su *.bin *.eps nep line* *.asci
-
diff --git a/marchenko/demo/test2/conv.gnp b/marchenko/demo/test2/conv.gnp
deleted file mode 100644
index 119341bef971d8e8dc3e7e4123c32a64f674f5d3..0000000000000000000000000000000000000000
--- a/marchenko/demo/test2/conv.gnp
+++ /dev/null
@@ -1,15 +0,0 @@
-set style data linespoints
-set mytics 10
-set xlabel 'number of iterations'
-set ylabel 'convergence rate'
-set size 2.0,2.0
-set size ratio 0.6
-set grid
-
-set log y
-set nolog x
-
-set term postscript eps font 'Helvetica,12' linewidth 4 fontscale 3
-set output 'convergence.eps'
-plot 'conv.txt' using 1:($2) lw 3 notitle
-
diff --git a/marchenko/demo/test2/conv.txt b/marchenko/demo/test2/conv.txt
deleted file mode 100644
index f2e99f36e05c42e47cbbae161bfad8327bae2530..0000000000000000000000000000000000000000
--- a/marchenko/demo/test2/conv.txt
+++ /dev/null
@@ -1,16 +0,0 @@
-0 1.000000e+00
-1 8.104102e-01
-2 2.776407e-01
-3 1.775258e-01
-4 1.278046e-01
-5 8.376110e-02
-6 6.221900e-02
-7 4.089906e-02
-8 3.275844e-02
-9 2.070254e-02
-10 1.920658e-02
-11 1.091778e-02
-12 1.282995e-02
-13 6.060715e-03
-14 9.706275e-03
-15 3.603180e-03
diff --git a/marchenko/demo/test2/epsBackprop.scr b/marchenko/demo/test2/epsBackprop.scr
deleted file mode 100755
index 5c2ecbcc92358b2cb0fe58a5914f54aabf1b4dc6..0000000000000000000000000000000000000000
--- a/marchenko/demo/test2/epsBackprop.scr
+++ /dev/null
@@ -1,68 +0,0 @@
-#!/bin/bash
-
-export PATH=$HOME/src/OpenSource/bin/:$PATH:
-
-# Add interface line to postscript file of model
-cat << EOF1 > line1
-400 -2500
-400 2500
-EOF1
-
-cat << EOF2 > line2
-700 -2500
-700 2500
-EOF2
-
-cat << EOF3 > line3
-1100 -2500
-1100 2500
-EOF3
-
-dx=5
-file_snap="backpropf2"
-dtsnap=0.01
-nsnap=101
-
-sumax < ${file_snap}_sp.su mode=abs outpar=nep
-clip=`cat nep | awk '{print $1/2}'`
-
-#first snap-shot with labels
-# fldr=71
-# 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 \
-# label1="depth (m)" label2="lateral distance (m)" \
-# x1beg=0 x1end=1250.0 clip=${clip} \
-# 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}_labels.eps
-
-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=${clip} \
- 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 99 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=${clip} \
- 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/test2/epsCompare.scr b/marchenko/demo/test2/epsCompare.scr
deleted file mode 100755
index e9dae68ee730f2cceb5b36018040ceca5b1f49fa..0000000000000000000000000000000000000000
--- a/marchenko/demo/test2/epsCompare.scr
+++ /dev/null
@@ -1,37 +0,0 @@
-#!/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 \
- label1="time (s)" label2="lateral distance (m)" \
- 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=#F \
- hbox=4 wbox=8 labelsize=10 linewidth=2.0 \
- label1="time (s)" label2="lateral distance (m)" \
- n1tic=2 d2=250 f1=0.0 x1beg=0 x1end=2.004 d1num=0.4 \
- f2=-2250 f2num=-2000 d2num=500 > ref.eps
-
-sed -i.old -e "s/%%EndProlog/[ 1 1 ] 0 setdash %%EndProlog/" green.eps
-sed -i.old -e "s/0.5 0.5 0.5 setrgbcolor/0.65 0.65 0.65 setrgbcolor /" ref.eps
-
-psmerge in=ref.eps in=green.eps > mergeGreenRef.eps
-
diff --git a/marchenko/demo/test2/epsIterwithLabels.scr b/marchenko/demo/test2/epsIterwithLabels.scr
deleted file mode 100755
index cfb5a6a0e6e0a2515986007b94d8ee53faa4a74d..0000000000000000000000000000000000000000
--- a/marchenko/demo/test2/epsIterwithLabels.scr
+++ /dev/null
@@ -1,76 +0,0 @@
-#!/bin/bash
-
-export PATH=$HOME/src/OpenSource/bin/:$PATH:
-
-#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/8}'`
-
-#set same clip factor for Green;s function updates
-file=pgreen_004.su
-file_base=${file%.su}
-sumax < $file mode=abs outpar=nep
-clipgreen=`cat nep | awk '{print $1/4}'`
-
-#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 \
- label1="time (s)" label2="lateral distance (m)" \
- 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}_labels.eps
-
-file=f1min_$piter.su
-file_base=${file%.su}
-sumax < $file mode=abs outpar=nep
-clipref=`cat nep | awk '{print $1/5}'`
-supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file\
- label1="time (s)" label2="lateral distance (m)" \
- n1tic=2 d2=5 x1beg=-1.504 x1end=1.5 d1num=0.4 \
- f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > ${file_base}_labels.eps
-
-file=f1plus_$piter.su
-file_base=${file%.su}
-sumax < $file mode=abs outpar=nep
-clipref=`cat nep | awk '{print $1/5}'`
-supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file\
- label1="time (s)" label2="lateral distance (m)" \
- n1tic=2 d2=5 x1beg=-1.504 x1end=1.5 d1num=0.4 \
- f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > ${file_base}_labels.eps
-
-file=pgreen_$piter.su
-file_base=${file%.su}
-sumax < $file mode=abs outpar=nep
-clipref=`cat nep | awk '{print $1/4}'`
-suwind key=gx min=-2250000 max=2250000 < $file | \
- supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 \
- label1="time (s)" label2="lateral distance (m)" \
- n1tic=2 d2=5 f1=0.0 x1beg=0 x1end=2.004 d1num=0.4 \
- f2=-2250 f2num=-2000 d2num=1000 clip=$clipgreen > ${file_base}_labels.eps
-
-done
-
-
-#special treatment of f1+ zero-iteration: which is zero, to make a nice gray plot (and not black)
-file=f1plus_001.su
-file_base=${file%.su}
-supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file\
- label1="time (s)" label2="lateral distance (m)" \
- n1tic=2 d2=5 x1beg=-1.504 x1end=1.5 d1num=0.4 \
- f2=-2250 f2num=-2000 d2num=1000 bclip=1 wclip=-1> ${file_base}_labels.eps
-
diff --git a/marchenko/demo/test2/epsMarchenkoIter.scr b/marchenko/demo/test2/epsMarchenkoIter.scr
deleted file mode 100755
index b2a417474810933105a76d01eb2c37168367ffda..0000000000000000000000000000000000000000
--- a/marchenko/demo/test2/epsMarchenkoIter.scr
+++ /dev/null
@@ -1,121 +0,0 @@
-#!/bin/bash
-
-export PATH=$HOME/src/OpenSource/bin/:$PATH:
-
-
-#Direct field of transmission repsponse
-file=p0plus.su
-file_base=${file%.su}
-sumax < $file mode=abs outpar=nep
-clipref=`cat nep | awk '{print $1/2}'`
-
-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 \
- label1="time (s)" label2="lateral distance (m)" \
- 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 \
- label1="time (s)" label2="lateral distance (m)" \
- 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/8}'`
-
-#set same clip factor for Green;s function updates
-file=pgreen_004.su
-file_base=${file%.su}
-sumax < $file mode=abs outpar=nep
-clipgreen=`cat nep | awk '{print $1/4}'`
-
-#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/5}'`
-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/5}'`
-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/4}'`
-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=$clipgreen > $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
-mg1=`cat nepmg | awk '{print $1}'`
-rf1=`cat neprf | awk '{print $1}'`
-value=${value/[eE][+][0]/*10^}
-mg=${mg1/[eE][+][0]/*10^}
-rf=${rf1/[eE][+][0]/*10^}
-rm nep*
-scale=$(echo "scale=3; ($rf)/($mg)" | bc -l)
-scale=2.0
-echo $scale
-
-(suwind key=gx min=0 max=0 < referenceP_rp.su; \
- suwind key=gx min=0 max=0 itmax=511 < pgreen_$piter.su | sugain scale=$scale ) | \
- supsgraph hbox=6 wbox=2 labelsize=10 linegray=0.5,0.0 style=seismic \
- lineon=1.0,1.0 lineoff=0.0,1.0 linewidth=1.0,1.0 x2beg=-$rf1 x2end=$rf1 > compare_$piter.eps
-
-done
-
-
-#special treatment of f1+ zero-iteration: which is zero, to make a nice gray plot (and not black)
-file=f1plus_001.su
-file_base=${file%.su}
-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 bclip=1 wclip=-1> $file_base.eps
-
diff --git a/marchenko/demo/test2/epsModel.scr b/marchenko/demo/test2/epsModel.scr
deleted file mode 100755
index 5ae0b460f468bf00cb8804d8882d6fa35a4f7885..0000000000000000000000000000000000000000
--- a/marchenko/demo/test2/epsModel.scr
+++ /dev/null
@@ -1,68 +0,0 @@
-#!/bin/bash
-
-export PATH=$HOME/src/OpenSource/bin/:$PATH:
-
-# Add interface line to postscript file of model
-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 lstyle=vertright legend=1 threecolor=0 \
- label1="depth (m)" label2="lateral distance (m)" \
- 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 lstyle=vertright legend=1 threecolor=0 \
- label1="depth (m)" label2="lateral distance (m)" \
- 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="frequency (1/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 \
- label1="time (s)" label2="lateral distance (m)" \
- 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 \
- label1="time (s)" label2="lateral distance (m)" \
- 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/test2/figAppendix.scr b/marchenko/demo/test2/figAppendix.scr
deleted file mode 100755
index 295d2cf53c10926bd222f8bc05b310763bf2f7d1..0000000000000000000000000000000000000000
--- a/marchenko/demo/test2/figAppendix.scr
+++ /dev/null
@@ -1,47 +0,0 @@
-#!/bin/bash
-
-export PATH=$HOME/src/OpenSource/bin/:$PATH:
-
-file=iter_002.su
-file_base=${file%.su}
-
-ns=`surange < $file | grep ns | awk '{print $2}'`
-dtrcv=`surange < $file | grep dt | awk '{print $2/1000000.0}'`
-shift=$(echo "scale=4; ($dtrcv*($ns/2.0-1))" | bc -l)
-suzero < $file itmax=$ns | suaddnoise | sushw key=f1 a=0 > noise.su
-file_base=noise
-sumax < ${file_base}.su mode=abs outpar=nep
-clipiter=`cat nep | awk '{print $1/6}'`
-clipref=$clipiter
-
-#basop choice=shift shift=$shift file_in=$file file_out=${file_base}_t0.su
-
-for above in 0 1 -1 2 4
-do
-fmute file_mute=iniFocus_rp.su file_shot=${file_base}.su file_out=nep.su above=${above} shift=8 verbose=1 check=1 hw=4
-
-basop choice=shift shift=-$shift file_in=nep.su file_out=nep_t0.su
-supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < nep.su \
- label1="time (s)" label2="lateral distance (m)" \
- n1tic=2 d2=5 x1beg=0 d1num=0.5 \
- 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}_above${above}.eps
-done
-
-for shift in 0 20 -20
-do
-fmute file_mute=iniFocus_rp.su file_shot=${file_base}.su file_out=nep.su above=${above} shift=$shift verbose=1 check=1 hw=4
-mv pslinepos.asci pslinepos${shift}.asci
-done
-
-suzero < $file itmax=$ns | sushw key=f1 a=0 > zero.su
-sumax < iniFocus_rp.su mode=abs outpar=nep
-clipiter=`cat nep | awk '{print $1/6}'`
-clipref=$clipiter
-supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < iniFocus_rp.su \
- label1="time (s)" label2="lateral distance (m)" \
- n1tic=2 d2=5 x1beg=0 d1num=0.5 \
- curve=pslinepos0.asci,pslinepos20.asci,pslinepos-20.asci npair=901,901,901 \
- curvewidth=1,1,1 curvecolor=white,black,black curvedash=3,3,3 \
- f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > iniFocus_shifts.eps
-
diff --git a/marchenko/demo/test2/initialFocus.scr b/marchenko/demo/test2/initialFocus.scr
deleted file mode 100755
index 6b1408fef5b71efe00987d2f892d83fa038aaef3..0000000000000000000000000000000000000000
--- a/marchenko/demo/test2/initialFocus.scr
+++ /dev/null
@@ -1,39 +0,0 @@
-#!/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=3000 dx=$dx dz=$dx cp0=2500 ro0=1000 \
- orig=-5000,0 file_base=modelup verbose=2 \
- intt=def x=-5000,5000 z=750,750 poly=0 cp=2500 ro=2000 \
- intt=def x=-5000,5000 z=1500,1500 poly=0 cp=2500 ro=1000
-
-makewave fp=25 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
-
-export OMP_NUM_THREADS=7
-
-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=10 \
- xrcv1=-4250 xrcv2=4250 \
- zrcv1=0 zrcv2=0 \
- xsrc=0 zsrc=2200 dxshot=20 nshot=5 \
- ntaper=101 \
- left=2 right=2 top=2 bottom=2
-
diff --git a/marchenko/demo/test2/marchenko.scr b/marchenko/demo/test2/marchenko.scr
deleted file mode 100755
index a6ddbffbcef776c38b7466ad19ce9f3c0172007a..0000000000000000000000000000000000000000
--- a/marchenko/demo/test2/marchenko.scr
+++ /dev/null
@@ -1,45 +0,0 @@
-#!/bin/bash -x
-
-export PATH=$HOME/src/OpenSource/bin:$PATH:
-export OMP_NUM_THREADS=8
-
-#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=8
-
-#marchenko file_shot=shotsdx5_rp.su file_tinv=p0plus.su verbose=2 \
-#marchenko file_shot=Joeri/shots.su file_tinv=p0plus.su verbose=2 \
-#apply the Marchenko algorithm
-#marchenko file_shot=Joeri/shots.su file_tinv=p0plus.su verbose=2 \
-marchenko file_shot=Joeri/shots.su file_tinv=Joeri/nep.su verbose=2 \
- tap=0 niter=8 hw=8 shift=12 smooth=3 scale=1 \
- file_green=pgreen.su file_gplus=Gplus0.su file_gmin=Gmin0.su \
- file_f1plus=f1plus.su file_f1min=f1min0.su file_f2=f2.su
-
-exit
-
-#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
-
-# plot for convergence rate, the values in conv.txt are collected from the output of the marhenko program with verbose=2
-# marchenko: - iSyn 0: Ni at iteration 0 has energy 6.234892e+02; relative to N0 1.000000e+00
-#a2b < conv.txt | \
-#psgraph n=16 style=normal hbox=2 wbox=6 labelsize=10 \
-#label2='convergence rate' label1='iteration number' > convergence.eps
-
-# If guplot is installed: the same plot can also be produced by gnuplot this figure is used in the paper
-#gnuplot conv.gnp
diff --git a/marchenko/demo/test2/marchenkoIter.scr b/marchenko/demo/test2/marchenkoIter.scr
deleted file mode 100755
index 401f97f7c2108e92e0ff5ca813d9fdfd2b4d183a..0000000000000000000000000000000000000000
--- a/marchenko/demo/test2/marchenkoIter.scr
+++ /dev/null
@@ -1,21 +0,0 @@
-#!/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/test2/model.scr b/marchenko/demo/test2/model.scr
deleted file mode 100755
index 57b76e5ee3efa2449e8f293bdd86653b7f41ca55..0000000000000000000000000000000000000000
--- a/marchenko/demo/test2/model.scr
+++ /dev/null
@@ -1,78 +0,0 @@
-#!/bin/bash
-
-#adjust this PATH to where the code is installed
-export PATH=$HOME/src/OpenSource/bin:$PATH:
-
-dx=2.5
-dt=0.0005
-
-makemod sizex=20000 sizez=3000 dx=$dx dz=$dx cp0=2500 ro0=1000 \
- orig=-10000,0 file_base=model10.su verbose=2 \
- intt=def x=-10000,10000 z=750,750 poly=0 cp=2500 ro=2000 \
- intt=def x=-10000,10000 z=1500,1500 poly=0 cp=2500 ro=1000 \
- intt=def x=-10000,10000 z=2375,2375 poly=0 cp=2500 ro=2000
-
-#define gridded model for FD computations
-
-#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=20
-
-#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=10.0 \
- xrcv1=-9000 xrcv2=9000 \
- zrcv1=0 zrcv2=0 \
- xsrc=0 zsrc=0 \
- npml=101 \
- left=2 right=2 top=2 bottom=2
-
-#define homogenoeus model to compute direct wave only
-makemod sizex=20000 sizez=3000 dx=$dx dz=$dx cp0=2500 ro0=1000 \
- orig=-10000,0 file_base=hom 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=10.0 \
- xrcv1=-9000 xrcv2=9000 \
- zrcv1=0 zrcv2=0 \
- xsrc=0 zsrc=0 \
- npml=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/test2/p5all.scr b/marchenko/demo/test2/p5all.scr
deleted file mode 100755
index d83fd2e270a26e19b0617fae2a553245ffb4052c..0000000000000000000000000000000000000000
--- a/marchenko/demo/test2/p5all.scr
+++ /dev/null
@@ -1,30 +0,0 @@
-#!/bin/bash -x
-
-export PATH=$HOME/src/OpenSource/bin:$PATH:
-
-# Generate the full R matrix for a fixed spread geometry.
-
-dxshot=10000 # with scalco factor of 1000
-ishot=0
-nshots=851
-
-echo $1
-
-rm shotsdx5_rp.su
-
-while (( ishot < nshots ))
-do
-
- (( xsrc = -4250000 + ${ishot}*${dxshot} ))
- (( tr1 = 901 - ${ishot} ))
- (( tr2 = ${tr1} + 850 ))
- 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,-4250000,$ishot,851 b=0,10000,0,0 j=0,851,0,0 >> shotsdx5_rp.su
-
-done
-
diff --git a/marchenko/demo/test2/referenceShot.scr b/marchenko/demo/test2/referenceShot.scr
deleted file mode 100755
index b7a2b771341b3115d71bfebe2ec06e308846cbc6..0000000000000000000000000000000000000000
--- a/marchenko/demo/test2/referenceShot.scr
+++ /dev/null
@@ -1,34 +0,0 @@
-#!/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/marchenko.c b/marchenko/marchenko.c
index f62852d4d3bd1d151e49da082360f5560dd65046..a0fc3e328ed27de41706943c00a31f3831e75dd7 100644
--- a/marchenko/marchenko.c
+++ b/marchenko/marchenko.c
@@ -839,7 +839,7 @@ nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int np
if (verbose>4) vmess("*** Shot gather %d processed ***", k);
} /* end of nshots (k) loop */
- } /* end of if reci
+ } /* end of if reci */
/* if reciprocal traces are enabled start a new loop over reciprocal shot positions */
if (reci != 0) {
diff --git a/raytime/JespersRayTracer.c b/raytime/JespersRayTracer.c
index ba6ae35845c1e6260999a0799339a8cc09dd66b5..255c4d87e1bb00b44b5a45e3c255b38fed4336ea 100644
--- a/raytime/JespersRayTracer.c
+++ b/raytime/JespersRayTracer.c
@@ -16,7 +16,7 @@
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
-static float H, L, W;
+static float H, L, W, iH, iL, iW;
typedef struct _icoord { /* 3D coordinate integer */
int z;
@@ -40,12 +40,12 @@ int yPointIndex(const float _y, int ny, float W);
fcoord getSlownessGradient(const float _x, const float _z, float *slowness, icoord size);
float qMulGradU1(const float _x, const float _z, const float _angle, float *slowness, icoord size);
float greenTwoP(const float _so, const float _slow, const float _sL, int nRay, fcoord s, fcoord r, float *slowness, icoord size);
-float qatso(const float _so, const float _angle, int nRay, fcoord s, fcoord r, fcoord *rayReference3D, float *slowness, icoord size);
+float qatso(const float _so, const float _angle, int nRay, fcoord s, fcoord r, fcoord *rayReference3D, float *slowness, icoord size, float uo);
float slownessA(float *slowness, icoord size, float _x, float _y, float _z);
-float getdT2(const float _x, const float _z, const float so, const float _angle, const float _ds, int nRay, fcoord s, fcoord r, fcoord *rayReference3D, float *slowness, icoord size);
+float getdT2(const float _x, const float _z, const float so, const float _angle, const float _ds, int nRay, fcoord s, fcoord r, fcoord *rayReference3D, float *slowness, icoord size, float uo);
float greenIntP(const float _so, const float _s, const float _sL, float *slowness, icoord size, int nRay, fcoord r, fcoord s);
float secondDerivativeU1(float *slowness, icoord size, const float _x, const float _z, const float _angle, fcoord s, fcoord r);
-int calculatePerturbedRay(fcoord *rayPerturbed3D, fcoord s, fcoord r, int nRay, fcoord *rayReference3D, float *slowness, icoord size);
+int calculatePerturbedRay(fcoord *rayPerturbed3D, fcoord s, fcoord r, int nRay, fcoord *rayReference3D, float *slowness, icoord size, float uo);
float angle2qx(const float _angle);
float angle2qz(const float _angle);
float ModelInterpolation_slowness2D(float *slowness, icoord size, const float _x, const float _z);
@@ -96,7 +96,7 @@ int getWaveParameter(float *slowness, icoord size, float dgrid, fcoord s, fcoord
return(-1);
uo = referenceSlowness(slowness, size, nRayTmp, r, s);
-
+
T0 = lengthRefRay*uo;
ds = lengthRefRay/(nRayTmp-1);
J = lengthRefRay;
@@ -121,14 +121,15 @@ int getWaveParameter(float *slowness, icoord size, float dgrid, fcoord s, fcoord
so = i*ds;
if (ray.useT2 != 0)
- T2 += getdT2(x, z, so, angle, ds, nRayTmp, s, r, rayReference3D, slowness, size);
+ T2 += getdT2(x, z, so, angle, ds, nRayTmp, s, r, rayReference3D, slowness, size, uo);
/*if (ray.geomspread != 0) {
if (so <= 0) {
dQdPhi = 0;
}
else {
- greentmp = greenIntP(lengthRefRay, so, lengthRefRay, slowness, size, nRayTmp, r, s);
+ greentmp = 0;
+ if (so <= lengthRefRay) greentmp = (lengthRefRay - so)/uo;
dQdPhi += greentmp*secondDerivativeU1(slowness, size, x, z, angle, r, s)*ds/so;
}
}*/
@@ -176,6 +177,10 @@ int getnRay(icoord size, fcoord s, fcoord r, float dx, int nRayStep)
L = (size.x-1)*dx;
W = (size.y-1)*dx;
+ if (H!=0.0) iH = 1.0/H;
+ if (L!=0.0) iL = 1.0/L;
+ if (W!=0.0) iW = 1.0/W;
+
if (size.y == 1) { // 2D model
dn = (size.x + size.z)/2;
dl = sqrt(pow(L, 2) + pow(H, 2))/dn;
@@ -216,7 +221,7 @@ int traceTwoPoint(fcoord s, fcoord r, int nRay, fcoord *rayReference3D)
}
-int calculatePerturbedRay(fcoord *rayPerturbed3D, fcoord s, fcoord r, int nRay, fcoord *rayReference3D, float *slowness, icoord size)
+int calculatePerturbedRay(fcoord *rayPerturbed3D, fcoord s, fcoord r, int nRay, fcoord *rayReference3D, float *slowness, icoord size, float uo)
{
float si, sl, deltaS, gso, angle, qx, qz;
int i;
@@ -236,7 +241,7 @@ int calculatePerturbedRay(fcoord *rayPerturbed3D, fcoord s, fcoord r, int nRay,
{
si = i*deltaS;
- gso = qatso(si, angle, nRay, s, r, rayReference3D, slowness, size);
+ gso = qatso(si, angle, nRay, s, r, rayReference3D, slowness, size, uo);
rayPerturbed3D[i].x = rayReference3D[i].x + qx*gso;
rayPerturbed3D[i].z = rayReference3D[i].z + qz*gso;
@@ -299,11 +304,16 @@ float angle2qz(const float _angle)
// Sofar used in 2D only
-float qatso(const float _so, const float _angle, int nRay, fcoord s, fcoord r, fcoord *rayReference3D, float *slowness, icoord size)
+float qatso(const float _so, const float _angle, int nRay, fcoord s, fcoord r, fcoord *rayReference3D, float *slowness, icoord size, float uo)
{
float slow, sl, deltaS, x, z;
float qatsol;
+ float greenTwoP = 0;
int i;
+ float qMulGradU1;
+ fcoord slownessGradient;
+ float gradu1x, gradu1z;
+ float qx, qz;
sl = sqrt(pow((r.x-s.x),2) + pow((r.z-s.z),2) + pow((r.y-s.y),2));
@@ -313,22 +323,49 @@ float qatso(const float _so, const float _angle, int nRay, fcoord s, fcoord r, f
}
deltaS = sl/(nRay-1);
-
+// uo = referenceSlowness(slowness, size, nRay, r, s);
+
qatsol = 0;
for (i = 0; i < nRay; i++)
{
slow = i*deltaS;
x = rayReference3D[i].x;
z = rayReference3D[i].z;
-// fprintf(stderr,"qatso: calling greenTwoP for iray %d (/%d)\n",i,nRay);
+
+ if (slow <= _so)
+ greenTwoP = -(1 - _so/sl)*slow/uo;
+ else
+ greenTwoP = -_so*(1-slow/sl)/uo;
+
+ slownessGradient = getSlownessGradient(x, z, slowness, size);
+ gradu1x = slownessGradient.x;
+ gradu1z = slownessGradient.z;
+
+ if ((_angle >= 0) && (_angle < PI/2)) {
+ qx = -cos(_angle);
+ qz = sin(_angle);
+ }
+ else if ((_angle >= PI/2) && (_angle < PI)) {
+ qx = sin(_angle - PI/2);
+ qz = cos(_angle - PI/2);
+ }
+ else if ((_angle >= PI) && (_angle < 3*PI/2)) {
+ qx = cos(_angle - PI);
+ qz = -sin(_angle - PI);
+ }
+ else if ((_angle >= 3*PI/2) && (_angle <= 2*PI)) {
+ qx = -sin(_angle - 3*PI/2);
+ qz = -cos(_angle - 3*PI/2);
+ }
- qatsol += greenTwoP(_so, slow, sl, nRay, s, r, slowness, size)*qMulGradU1(x, z, _angle, slowness, size)*deltaS;
+ qMulGradU1 = qx*gradu1x + qz*gradu1z;
+ qatsol += greenTwoP*qMulGradU1*deltaS;
}
return(qatsol);
}
-float getdT2(const float _x, const float _z, const float _so, const float _angle, const float _ds, int nRay, fcoord s, fcoord r, fcoord *rayReference3D, float *slowness, icoord size)
+float getdT2(const float _x, const float _z, const float _so, const float _angle, const float _ds, int nRay, fcoord s, fcoord r, fcoord *rayReference3D, float *slowness, icoord size, float uo)
{
float T2 = 0;
float qatsol;
@@ -336,7 +373,7 @@ float getdT2(const float _x, const float _z, const float _so, const float _angle
// fprintf(stderr,"getdT2: calling qatso nRay=%d\n",nRay);
- qatsol = qatso(_so, _angle, nRay, s, r, rayReference3D, slowness, size);
+ qatsol = qatso(_so, _angle, nRay, s, r, rayReference3D, slowness, size, uo);
// fprintf(stderr,"getdT2: calling qMulGradU1\n");
@@ -490,7 +527,7 @@ int xPointIndex(const float _x, int nx, float L)
else
{
if (0 < L)
- i = _x*nx/L;
+ i = _x*nx*iL;
else
i = 0;
}
@@ -509,7 +546,7 @@ int zPointIndex(const float _z, int nz, float H)
else
{
if (0 < H)
- i = _z*nz/H;
+ i = _z*nz*iH;
else
i = 0;
}
@@ -529,7 +566,7 @@ int yPointIndex(const float _y, int ny, float W)
else
{
if (0 < W)
- i = ny*(_y/W + 0.5);
+ i = ny*(_y*iW + 0.5);
else
i = 0;
}
@@ -551,14 +588,11 @@ float ModelInterpolation_slowness2D(float *slowness, icoord size, const float _x
nx = size.x;
nz = size.z;
- ixCoordinate = (int) _x*nx/L;
-
- if (ixCoordinate >= nx)
- ixCoordinate = nx;
+ ixCoordinate = (int)(_x*nx)*iL;
- if (ixCoordinate == nx)
+ if (ixCoordinate >= nx)
{
- x1 = (float) L*(ixCoordinate-1)/nx;
+ x1 = (float) L*(nx-1)/nx;
x2 = (float) L;
}
else if (ixCoordinate <= 0)
@@ -572,26 +606,11 @@ float ModelInterpolation_slowness2D(float *slowness, icoord size, const float _x
x2 = (float) L*(ixCoordinate+1)/nx;
}
- if (x1 < 0)
- x1 = 0;
-
- if (x1 > L)
- x1 = L;
-
- if (x2 < 0)
- x2 = 0;
-
- if (x2 > L)
- x2 = L;
-
- izCoordinate = (int) _z*nz/H;
+ izCoordinate = (int) _z*nz*iH;
if (izCoordinate >= nz)
- izCoordinate = nz;
-
- if (izCoordinate == nz)
{
- z1 = (float) H*(izCoordinate-1)/nz;
+ z1 = (float) H*(nz-1)/nz;
z2 = (float) H;
}
else if (izCoordinate <= 0)
@@ -605,20 +624,8 @@ float ModelInterpolation_slowness2D(float *slowness, icoord size, const float _x
z2 = (float) H*(izCoordinate+1)/nz;
}
- if (z1 < 0)
- z1 = 0;
-
- if (z1 > H)
- z1 = H;
-
- if (z2 < 0)
- z2 = 0;
-
- if (z2 > H)
- z2 = H;
-
- ix = xPointIndex(_x, size.x, L);
- iz = zPointIndex(_z, size.z, H);
+ ix = xPointIndex(_x, nx, L);
+ iz = zPointIndex(_z, nz, H);
if (ix == 0)
{
@@ -685,7 +692,7 @@ float ModelInterpolation_slowness3D(float *slowness, icoord size, const float _x
slow = f111 = f112 = f212 = f211 = f121 = f122 = f222 = f221 = 0;
- ixCoordinate = _x*nx/L;
+ ixCoordinate = _x*nx*iL;
if (ixCoordinate >= nx)
ixCoordinate = nx;
@@ -718,7 +725,7 @@ float ModelInterpolation_slowness3D(float *slowness, icoord size, const float _x
if (x2 > L)
x2 = L;
- izCoordinate = _z*nz/H;
+ izCoordinate = _z*nz*iH;
if (izCoordinate >= nz)
izCoordinate = nz;
@@ -751,7 +758,7 @@ float ModelInterpolation_slowness3D(float *slowness, icoord size, const float _x
if (z2 > H)
z2 = H;
- iyCoordinate = ny*(_y/W + 0.5);
+ iyCoordinate = ny*(_y*iW + 0.5);
if (iyCoordinate >= ny)
iyCoordinate = ny;
diff --git a/raytime/model.scr b/raytime/model.scr
index ac8479bb19953493dd8feb68d5b4bc8deb0adbe4..be518c400ab16df63d5cce36b31a21a6a7b0e3c0 100755
--- a/raytime/model.scr
+++ b/raytime/model.scr
@@ -23,7 +23,9 @@ export OMP_NUM_THREADS=4
zrcv1=0 zrcv2=0 \
xsrc=0 zsrc=1100 \
nxshot=1 dxshot=10 \
- nzshot=1 dzshot=10
+ nzshot=1 dzshot=10 useT2=1
+
+exit
./raytime \
file_cp=syncl_cp.su \