function stack=pseudostack(d,x,x0out,dxout,nxout,sigma,xmax)
% arguments:
%  d - data matrix assumed to be seismograms stored in columns
%  x - parallel array to data containing offset values for data column
%    members.  that is x(i) is taken to be the offset for seismic\
%    data in the column data(:,i). 
%  x0out - first pseudostation position
%  dxout - pseudostation spacing for output
%  nxout - number of output pseudostations
%  sigma - Gaussian width parameter (units of x) to define pseudostation weights
%  xmax - cutoff for x.  Outside this distance, data are dropped completely.
%
%  x and data must have matching sizes.  Data do not need to be sorted
% in x order.
dsize=size(d);
xsize=size(x);
if(xsize(1)~=dsize(2))
    'Columns in data array must equal rows in x'
    xsize(1)
    dsize(2)
    return
end
stack=zeros(dsize(1),nxout);
sigma2=2.0*sigma*sigma;
for ips=1:nxout
    s=zeros(dsize(1),1);
    sumwt=0.0;
    xp=x0out+(ips-1)*dxout;
    for i=1:dsize(2)
        delta=x(i)-xp;
        if(abs(delta)<xmax)
            wt=exp(-(delta*delta)/sigma2);
            s=s+wt.*d(:,i);
            sumwt=sumwt+wt;
        end
    end
    if(sumwt>0)
        s=s./sumwt;
        stack(:,ips)=s;
    end
end