if false
	axis_a = m1;
	support_a = normalize(cross(m1,m2));
	[thetas_a, values_a] = h3d_compensate_one_rotation(ht1u, ht2u, axis_a, support_a, r1, 1);

	r2 = rot_axis_angle(r1*axis_a, thetas_a(1));

	fprintf('After compensating %f deg, max error is %f deg.\n', rad2deg(thetas_a(1)), h3d_rotation_error(R, r2*r1));
	fprintf(' axis_a %1.2f %1.2f %1.2f \n', axis_a(1),axis_a(2),axis_a(3));
	fprintf(' quality = %f \n', values_a(1));

	axis_b = support_a;
	support_b = normalize(cross(axis_b, [0;0;1]));
	[thetas_b, values_b] = h3d_compensate_one_rotation(ht1u, ht2u, axis_b, support_b, r2*r1, 1);

	r3 = rot_axis_angle(r2*r1*axis_b, thetas_b(1));

	fprintf('After compensating %f deg, max error is %f deg.\n', rad2deg(thetas_b(1)), h3d_rotation_error(R, r3*r2*r1));
	fprintf(' axis_b %1.2f %1.2f %1.2f \n', axis_b(1),axis_b(2),axis_b(3));
	fprintf(' quality = %f \n', values_b(1));

	axis_c = normalize(cross(axis_a, axis_b));
	support_c = normalize(cross(axis_c, [0;0;1]));
	
	[thetas_c, values_c] = h3d_compensate_one_rotation(ht1u, ht2u, axis_c, support_c, r3*r2*r1, 1);

	r4 = rot_axis_angle(r3*r2*r1*axis_c, thetas_c(1));

	fprintf('After compensating %f deg, max error is %f deg.\n', rad2deg(thetas_c(1)), h3d_rotation_error(R, r4*r3*r2*r1));
	fprintf(' axis_c %1.2f %1.2f %1.2f \n', axis_c(1),axis_c(2),axis_c(3));
	fprintf(' quality = %f \n', values_c(1));

	axis_d = support_c;
	support_d = normalize(cross(support_c, r4*r3*r2*r1*support_c));
	[thetas_d, values_d] = h3d_compensate_one_rotation(ht1u, ht2u, axis_d, support_d, r4*r3*r2*r1, 1);

	r5 = rot_axis_angle(r4*r3*r2*r1*axis_d, thetas_d(1));

	fprintf('After compensating %f deg, max error is %f deg.\n', rad2deg(thetas_d(1)), h3d_rotation_error(R, r5*r4*r3*r2*r1));
	fprintf(' axis_d %1.2f %1.2f %1.2f \n', axis_d(1),axis_d(2),axis_d(3));
	fprintf(' quality = %f \n', values_d(1));

	axis_e = support_d;
	support_e = normalize(cross(support_d, Rest*r4*r3*r2*r1*support_d));
	[thetas_e, values_e] = h3d_compensate_one_rotation(ht1u, ht2u, axis_e, support_e, r5*r4*r3*r2*r1, 1);

	r6 = rot_axis_angle(r5*r4*r3*r2*r1*axis_d, thetas_e(1));

	fprintf('After compensating %f deg, max error is %f deg.\n', rad2deg(thetas_e(1)), h3d_rotation_error(R, r6*r5*r4*r3*r2*r1));
	fprintf(' axis_e %1.2f %1.2f %1.2f \n', axis_e(1),axis_e(2),axis_e(3));
	fprintf(' quality = %f \n', values_e(1));


	Rest = r6 * r5*r4*r3*r2*r1;
end