
% clear all
clf
clear

% 1D wave equation with Chebyshev
% requires get_cheby_matrix

% parameters

nx=125;    % space dimension
nt=5000;   % time steps
c=20.;   % ac. velocity
dt=.00001;   % time increment
it0=50;    % source time
ncheck=50;  % checkpoint
idisp=25;   % display frequency
ag=.05;      % gaussian in space for source


% space coordinates
for ix=0:nx, x(ix+1)=cos(pi*ix/nx);  end

% Stability
disp(sprintf(' Stability criterium %g (pref. < 0.2) ',c*dt/abs(max(diff(x)))));

% Frequenz
s=zeros(1,nt);     % source
seis=(1:nt)*0;     % seismogram
a=input(' Give source frequency (e.g. 1000) (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=s/max(max(s));s(nt+1)=0;

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

% initialization of space dependent fields

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


gauss=(1:nx)*0;
x0=-.25;
if ag ~= 0,
gauss=exp(-1/ag^2*(x-x0).*(x-x0));gauss=gauss';
end


subplot(222),plot(x,gauss,'b-'),title(' Gauss in space'),pause
pause

disp(' Begin time extrapolation ... ')

t0=clock;

pause off

D=get_cheby_matrix(nx);

for it=1:nt,

% loop over time order

    % 2nd derivative (D twice operated)
    lp=D*D*p';
	% add source
	lp=lp+gauss*s(it);
	pnew=2*p-pold+c*c*dt*dt*lp';
	pold=p;p=pnew;

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

end

