function t=GenTestVector(t1,thetamax,N)

if nargin<3
    N=1;
end

%generate test vector within angle theta_max (stnd spherical coordinates)
t_theta=acos(rand(N,1).*(cos(thetamax)-1)+1);
t_phi=2*pi*rand(N,1);
tst=[sin(t_theta).*cos(t_phi) sin(t_theta).*sin(t_phi) cos(t_theta)]';

%decompose reference vector to its orthogonal directions
xp=cross(t1,[0 1 0]'); % these define orthogonal coordinates
yp=cross(xp,t1);
zp=t1;
%Make these unit vectors
xp=xp/sqrt(sum(xp.^2));
yp=yp/sqrt(sum(yp.^2));
zp=zp/sqrt(sum(zp.^2));

%repmat for vectors
xp=repmat(xp,[1 N]);
yp=repmat(yp,[1 N]);
zp=repmat(zp,[1 N]);

%rotate the vector to the new axis (repmat for dealing with N vectors)
t=repmat(tst(1,:),[3 1]).*xp+...
    repmat(tst(2,:),[3 1]).*yp+...
    repmat(tst(3,:),[3 1]).*zp;

