function [beta] = ouo_theta_M(gamma, theta, M)
% angoli in radianti
%
beta_min=asin(1/M);
beta_max=pi/2;

beta_list=linspace(beta_min,beta_max,1000);

theta_list = ouo_beta_M(gamma, beta_list, M);
[~,closestIndex] = min(abs(theta_list-theta));
beta1a=beta_list(closestIndex-1);
beta1b=beta_list(closestIndex+1);

theta_list([closestIndex-1,closestIndex,closestIndex+1])=[];
[~,closestIndex] = min(abs(theta_list-theta));
beta2a=beta_list(closestIndex-1);
beta2b=beta_list(closestIndex+1);


beta_list=linspace(beta1a,beta1b,1000);
theta_list = ouo_beta_M(gamma, beta_list, M);
[~,closestIndex] = min(abs(theta_list-theta));
beta1=beta_list(closestIndex);

beta_list=linspace(beta2a,beta2b,1000);
theta_list = ouo_beta_M(gamma, beta_list, M);
[~,closestIndex] = min(abs(theta_list-theta));
beta2=beta_list(closestIndex);

beta = sort([beta1,beta2]);
end