function s1 = stabrnd(a1,n)
u=rand(n,1);  %comment out this if you want u as an argument
% rewrite so that we can permit vectors for a.
u = (u*pi)-.5*pi;
e = exprnd(1,n,1);
s1=[];
for a=a1
s = sin(a*u).*(cos((1-a)*u)./e).^(1/a-1).*(cos(u)).^(-1/a);
s1=[s1 s];
end
