function [Beta,Gamma_a,Gamma_b] = PrismRefraction_JMS(Phi,A,B,Lambda,Alpha)
% Find refraction angle by a prism
% Input:
%   Phi       : prism angle (rad)
%   A,B       : refraction index parameters: n=A+B/lambda^2, B in microns^2
%   Lambda(N) : wavelength (microns)
%   Alpha(L)  : incidence angle, relative to normal to the incidence surface (rad)
% Output:
%   Beta(N,L)    : angle of the refracted ray, relative to normal to the output surface (rad)
%   Gamma_a(N,L) : angle of the internal ray, relative to normal to incidence surface (rad)
%   Gamma_b(N,L) : idem, relative to the normal to the output surface

% Reshape input arrays to expected shapes
N = numel(Lambda);
L = numel(Alpha);
Lambda = reshape(Lambda,N,1);
Alpha = reshape(Alpha,1,L);

% Find refraction angles
n = A+B./Lambda.^2;
Gamma_a = asin(sin(Alpha)./n);
Gamma_b = Phi-Gamma_a;
Beta = asin(n.*sin(Gamma_b));

end
