function y = rsn(n,xi,omega,alpha)
    delta = alpha/sqrt(1+alpha^2);
    tn    = randn(2,n);
    chi   = abs(tn(1,:));
    nrv   = tn(2,:);
    z     = delta*chi+sqrt(1-delta^2)*nrv;
    y     = xi+omega*z;
end
