
% clear all
clf
clear

% 1D wave equation
% to investigate accuracy of FD operators in space and time


% parameters

nx=512;    % space dimension
isx=256;     % source location
nt=1000;   % time steps
c=3000.;   % ac. velocity
dx=20.;   % grid spacing
dt=.0015;   % time increment
a=1./5;    % dominant source period
it0=50;    % source time
ir=400;     % receiver location
ncheck=50;  % checkpoint
idisp=50;   % display frequency
ag=0;      % gaussian in space for source
nop=3;      % length of operator in space (one side)
tord=1;     % 2^tord is order of time extrapolation

% Stability

%disp(sprintf(' Stability criterium %g (pref. < 2.09439) ',c*dt/dx));
disp(sprintf(' Stability criterium %g (pref. < 0.2) ',c*dt/dx));


% spectral or local
disp(' Give numerical scheme : ');

disp(' 1) 3-point ');
disp(' 2) 5-point ');
disp(' 3) Fourier ');
im=input(' Which ? '); 

if im==1, nop=3; end
if im==2, nop=5; end

% Frequenz


disp(' Give source type : ')
disp(' 1) Ricker ')
disp(' 2) Delta ')
disp(' 3) Heaviside ')
ans=input(' Which ? ');

% initialization of the source time function

s=zeros(1,nt);     % source
seis=(1:nt)*0;     % seismogram


if ans == 1,

a=input(' Give frequency (Hz) ');a=1/a;
% ricker wavelet period a
t0=4*a/pi;
for it=1:nt,
t=it*dt;
s(it)= sqrt(pi)/2* (1/a*pi*(t-t0))*exp(-(1/a*pi*(t-t0))^2);
end
s=diff(s);
s=s/max(max(s));s(nt+1)=0;

elseif ans ==2,

s(10)=1;

elseif ans ==3,

s(10:nt)=1;

end



subplot(221),plot(s(1:nt/2),'b-')
title('Source time function')

% initialization of space dependent fields

p=(1:nx)*0.;        % pressure
lp=(1:nx)*0.;       % pressure 1s derivative
lp1=(1:nx)*0.;      % pressure 2nd derivative
pold=(1:nx)*0.;     % p at previous time step
pnew=(1:nx)*0.;     % p at next time step


gauss=(1:nx)*0;
if ag ~= 0,
gauss=exp(-ag*(gauss-isx).*(gauss-isx));
end
if ag == 0, gauss(isx)=1; end

subplot(222),plot(gauss(isx-10:isx+10),'b-'),title(' Gauss in space'),pause

disp(' Begin time extrapolation ... ')

t0=clock;

pause off
x=(1:nx)*dx/1000;

for it=1:nt,

% loop over time order


	% laplacian 4th order


	if im == 3;	
	lp1=sder1d(p,dx); 
	lp=sder1d(lp1,dx);
	else
	lp=fder1d(p,dx,nop);
	end

	% add source


	
	lp=lp+gauss*s(it);
	pnew=2*p-pold+c*c*dt*dt*lp;
	


pold=p;p=pnew;

% output
hseis(it)=p(ir);
if(rem(it,idisp)==0), subplot(212),
plot(x,p), title(sprintf(' Time step : %g ', it)),xlabel(' Distance (km) '),hold on
plot(x(ir),p(ir),'+')
hold off
drawnow, end


% timing
if(rem(it,ncheck)==0),
t1=clock;
elapsed=etime(t1,t0);
t0=t1;elapsed=elapsed/ncheck*nt;
disp(sprintf(' Predicted total : %g seconds',elapsed))
end


end


% seismogram compared with analytical solution 
src=int(s(1,1:nt-1)); 
green=(1:nt)*.0;
isig=((ir-isx)*dx/c)/dt; 
green(isig)=1.;
aseis=mconv(green,src);

aseis=aseis/max(abs(aseis));
hseis=hseis/max(abs(hseis));

plot((1:nt)*dt,aseis,'r-',(1:nt)*dt,hseis,'b-',(1:nt)*dt,hseis-aseis,'g-')
xlabel(' Time (sec) ') 
axis([.8 1.5, -Inf Inf])

% Calculation of error of solution
err=sqrt(sum((aseis(nt/2:nt)-hseis(nt/2:nt)).*(aseis(nt/2:nt)-hseis(nt/2:nt)))/sum(aseis(nt/2:nt).*aseis(nt/2:nt)))*100;
if im<3, title(sprintf(' %g point - %g order; T = %g s, Error = %g%% ',nop,tord*2,elapsed,err)), end
if im==3, title(sprintf(' Fourier - %g order; T = %g s, Error = %g%% ',tord*2,elapsed,err)), end


% End