function [Z, S, xi, Axes, B] = biaxial(a,b,c);
%% biaxial takes as inputs three rotational parameters:
%% The molecule  is rotated about its z-axis by a
%% random amount depending on c. If c=0, the
%% amount is 0, if c=1, all angles are equally likely.
%% 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.
%% If c=0, a and b do the same job as in nemord.
%%
%% Ouputs: Z S xi Axes B
%% S, Z and xi use the computed principal director,
%% Axes gives X and Y, B is the biaxiality measure.

%% Set up a system of 1000 molecules.
Qxx=[0 0 0;0 0 0;0 0 0];
Qyy=[0 0 0;0 0 0;0 0 0];
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);
   psi=c*2*pi*(rand-0.5);
   M1inf=theta*[0 0 (-sin(phi));0 0 cos(phi); sin(phi) (-cos(phi)) 0];
   M1 = expm(M1inf);
   M2=[cos(psi) -sin(psi) 0;sin(psi) cos(psi) 0;0 0 1];
   M = M1*M2;
   x=M(:,1); y=M(:,2); z=M(:,3);
   Qxx = Qxx+x*x';Qyy=Qyy+y*y';Qzz=Qzz+z*z';
end
Qzz=(3*Qzz/1000-eye(3))/2;
Qxx=(3*Qxx/1000-eye(3))/2;
Qyy=(3*Qyy/1000-eye(3))/2;

%Get the eigenvectors etc for each Q
[Ux Vx]=eig(Qxx); Vx=eig(Vx);
[Uy Vy]=eig(Qyy); Vy=eig(Vy);
[Uz Vz]=eig(Qzz); Vz=eig(Vz);

%Now pick out appropriate maximal eigenvectors
%to obtain the system director.
%Keep also the eigenvector associated with the biggest
%positive eigenvalue, and not in a set that's been used.
Qs=[Qxx Qyy Qzz];
Vectors = [Ux Uy Uz];
Values=[Vx; Vy; Vz];
[maxz, imaxz]= max(abs(Values));
Qzindex =floor((imaxz-1)/3)+1;
Z=Vectors(:,imaxz);
Vz=Values(3*Qzindex-2:3*Qzindex);
Values(3*Qzindex-2:3*Qzindex)=0;
[maxy, imaxy]= max(Values);
Y=Vectors(:,imaxy);
Y=(Y-Z'*Y)/norm(Y-Z'*Y);
X=cross(Y,Z);
Qyindex =floor((imaxy-1)/3)+1;
Qxindex = 6-(Qzindex+Qyindex);

%Shuffle the Q tensors to match the axes.
Qxx=Qs(:,3*Qxindex-2:3*Qxindex);
Qyy=Qs(:,3*Qyindex-2:3*Qyindex);
Qzz=Qs(:,3*Qzindex-2:3*Qzindex);

%Provide the various order parameters.
S = Z'*Qzz*Z;

%% Sort the eigenvalues of (new) Qzz and take lambda2-lambda3.
[L,Lindex]=sort(-abs(Vz));
L=Vz(Lindex);
xi=(L(2)-L(3))/2;

%% Find the B biaxiality parameter.
Axes=[X Y];
B =(X'*Qxx*X + Y'*Qyy*Y - X'*Qyy*X - Y'*Qxx*Y)/3;