function [Z, S, xi, Plane] = nemord(a,b);

%% nemord takes as inputs two rotational parameters:
%% An axis is chosen randomly in the xy plane,
%% the angle this makes with the x-axis depending
%% on b; if b=0, the axis is always the x-axis,
%% when b=1, all angles are equally probably.
%% The molecule's director is then rotated about
%% this axisby a random amount depending on a;
%% if a=0, the angle is zero, if a=1 all angles
%% are equally probable.
%%
%% Outputs: Z S xi Plane
%% Z and S are the main director and order parameter respectively.
%% Plane is the two other eigenvectors of Qzz, xi is the
%% biaxiality measure.

%% Set up a system of 1000 molecules.
Qzz=[0 0 0;0 0 0;0 0 0];
for i=1:1000
   theta=a*2*pi*(rand-0.5);
   phi=b*2*pi*(rand-0.5);
   Minf=theta*[0 0 (-sin(phi));0 0 cos(phi); sin(phi) (-cos(phi)) 0];
   Mrot=expm(Minf);
   z = Mrot(:,3);
   Qzz=Qzz+z*z';
end
Qzz=(3*Qzz/1000-eye(3))/2;

%Get the eigenvectors etc for Qzz  and sort them.
[Vec,V]=eig(Qzz); V=[V(1,1);V(2,2);V(3,3)];
[L,Lindex]=sort(-abs(V));
L=V(Lindex);
Eig=Vec(:,Lindex');

% Compute the order parameters.
Z = Eig(:,1);
Plane=Eig(:,2:3);
S=L(1);
xi=(L(2)-L(3))/2;