function wobbling_sh(trajectory)
% This code calculated wobble angle
trj = csvread(trajectory);

[~,~,ic] = unique(trj(:,5), 'rows');
trj(:,5) = ic;

T = cell(max(trj(:,5)),1);
W_rms = zeros(max(trj(:,5)),1);
W_abs = zeros(max(trj(:,5)),1);


for i = 1:1:max(trj(:,5))
   temp1 = []; temp2 = [];
   ind = find(trj(:,5) == i);
   T{i} = trj(ind,1:3);
   temp1 = T{i}(2:end,1:2) - T{i}(1:end-1,1:2);
   T{i}(1:end-1,1:2) = temp1;
   T{i}(end,:) = [];
   temp2 = atan(-T{i}(:,2)./T{i}(:,1));
   T{i}(:,4) = temp2;
   T{i}(:,5) = acos(cos(T{i}(:,3) - T{i}(:,4))); 
   T{i}(:,3:5) = 180/pi * T{i}(:,3:5);
   W_rms(i,1) = min(180-(mean((T{i}(:,5)).^2)).^0.5,(mean((T{i}(:,5)).^2)).^0.5); %rms wobbling angle
   W_abs(i,1) = min(180-(mean(abs(T{i}(:,5)))),mean(abs((T{i}(:,5))))); %abs wobbling angle
end

W_avg_rms = nanmean(W_rms(1:end,1))
W_std_rms = nanstd(W_rms(1:end,1))
W_avg_abs = nanmean(W_abs(1:end,1))
W_std_abs = nanstd(W_abs(1:end,1))

figure
histogram(W_rms,'normalization','pdf')

figure
histogram(W_abs,'normalization','pdf')

figure
plot(W_rms,'r')
hold on
plot(W_abs,'b')
hold off

end