%{
Function getTacEnergy returns the energy for Maier-Saupe energy
for the auxiliary variables eta, mu, nu.
This version of the function includes up to cubic gradient terms
%}

function E = getTacEnergyMS(x,u,a,kappa2,kappa3)
    size = [length(x),length(x)];
    dx = x(2) - x(1);
    fb = zeros(length(x));
    %Calculate derivative fields
	dux = zeros(3,length(x),length(x));
	duy = zeros(3,length(x),length(x));
	for i=2:length(x)-1
		for j=2:length(x)-1
			dux(1,i,j) = xderiv(u,1,i,j,dx);
			duy(1,i,j) = yderiv(u,1,i,j,dx);
			dux(2,i,j) = xderiv(u,2,i,j,dx);
			duy(2,i,j) = yderiv(u,2,i,j,dx);
			dux(3,i,j) = xderiv(u,3,i,j,dx);
			duy(3,i,j) = yderiv(u,3,i,j,dx);
		end
	end
    %Calculate free energy densities
    fb = -2*a*((1/3)*reshape(u(1,:,:),size).^2 + reshape(u(2,:,:),size).^2 + reshape(u(3,:,:),size).^2);
    fg = 2*((1/3)*(reshape(dux(1,:,:),size).^2 + reshape(duy(1,:,:),size).^2) + reshape(dux(2,:,:),size).^2 + reshape(duy(2,:,:),size).^2 + reshape(dux(3,:,:),size).^2 + reshape(duy(3,:,:),size).^2) + kappa2*(((2/3)*reshape(dux(1,:,:),size) + reshape(duy(3,:,:),size)).^2 + ((1/3)*reshape(duy(1,:,:),size) - (reshape(duy(2,:,:),size) + reshape(dux(3,:,:),size))).^2) + (2/9)*kappa3*(3*(reshape(u(1,:,:),size).*(reshape(duy(1,:,:),size).^2 + 3*(reshape(duy(2,:,:),size).^2 + reshape(duy(3,:,:),size).^2)) + 2*reshape(u(3,:,:),size).*(reshape(duy(1,:,:),size).*reshape(dux(1,:,:),size) + 3*reshape(duy(2,:,:),size).*reshape(dux(2,:,:),size) + 3*reshape(duy(3,:,:),size).*reshape(dux(3,:,:),size))) - reshape(u(1,:,:),size).*(reshape(duy(1,:,:),size).^2 + 3*reshape(duy(2,:,:),size).^2 + 3*reshape(duy(3,:,:),size).^2 - 2*(reshape(dux(1,:,:),size).^2 + 3*(reshape(dux(2,:,:),size).^2 + reshape(dux(3,:,:),size).^2))));
    f = fb + fg;
    E = trapz(x,trapz(x,f));
end