diff --git a/fdelmodc/applySource.c b/fdelmodc/applySource.c
index 74d19650569cef760b5b5c1649195cbdba25c69e..8e4574f299c49a15ad172dee6e9273edd1de81ae 100644
--- a/fdelmodc/applySource.c
+++ b/fdelmodc/applySource.c
@@ -114,6 +114,7 @@ int applySource(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int i
if (verbose>=4 && itime==0) {
vmess("Source %d positioned at grid ix=%d iz=%d",isrc, ix, iz);
}
+
/* cosine squared windowing to reduce edge effects on shot arrays */
if ( (src.n>1) && src.window) {
scl = 1.0;
@@ -139,7 +140,6 @@ int applySource(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int i
vmess("Source %d at grid [ix=%d,iz=%d] at itime %d has value %e",isrc, ix,iz, itime, src_ampl);
}
-
/* Force source */
if (src.type == 6) {
diff --git a/fdelmodc/demo/matlab/FD_mod_grid.m b/fdelmodc/demo/matlab/FD_mod_grid.m
new file mode 100644
index 0000000000000000000000000000000000000000..fcc86b9cc11f109aa05c9c5a5ae4f0c63a378c45
--- /dev/null
+++ b/fdelmodc/demo/matlab/FD_mod_grid.m
@@ -0,0 +1,138 @@
+function [ Ptsct, Pfinc, Pfsct, f_out] = FD_mod_grid( xS, xR, ntrcv, dtrcv, dx, cgrid, rhogrid, f )
+%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=-(dims(1)-1)*dx/2
+origx=-(dims(2)-1)*dx/2
+zrcv1=xR(1)-origz;
+zrcv2=xR(1)-origz;
+xrcv1=0.0;
+xrcv2=xrcv1+(dims(2)-1)*dx;
+zsrc=xS(1)-origz;
+xsrc=xS(2)-origx;
+tmod=(ntrcv-1)*dtrcv;
+
+%compute dt for modeling dt < 0.606*dx/Cmax
+Cmax=max(cgrid(:));
+dxmod=dx;
+dtmax=0.606*dxmod/Cmax;
+dtmod=dtrcv/(ceil(dtrcv/dtmax));
+ntwave=16384;
+ntfft=4*ntrcv;
+
+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,'which surange\n');
+fprintf(fileID,'set -x\n');
+
+fprintf(fileID,'df=0.25\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=d1,d2,gx,scalco a=%f,%f,%d,-1000 b=0,0,%d,0 > mod_ro.su\n',dx, dx, int32(xrcv1*1000), int32(dx*1000));
+fprintf(fileID,'suaddhead < mod_cp.bin ntrpr=%d ns=%d | \\\n',dims(2), dims(1));
+fprintf(fileID,'sushw key=d1,d2,gx,scalco a=%f,%f,%d,-1000 b=0,0,%d,0 > mod_cp.su\n',dx, dx, int32(xrcv1*1000), int32(dx*1000));
+fprintf(fileID,'makewave w=fw fmin=0 flef=6 frig=94 fmax=100 dt=$dt file_out=wavefw.su nt=%d t0=0.4 scale=0 scfft=1 verbose=1\n', ntwave);
+fprintf(fileID,'makewave w=fw fmin=0 flef=6 frig=94 fmax=100 dt=%f file_out=wavefwdt.su nt=%d t0=0.4 scale=0 scfft=1 verbose=1\n', dtrcv, ntfft);
+fprintf(fileID,'export filecp=mod_cp.su\n');
+fprintf(fileID,'export filero=mod_ro.su\n');
+fprintf(fileID,'export OMP_NUM_THREADS=2\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,'file_rcv=recvgrid.su \\\n');
+fprintf(fileID,'rec_type_vz=0 rec_type_p=1 rec_int_vz=2 \\\n');
+fprintf(fileID,'xrcv1=%d xrcv2=%d zrcv1=%d zrcv2=%d \\\n',xrcv1,xrcv2,zrcv1,zrcv2);
+fprintf(fileID,'dxrcv=%d \\\n', dx);
+fprintf(fileID,'dtrcv=%e \\\n', dtrcv);
+fprintf(fileID,'xsrc=%d zsrc=%d\\\n', xsrc, zsrc);
+fprintf(fileID,'src_type=1 tmod=%e \\\n', tmod);
+fprintf(fileID,'ntaper=100 \\\n');
+fprintf(fileID,'left=2 right=2 bottom=2 top=2\n\n');
+fprintf(fileID,'\n');
+fprintf(fileID,'makemod file_base=hom.su \\\n');
+fprintf(fileID,'cp0=1500 ro0=3000 sizex=%d sizez=%d dx=%.1f dz=%.1f orig=0,0 verbose=1\n', sizex,sizez, dxmod,dxmod);
+fprintf(fileID,'export filecp=hom_cp.su\n');
+fprintf(fileID,'export filero=hom_ro.su\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,'file_rcv=recvhom.su \\\n');
+fprintf(fileID,'rec_type_vz=0 rec_type_p=1 \\\n');
+fprintf(fileID,'xrcv1=%d xrcv2=%d zrcv1=%d zrcv2=%d \\\n',xrcv1,xrcv2,zrcv1,zrcv2);
+fprintf(fileID,'dxrcv=%d \\\n', dx);
+fprintf(fileID,'dtrcv=%e \\\n', dtrcv);
+fprintf(fileID,'xsrc=%d zsrc=%d\\\n', xsrc, zsrc);
+fprintf(fileID,'src_type=1 tmod=%e \\\n', tmod);
+fprintf(fileID,'ntaper=100 \\\n');
+fprintf(fileID,'left=2 right=2 bottom=2 top=2\n\n');
+fprintf(fileID,'suop2 recvgrid_rp.su recvhom_rp.su op=diff > recv_rp.su\n');
+fprintf(fileID,'sustrip < recv_rp.su > recv_rp.bin\n');
+fprintf(fileID,'sustrip < recvhom_rp.su > recvhom_rp.bin\n');
+fprintf(fileID,'sustrip < wavefwdt.su > wavefwdt.bin\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);
+
+% Scattered field
+ns=ntrcv;
+ntr=dims(2);
+file='recv_rp.bin';
+fid=fopen(file,'r');
+Ptsct=fread(fid,[ns,ntr],'float32');
+fclose(fid);
+
+% Direct field
+file='recvhom_rp.bin';
+fid=fopen(file,'r');
+Ptinc=fread(fid,[ns,ntr],'float32');
+fclose(fid);
+
+df=double(1.0/(ntfft*dtrcv));
+a=round(f/df)+1;
+f_out=(a-1)*df
+fprintf('selected discrete frequency P %f\n', (a-1)*df)
+% f2=a*df %these are the selected discrete frequencies
+Pfft=fft(Ptsct,ntfft);
+Pfsct=Pfft(a,:);
+%
+Pfft=fft(Ptinc,ntfft);
+Pfinc=Pfft(a,:);
+
+ns=ntfft;
+ntr=1;
+file='wavefwdt.bin';
+fid=fopen(file,'r');
+Wt=fread(fid,[ns,ntr],'float32');
+fclose(fid);
+
+df=double(1.0/(ntfft*dtrcv));
+c=round(f/df)+1; % select frequencies as close as the one's in Pf
+fprintf('selected discrete frequency W %f\n', (c-1)*df)
+Wfft=fft(Wt,ntfft);
+Wf=Wfft(c,1);
+
+Pfsct = (Pfsct)./(1.0*Wf); % deconvolve for the wavelet
+Pfinc = (Pfinc)./(1.0*Wf); % deconvolve for the wavelet
+
+end
+
diff --git a/fdelmodc/demo/matlab/ForwardCircle.m b/fdelmodc/demo/matlab/ForwardCircle.m
new file mode 100644
index 0000000000000000000000000000000000000000..dd67cc7d51751acbcfebbc2e76f102e621a2a4c4
--- /dev/null
+++ b/fdelmodc/demo/matlab/ForwardCircle.m
@@ -0,0 +1,152 @@
+function [P_inc,P_sct,xR] = ForwardCircle( xS, xR, gsize, dxR, c_0, c_s, rho_0, rho_s, f )
+
+% xS = source position
+% dxR = receiver grid size
+% c_0 = wave speed in embedding
+% c_s = wave speed of scattering object
+% rho_0 = mass density in embedding
+% rho_s = mass density of scattering object
+% f = temporal frequency;
+
+wavelength = c_0 / f; % wavelength
+s = 1e-16 + 1i*2*pi*f; % LaPlace parameter
+gam_0 = s/c_0; % propagation coefficient 0
+gam_s = s/c_s; % propagation coefficient s
+a = 40; % radius circle cylinder
+depth = 20; % distance betwen circle and origin
+
+disp(['wavelength = ' num2str(wavelength)]);
+
+% Thorbecke coordinates
+ xSource1 = xS(1); % source position
+ xSource2 = xS(2);
+ xReceiver1 = xR(1); % receiver positions
+ nr=(gsize(2)-1)/2;
+ xReceiver2 = (-nr:nr) * dxR;
+ NR = size(xReceiver2,2);
+
+% Van den Berg coordinates
+ xS(1) = xSource1;
+ xS(2) = xSource2;
+
+ xR(1,1:NR) = xReceiver1;
+ xR(2,1:NR) = xReceiver2(1:NR);
+
+% Compute 2D incident field ---------------------------------------------
+ DIS = sqrt( (xR(1,:)-xS(1)).^2 + (xR(2,:)-xS(2)).^2 );
+ p_inc = 1/(2*pi).* besselk(0,gam_0*DIS);
+
+ % multiply by rho_0 and plot
+ P_inc = rho_0 * p_inc ; % Take care of factor \rho
+
+ % Make grid
+ N1 = gsize(1); % number of samples in x_1
+ N2 = gsize(2); % number of samples in x_2
+ dx = dxR; % with meshsize dx
+ x1 = -(N1+1)*dx/2 + (1:N1)*dx;
+ x2 = -(N2+1)*dx/2 + (1:N2)*dx;
+ [X1,X2] = ndgrid(x1,x2);
+ % Now array subscripts are equivalent with Cartesian coordinates
+ % x1 axis points downwards and x2 axis is in horizontal direction
+ % x1 = X1(:,1) is a column vector in vertical direction
+ % x2 = X2(1,:) is a row vector in horizontal direction
+ R = sqrt(X1.^2 + X2.^2);
+ cgrid = c_s * (R < a) + c_0 * (R >= a);
+ rhogrid = rho_s * (R < a) + rho_0 * (R >= a);
+
+ x1 = X1(:,1); x2 = X2(1,:);
+ set(figure(8),'Units','centimeters','Position',[5 5 18 12]);
+ subplot(1,2,1)
+ imagesc(x2,x1,cgrid);
+ title('\fontsize{13} c-grid');
+ xlabel('x_1 \rightarrow');
+ ylabel('\leftarrow x_3');
+ axis('equal','tight');
+ colorbar('hor'); colormap jet;
+ subplot(1,2,2)
+ imagesc(x2,x1,rhogrid);
+ title('\fontsize{13} \rho-grid');
+ xlabel('x_1 \rightarrow');
+ ylabel('\leftarrow x_3');
+ axis('equal','tight');
+ colorbar('hor'); colormap jet;
+
+
+
+
+
+%-------------------------------------------------------------------------
+% CHECK EXACT INCIDENT FIELD AND BESSEL FUNCTION REPRESENTATION for r = a
+%-------------------------------------------------------------------------
+
+% Transform Cartesian coordinates to polar ccordinates
+ rS = sqrt(xS(1)^2+xS(2)^2); phiS = atan2(xS(2),xS(1));
+ r = a; phi = 0:.01:2*pi;
+
+% (1) Compute incident wave in closed form --------------------------------
+ DIS = sqrt(rS^2 + r.^2 - 2*rS*r.*cos(phiS-phi));
+ p_inc_exact = 1/(2*pi) .* besselk(0,gam_0*DIS);
+ dp_inc_exact = - gam_0 *(r-rS*cos(phiS-phi))./DIS ...
+ .* 1/(2*pi) .* besselk(1,gam_0*DIS);
+
+% (2) Compute incident wave as Bessel series with M+1terms --------------
+ M = 100; % increase M for more accuracy
+
+ p_inc = besselk(0,gam_0*rS) .* besseli(0,gam_0*r);
+ dp_inc = gam_0 * besselk(0,gam_0*rS) .* besseli(1,gam_0*r);
+ for m = 1 : M;
+ Ib0 = besseli(m,gam_0*r);
+ dIb0 = gam_0 * (besseli(m+1,gam_0*r) + m/(gam_0*r) * Ib0);
+ p_inc = p_inc + 2 * besselk(m,gam_0*rS) * Ib0 .* cos(m*(phiS-phi));
+ dp_inc = dp_inc + 2 * besselk(m,gam_0*rS) .* dIb0 .* cos(m*(phiS-phi));
+ end % m_loop
+ p_inc = 1/(2*pi) * p_inc;
+ dp_inc = 1/(2*pi) * dp_inc;
+
+% (3) Determine mean error and plot error in domain -----------------------
+ Error_p = p_inc - p_inc_exact;
+ disp(['normalized norm of error = ' ...
+ num2str(norm(Error_p(:),1)/norm(p_inc_exact(:),1))]);
+ Error_dp = dp_inc - dp_inc_exact;
+ disp(['normalized norm of error = ' ...
+ num2str(norm(Error_dp(:),1)/norm(dp_inc_exact(:),1))]);
+
+ set(figure(9),'Units','centimeters','Position',[5 5 18 14]);
+ subplot(2,1,1)
+ angle = phi * 180 / pi;
+ semilogy(angle,abs(Error_p)./abs(p_inc_exact)); axis tight;
+ xlabel('observation angle in degrees \rightarrow');
+ ylabel('abs(p_{inc}-p_{inc}^{exact}) / abs(p_{inc}^{exact}) \rightarrow');
+ subplot(2,1,2)
+ semilogy(angle,abs(Error_dp)./abs(dp_inc_exact)); axis tight;
+ title('\fontsize{12} relative error on circle boundary');
+ xlabel('observation angle in degrees \rightarrow');
+ ylabel('abs(dp_{inc}-dp_{inc}^{exact}) / abs(dp_{inc}^{exact}) \rightarrow');
+
+%--------------------------------------------------------------------------
+% COMPUTE SCATTERED FIELD WITH BESSEL FUNCTION REPRESENTATION for r > a
+%--------------------------------------------------------------------------
+
+% (4) Compute coefficients of series expansion ----------------------------
+ Z_0 = c_0 * rho_0; Z_s = c_s * rho_s;
+ arg0 = gam_0 * a; args = gam_s *a;
+ A = zeros(1,M+1);
+ for m = 0 : M;
+ Ib0 = besseli(m,arg0); dIb0 = besseli(m+1,arg0) + m/arg0 * Ib0;
+ Ibs = besseli(m,args); dIbs = besseli(m+1,args) + m/args * Ibs;
+ Kb0 = besselk(m,arg0); dKb0 = -besselk(m+1,arg0) + m/arg0 * Kb0;
+ A(m+1) = - ((1/Z_s) * dIbs*Ib0 - (1/Z_0) * dIb0*Ibs) ...
+ /((1/Z_s) * dIbs*Kb0 - (1/Z_0) * dKb0*Ibs);
+ end
+
+% (5) Compute scattered field at receivers (data) -------------------------
+ rR = sqrt(xR(1,:).^2 + xR(2,:).^2); phiR = atan2(xR(2,:),xR(1,:));
+ rS = sqrt(xS(1)^2 + xS(2)^2); phiS = atan2(xS(2),xS(1));
+ p_sct = A(1) * besselk(0,gam_0*rS).* besselk(0,gam_0*rR);
+ for m = 1 : M;
+ factor = 2 * besselk(m,gam_0*rS) .* cos(m*(phiS-phiR));
+ p_sct = p_sct + A(m+1) * factor .* besselk(m,gam_0*rR);
+ end % m_loop
+ p_sct = 1/(2*pi) * p_sct;
+ P_sct = rho_0 * p_sct; % Take care of factor \rho
+end
\ No newline at end of file
diff --git a/fdelmodc/demo/matlab/comparison.m b/fdelmodc/demo/matlab/comparison.m
new file mode 100644
index 0000000000000000000000000000000000000000..796e68149c957dade92e5c00dc3dce0f20e6e361
--- /dev/null
+++ b/fdelmodc/demo/matlab/comparison.m
@@ -0,0 +1,83 @@
+clear all; close all; clc;
+
+display('Running test');
+
+xS = [-100,0]; % source position: 100 m above center
+xR = [-60,0]; % central point of receiver array (-50*dxR:50*dxR)
+c_0 = 1500; % wave speed in embedding
+c_s = 3000; % wave speed in scattering object
+rho_0 = 3000; % mass density of enbedding
+rho_s = 1500; % mass density of scattering object
+
+dxR = 0.5; % receiver grid size
+gsize = [1+240/dxR,1+200/dxR]; % gridsize in [z,x], center of model at coordinate (0,0)
+f_in = 50; % selected frequency to compare, returned in Pf
+ntrcv = 256; % number of time samples in FD modeling
+dtrcv = 0.004; % dt in receivers
+
+
+% Make grid
+a = 40; % radius circle cylinder
+N1 = gsize(1); % number of samples in x_1
+N2 = gsize(2); % number of samples in x_2
+dx = dxR; % with meshsize dx
+x1 = -(N1+1)*dx/2 + (1:N1)*dx;
+x2 = -(N2+1)*dx/2 + (1:N2)*dx;
+[X1,X2] = ndgrid(x1,x2);
+% Now array subscripts are equivalent with Cartesian coordinates
+% x1 axis points downwards and x2 axis is in horizontal direction
+% x1 = X1(:,1) is a column vector in vertical direction
+% x2 = X2(1,:) is a row vector in horizontal direction
+R = sqrt(X1.^2 + X2.^2);
+cgrid = c_s * (R < a) + c_0 * (R >= a);
+rhogrid = rho_s * (R < a) + rho_0 * (R >= a);
+
+% DATA from Thorbecke's finite difference code
+[Ptsct, Pfinc, Pfsct, f_out]=FD_mod_grid( xS, xR, ntrcv, dtrcv, dxR, cgrid, rhogrid, f_in );
+
+f=f_out % nearest computed discrete frequency
+
+% Compare with analytical solution ---------------------------------------
+[P_inc,P_sct,xR] = ForwardCircle( xS, xR, gsize, dxR, c_0, c_s, rho_0, rho_s, f );
+
+set(figure(1),'Units','centimeters','Position',[1 1 30 10]);
+subplot(1,3,1)
+ plot(xR(2,:),real(Pfinc),'LineWidth',1.2);
+ title('Real part of P^{inc}'); axis tight; hold on;
+ plot(xR(2,:),real(P_inc),'--r','LineWidth',1.2);
+ axis tight; hold off;
+subplot(1,3,2)
+ plot(xR(2,:),imag(Pfinc),'LineWidth',1.2);
+ title('Imaginary part of P^{inc}'); axis tight; hold on;
+ plot(xR(2,:),imag(P_inc),'--r','LineWidth',1.2);
+ axis tight; hold off;
+subplot(1,3,3)
+ plot(xR(2,:), abs(Pfinc),'LineWidth',1.2);
+ title('Absolute value of P^{inc}'); axis tight; hold on;
+ plot(xR(2,:), abs(P_inc),'--r','LineWidth',1.2);
+ axis tight; hold off
+legendtitle1 = sprintf('f=%.2fHz FiniteDiff', f);
+legendtitle2 = sprintf('f=%.2fHz Analytic ', f);
+legend(legendtitle1,legendtitle2,'Location','Best');
+
+set(figure(2),'Units','centimeters','Position',[9 12 10 10]);
+error = abs(P_sct(:)-Pfsct(:))./abs(Pfsct(:));
+plot(error,'LineWidth',1.2); title('Relative error'); axis tight;
+
+set(figure(3),'Units','centimeters','Position',[1 1 30 10]);
+subplot(1,3,1)
+ plot(xR(2,:),real(Pfsct),'LineWidth',1.2);
+ title('Real part of P_{sct}'); axis tight; hold on;
+ plot(xR(2,:),real(P_sct),'LineWidth',1.2);
+subplot(1,3,2)
+ plot(xR(2,:),imag(Pfsct),'LineWidth',1.2);
+ title('Imaginary part of P^{sct}'); axis tight; hold on;
+ plot(xR(2,:),imag(P_sct),'LineWidth',1.2);
+subplot(1,3,3)
+ plot(xR(2,:), abs(Pfsct),'LineWidth',1.2);
+ title('Absolute value of P^{sct}'); axis tight; hold on;
+ plot(xR(2,:), abs(P_sct),'LineWidth',1.2);
+ axis tight; hold off
+legendtitle1 = sprintf('f=%.2fHz FiniteDiff', f);
+legendtitle2 = sprintf('f=%.2fHz Analytic ', f);
+legend(legendtitle1,legendtitle2,'Location','Best');