function [thetas, values] = h3d_compensate_one_rotation(p, ht1, ht2, axis, current_guess,  max_correction, spice)

	support = h3d_get_one_perpendicular(axis);

	num_alpha_steps = p.cyl_num_alpha_steps;
	num_beta_steps =  p.cyl_num_beta_steps;

	amplitude = deg2rad(p.cyl_amplitude_deg);
	weight_by_alpha = 1;
	cyl1 = h3d_project_hs_on_cylinder(ht1, axis, support, num_alpha_steps,num_beta_steps, amplitude, weight_by_alpha, spice);
	cyl2 = h3d_project_hs_on_cylinder(ht2, current_guess*axis, current_guess*support, num_alpha_steps,num_beta_steps,amplitude,weight_by_alpha,spice);
	
	% normalization
	cyl1 = cyl1 / matrix_norm(cyl1);
	cyl2 = cyl2 / matrix_norm(cyl2);

	% cerchiamo correlazione su beta
	
	max_correction_steps = ceil(max_correction * (num_beta_steps/(2*pi)));
	lags = -max_correction_steps:max_correction_steps;

	corr = zeros(1,numel(lags));
	for alpha=1:num_alpha_steps
		% If many, use the FFT algorithm
		if numel(lags)>50
			row_corr = circular_correlation2(cyl2(alpha,:),cyl1(alpha,:),lags);
		else
			row_corr = circular_correlation(cyl2(alpha,:),cyl1(alpha,:),lags);
		end
		corr = corr + row_corr;
	end

	% Find local maxima of correlation
	[coords, values] = h3d_matrix_local_maxima(corr);
	indices = coords(2,:); % take only column
	
	% Sort
	[unused,I] = sort(values, 'descend');
	values = values(I);
	indices = indices(I);
	

	% Convert indices to degrees
	thetas = lags(indices) * 2 * pi / num_beta_steps;


	% if near bounds, discard
	for i=1:numel(indices)
		if( (indices(i) == 1) || (indices(i) == numel(corr)) )
			fprintf('detected max bounds\n');
			thetas(i) = 1000;
		end
	end
	
	valid = not(thetas > 999);
	
	values = values(valid);
	indices = indices(valid);
	thetas = thetas(valid);
	
	
	if p.debug_compensation_show
		lags_degrees = lags*(360/num_beta_steps);
		n = 3;
		subplot(n,1,1);
		imagesc(cyl1);
		xlabel('\beta')
		ylabel('\alpha')
		subplot(n,1,2);
		imagesc(cyl2);
		xlabel('\beta')
		ylabel('\alpha')

		subplot(n,1,3); 
		plot(lags_degrees,corr,'x-'); hold on;
		title('xcorrelation');
		xlabel('\theta (degrees)');
		
		for theta=rad2deg(thetas)
			plot([theta theta],[0, max(values)],'k-');
		end
		hold off;
		
		pause(0.01)
%		pause(0.5)
%		pause
	end
	
function res= matrix_norm(m)
	n = 0;
	for i=1:size(m,1)
	for j=1:size(m,2)
		n = n + m(i,j)^2;
	end
	end
	
	res =sqrt(n);

function c = circular_correlation(a,b,lags)
	aa = [a a];
	n = numel(a);
	for k=1:numel(lags)
		lag = mod(lags(k),n);
		aa_part = aa(1+lag:1+lag+(n-1));
		c(k) =  aa_part * b';
	end
	

function res=circular_correlation2(a,b,lags)
	cc=ifft(fft(a).*fft(b(length(b):-1:1)));

	for i=1:numel(lags)
		l = mod(lags(i)-1, numel(a));
		res(i) = cc(l+1);
	end