function [d,s0]=psstutorial(nsta,aperture,snr,sigma,xmax,nxout)
%  Pseudostation stacking tutorial using sequence of 
% dipping layers.  Creates a suite of simulated dipping
% layers with the impulse response convolved with a ricker
% wavelet intended to approximate upper frequency limit
% possible with teleseismic data.  White noise is added
% with a scale controlled by snr parameter.  snr here is 
% defined as peak of impulse divided by range of random
% numbers added to data.  
%
% Tutorial produces two types of plots.  First it uses
% imagesc to attempt to display input, randomly spaced
% station data.  It then produces one plot for each of 
% the defined dipping layers using pseudostation stacking.
% For each pseudostation stack result the data are first
% moveout corrected to flatten each successive dipping 
% feature, stack, and then the moveout is restored.  
% Tutorial shows that pseudostation stacking interpolates
% data AND can be used to focus dipping features by varying
% the plane wave stacking velocity.  This is the essence
% of potential benefits of the plane wave method of 
% Poppeliers and Pavlis (2003).  
%
% Arguments:
%  nsta - number of actual stations simulated
%  aperture - 1D array aperture (km) simulated.  nsta simulated
%   stations will be randomly distributed from 0 to aperture
%  snr - signal to noise ratio of simulated data.  snr is
%   defined here as ratio of peak amplitude of dipping 
%   feature impulse and the range of random numbers added
%   to simulate noise.  Note noise is simple random numbers
%   from a uniform distribution (matlab rand function) 
%   offset to a zero mean process.
%  sigma - width (km) of Gaussian smoother used in pseudostation
%   stack computation.
%  xmax - outer cutoff parameter for pseudostation stacking.
%   Stations more than xmax distance from the pseudostation 
%   point will be completely dropped from the stack.
%  nxout - number of evenly spaced output pseudostations.
%   These points will be uniformly distributed between 0 and
%   aperture.
% Global parameter definitions
%nsta=50;
%aperture=100.0;
%snr=10;
dt=0.025;
nt=1000;
tmin=0.0;
tmax=(nt-1)*dt;
t=zeros(1000,1);
for i=1:nt
    t(i)=(i-1)*dt;
end
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;
dxout=aperture/nxout;

% 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;
offset=nrw/2;
offset=offset+1;
for i=1:nsta
    c=conv(d(:,i),rw);
    d(:,i)=c(offset:offset+nt-1);
end
d=d+n;
plotirregular(d,x,nbins);
title('Irregularly spaced input data - image format');
figure;
wigb(d,1.0,x,t);
s0=pseudostack(d,x,x0out,dxout,nxout,sigma,xmax);
%figure;
%imagesc(stack);
%title('pseudostation stack with no moveout');

%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

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);
    % Remove these to plot with wiggle overlay
    % hold;
    % wigb(stack);
    tit=sprintf('pseudostation stack with u=%f',u(i));
    title(tit);
end