% Global parameter definitions
nsta=50;
aperture=100.0;
snr=10;
dt=0.025;
nt=1000;
rwsigma=1.0;
nrw=128;
lag0=[100,200,300,400];  %lag at x=0 in samples
u=[0.0,0.02,0.05,0.1]; % slowness values (must match lag0 size)
nbins=400;  % used for irregular plot display
nplanes=max(size(u));
% pseudostation stack parameters
sigma=10.0;
xmax=30.0;
x0out=0.0;
dxout=1.0;
nxout=100;

% Done with parameter definitions.
% First compute a ricker wavelet 
rw=ricker(nrw,dt,rwsigma);

x=rand(nsta,1).*aperture;
x=sort(x);
n=randn(nt,nsta);
% This is required for rand, but randn produces zero mean normal deviates
%n=n-0.5; 
n=n/snr;
d=zeros(nt,nsta);
for np=1:nplanes
for i=1:nsta
    lag=fix(lag0(np)+(u(np)*x(i))/dt);
    d(lag,i)=1.0;
end
end
% Use this to add noise filtered to bandwidth of Ricker wavelet
% For this tutorial we use Gaussian white noise
%d=d+n;
for i=1:nsta
    c=conv(d(:,i),rw);
    d(:,i)=c(1:nt);
end
d=d+n;
plotirregular(d,x,nbins);
title('Irregularly spaced input data');
stack=pseudostack(d,x,x0out,dxout,nxout,sigma,xmax);
figure;
imagesc(stack);
%We need this array of pseudostation positions to restore moveout
% after stack
xout=zeros(nxout,1);
for i=1:nxout
    xout(i)=x0out+(i-1)*dxout;
end
title('pseudostation stack with no moveout');
nu=max(size(u));
for i=1:nu
    dpwmo=pwmoveout(d,x,dt,-u(i));
    stack=pseudostack(dpwmo,x,x0out,dxout,nxout,sigma,xmax);
    stack=pwmoveout(stack,xout,dt,u(i));
    figure;
    imagesc(stack);
    tit=sprintf('pseudostation stack with u=%f',u(i));
    title(tit);
end