function V = gradient_fun(V, operator, options)
	% calculate velocity gradients and related quantities from
	% the velocity field provided
	% 
	% the particular operator used to calculate the gradients is given by
	% operator.
	%
	% options is a cell array which determines which gradient quantities
	% need to be calculated, in addition to the velocity gradient tensor
	% 
	% these are:
	% - deviatoric		remove divergence from VGT isotropically for all
	%						calculations
	% - QA				2nd invariant of VGT
	% - RA				3rd invariant of VGT
	% - QS				2nd invariant of RST
	% - RS				3rd invariant of RST
	% - QW				2nd invariant of RRT
	% - sigma			eigenvalues of RST
	% - lab_to_str		eigenframe of RST
	%						defined with the convention that in the strain frame
	%						w1*w2 >= 0 and w1*w3 >= 0
	%
	% the result is returned in V, in the following fields
	% Aij					velocity gradient tensor
	% Sij					rate of strain tensor
	% Wij					rate of rotation tensor
	% omega				vorticity vector
	% 
	% the following are optional:
	% Lij					SikSkj + WikWkj = S^2 + W^2
	% QA					2nd invariant of VGT
	% RA					3rd invariant of VGT
	% QS					2nd invariant of RST
	% RS					3rd invariant of RST
	% QW					2nd invariant of RRT
	% sigma				eigenvalues of RST
	% lambda				eigenvalues of lambda tensor
	% lab_to_str		eigenvectors of RST
	% lab_to_lda		eigenvectors of lambda tensor
	% lab_to_ssf		eigenvectors of S^2, same as eigenvectors of S,
	%						ut ordered by magnitude instead
	% lab_to_wwf		eigenvectors of W^2, frame composed of vorticity and
	%						two randomly oriented orthogonal vectors
	
	%% warn if user has supplied a misspelled or unrecognised keyword
	if nargin < 3
		options			= {};
	end
	recognised_words	= {'deviatoric', 'Lij', 'QA', 'RA', 'QS', 'RS', 'QW', 'RW', 'VS', 'sigma', 'lambda', 'lab_to_str', 'lab_to_lda', 'lab_to_ssf', 'lab_to_wwf'};
	n_options			= length(options);
	for i = 1 : n_options
		keyword			= options{i};
		if ~any(strcmpi(keyword, recognised_words))
			warning('unrecognised option ''%s''\n', keyword);
		end
	end
	
	%% velocity gradient tensor
	dx				= V.winCtrsX(2) - V.winCtrsX(1);
	dy				= V.winCtrsY(2) - V.winCtrsY(1);
	dz				= V.winCtrsZ(2) - V.winCtrsZ(1);
	n_vectors	= size(V.ux);

	% velocity gradients
	V.Aij			= zeros([3 3 n_vectors]);
	[V.Aij(1,1,:,:,:), V.Aij(1,2,:,:,:), V.Aij(1,3,:,:,:)] = operator(V.ux, dx, dy, dz);
	[V.Aij(2,1,:,:,:), V.Aij(2,2,:,:,:), V.Aij(2,3,:,:,:)] = operator(V.uy, dx, dy, dz);
	[V.Aij(3,1,:,:,:), V.Aij(3,2,:,:,:), V.Aij(3,3,:,:,:)] = operator(V.uz, dx, dy, dz);
	
	% deviatoric correction to Aij and Sij
	if any(strcmpi('deviatoric', options))
		div					= V.Aij(1,1,:,:,:) + V.Aij(2,2,:,:,:) + V.Aij(3,3,:,:,:);
		V.Aij(1,1,:,:,:)	= V.Aij(1,1,:,:,:) - div/3;
		V.Aij(2,2,:,:,:)	= V.Aij(2,2,:,:,:) - div/3;
		V.Aij(3,3,:,:,:)	= V.Aij(3,3,:,:,:) - div/3;
	end
	
	% strain and rotation tensors
	V.Sij			= 1/2*(V.Aij + permute(V.Aij, [2 1 3 4 5]));
	V.Wij			= 1/2*(V.Aij - permute(V.Aij, [2 1 3 4 5]));
	V.omega		= cat(4, squeeze(V.Wij(3,2,:,:,:)*2), ...
								squeeze(V.Wij(1,3,:,:,:)*2), ...
								squeeze(V.Wij(2,1,:,:,:)*2));
	V.omega		= permute(V.omega, [4 1 2 3]);
	
	%% lambda tensor
	if any(strcmpi('Lij', options)) || any(strcmpi('lab_to_lda', options)) || any(strcmpi('lambda', options))
		V.Lij		= mtimesx(V.Sij,V.Sij) + mtimesx(V.Wij,V.Wij);
	end
	
	%% QA
	if any(strcmpi('QA', options))
		V.QA		=  -1/2 * (V.Aij(1,1,:).^2 ...
								+ V.Aij(2,2,:).^2 ...
								+ V.Aij(3,3,:).^2 ...
								+ 2*V.Aij(1,2,:).*V.Aij(2,1,:) ...
								+ 2*V.Aij(1,3,:).*V.Aij(3,1,:) ...
								+ 2*V.Aij(2,3,:).*V.Aij(3,2,:));
		V.QA		= reshape(V.QA, n_vectors);
	end
	%% RA
	if any(strcmpi('RA', options))
		V.RA		=  -(	V.Aij(1,1,:).*V.Aij(2,2,:).*V.Aij(3,3,:) ...
						-  V.Aij(1,1,:).*V.Aij(2,3,:).*V.Aij(3,2,:) ...
						-  V.Aij(1,2,:).*V.Aij(2,1,:).*V.Aij(3,3,:) ...
						+  V.Aij(1,2,:).*V.Aij(2,3,:).*V.Aij(3,1,:) ...
						+  V.Aij(1,3,:).*V.Aij(2,1,:).*V.Aij(3,2,:) ...
						-  V.Aij(1,3,:).*V.Aij(2,2,:).*V.Aij(3,1,:));
		V.RA		= reshape(V.RA, n_vectors);
	end 
	%% QS
	if any(strcmpi('QS', options))
		V.QS		= -1/2 * ( V.Sij(1,1,:).^2 ...
								+ V.Sij(2,2,:).^2 ...
								+ V.Sij(3,3,:).^2 ...
								+ 2*V.Sij(1,2,:).^2 ...
								+ 2*V.Sij(1,3,:).^2 ...
								+ 2*V.Sij(2,3,:).^2);
		V.QS		= reshape(V.QS, n_vectors);
	end
	%% RS
	if any(strcmpi('RS', options))
		V.RS		=  -(	V.Sij(1,1,:).*V.Sij(2,2,:).*V.Sij(3,3,:) ...
						-  V.Sij(1,1,:).*V.Sij(2,3,:).*V.Sij(3,2,:) ...
						-  V.Sij(1,2,:).*V.Sij(2,1,:).*V.Sij(3,3,:) ...
						+  V.Sij(1,2,:).*V.Sij(2,3,:).*V.Sij(3,1,:) ...
						+  V.Sij(1,3,:).*V.Sij(2,1,:).*V.Sij(3,2,:) ...
						-  V.Sij(1,3,:).*V.Sij(2,2,:).*V.Sij(3,1,:));
		V.RS		= reshape(V.RS, n_vectors);
	end
	%% QW
	if any(strcmpi('QW', options))
		V.QW		= sum(V.omega.^2, 1) / 4;
		V.QW		= squeeze(V.QW);
	end
	%% RW
	% enstrophy production term, wi Sij wj
	if any(strcmpi('RW', options))
		omega		= reshape(V.omega, [3 1 n_vectors]);
		V.RW		= mtimesx(shiftdim(V.omega,-1), mtimesx(V.Sij, omega));
		V.RW		= squeeze(V.RW);
	end
	%% VS
	% norm of vortex stretching term
	% Sij wj Sik wk
	if any(strcmpi('VS', options))
		vs			= mtimesx(V.Sij, reshape(V.omega, [3 1 n_vectors]));
		V.VS		= squeeze(sum(vs.^2, 1));
	end
	
	%% strain eigenvalues and eigenvectors
	if any(strcmpi('lab_to_str', options))
		% strain eigenvalues and eigenvectors
		% systematically choose RH basis
		% need to compute determinant of str_to_dns, little bit cumbersome but it
		% works
		[V.sigma, str_to_lab]	= ndsyev(single(V.Sij));
		V.sigma						= reshape(V.sigma, [3 n_vectors]);
		lab_to_str					= permute(reshape(str_to_lab, [3 3 prod(n_vectors)]), [2 1 3]);

		det_lab_to_str				= shiftdim(nddet(lab_to_str), -2);
		lab_to_str					= bsxfun(@times, lab_to_str, det_lab_to_str);

		% orient strain basis so w1w2>0 and w1w3 > 0
		lab_omega					= reshape(V.omega, [3 1 prod(n_vectors)]);
		str_omega					= mtimesx(lab_to_str, lab_omega);
		str_omega					= squeeze(str_omega);
		perm_vec						= [sign(str_omega(2,:) .* str_omega(3,:));
											sign(str_omega(1,:) .* str_omega(3,:));
											sign(str_omega(1,:) .* str_omega(2,:))];
		perm_vec						= reshape(perm_vec, [3 1 prod(n_vectors)]);
		lab_to_str					= bsxfun(@times, lab_to_str, perm_vec);
		str_omega					= mtimesx(lab_to_str, lab_omega);
		V.lab_to_str				= reshape(lab_to_str, [3 3 n_vectors]);
		V.str_omega					= reshape(str_omega, [3 n_vectors]);
	elseif any(strcmpi('sigma', options))
		%% just strain eigenvalues
		V.sigma						= ndsyev(single(V.Sij));
	end
	
	%% eigenbasis of S^2
	% pretty much the same as eigenbasis of S, but with ordering by
	% descending magnitude, not descending value
	if any(strcmpi('lab_to_ssf', options))
		% strain eigenvalues and eigenvectors
		% systematically choose RH basis
		% need to compute determinant of str_to_dns, little bit cumbersome but it
		% works
		[~, ssf_to_lab]			= ndsyev(single(V.Sij), 'u-');
		lab_to_ssf					= permute(reshape(ssf_to_lab, [3 3 prod(n_vectors)]), [2 1 3]);

		det_lab_to_ssf				= shiftdim(nddet(lab_to_ssf), -2);
		lab_to_ssf					= bsxfun(@times, lab_to_ssf, det_lab_to_ssf);

		% orient S^2 basis so w1w2>0 and w1w3 > 0
		lab_omega					= reshape(V.omega, [3 1 prod(n_vectors)]);
		ssf_omega					= mtimesx(lab_to_ssf, lab_omega);
		ssf_omega					= squeeze(ssf_omega);
		perm_vec						= [sign(ssf_omega(2,:) .* ssf_omega(3,:));
											sign(ssf_omega(1,:) .* ssf_omega(3,:));
											sign(ssf_omega(1,:) .* ssf_omega(2,:))];
		perm_vec						= reshape(perm_vec, [3 1 prod(n_vectors)]);
		lab_to_ssf					= bsxfun(@times, lab_to_ssf, perm_vec);
		V.lab_to_ssf				= reshape(lab_to_ssf, [3 3 n_vectors]);
		
		ssf_omega					= mtimesx(lab_to_ssf, lab_omega);
		V.ssf_omega					= reshape(ssf_omega, [3 n_vectors]);
	end
	
	%% lambda eigenvalues and eigenvectors
	if any(strcmpi('lab_to_lda', options))
		% lambda eigenvalues and eigenvectors
		% these are eigenvalues/vectors of the Lambda tensor
		% L = S^2 + W^2
		% its intermediate eigenvalue lambda2 defines the lambda2 criterion
		% of vortex identification
		
		% systematically choose RH basis
		% need to compute determinant of str_to_dns, little bit cumbersome but it
		% works
		[V.lambda, lda_to_lab]	= ndsyev(single(V.Lij));
		V.lambda						= reshape(V.lambda, [3 n_vectors]);
		lab_to_lda					= permute(reshape(lda_to_lab, [3 3 prod(n_vectors)]), [2 1 3]);

		det_lab_to_lda				= shiftdim(nddet(lab_to_lda), -2);
		lab_to_lda					= bsxfun(@times, lab_to_lda, det_lab_to_lda);

		% orient strain basis so w1w2>0 and w1w3 > 0
		lab_omega					= reshape(V.omega, [3 1 prod(n_vectors)]);
		lda_omega					= mtimesx(lab_to_lda, lab_omega);
		lda_omega					= squeeze(lda_omega);
		perm_vec						= [sign(lda_omega(2,:) .* lda_omega(3,:));
											sign(lda_omega(1,:) .* lda_omega(3,:));
											sign(lda_omega(1,:) .* lda_omega(2,:))];
		perm_vec						= reshape(perm_vec, [3 1 prod(n_vectors)]);
		lab_to_lda					= bsxfun(@times, lab_to_lda, perm_vec);
		lda_omega					= mtimesx(lab_to_lda, lab_omega);
		V.lab_to_lda				= reshape(lab_to_lda, [3 3 n_vectors]);
		V.lda_omega					= reshape(lda_omega, [3 n_vectors]);
	elseif any(strcmpi('lambda', options))
		%% just lambda eigenvalues
		V.lambda						= ndsyev(single(V.Lij));
	end
	
	%% W^2 eigenvaluse and eigenvectors
	if any(strcmpi('lab_to_wwf', options))
		% W^2 eigenvalues and eigenvectors
		[~, wwf_to_lab]			= ndsyev(single(mtimesx(V.Wij, V.Wij)));
		V.lab_to_wwf				= permute(wwf_to_lab, [2 1 3 4 5]);
	end
end