N = 15 ;
beta_a = 10;
alpha_n = 10 ;
eta = zeros(1,N) ;
eta(1)=15;
for i=1:N-1
theta_i = asind(sind(beta_a)*sind(eta(i)));
theta_n = atand(tand(beta_a)*cosd(eta(i)))- alpha_n;
phi_i = asind(sqrt(2)*sind(theta_i));
phi_n = acosd(tand(theta_i)/tand(phi_i))-theta_n ;
eta(i+1) = atand((tand(i)*cosd(phi_n-alpha_n)-cosd(alpha_n)*tand(phi_i))/sind(phi_n));
end
plot(1:N,eta)