% Control3_1_JMS.m
% Solution to problem 1 of control 3. 
% Jose M Soler, May.2022

% Set parameters
clear
Phi = pi/6;        % prism angle (rad)
A = 1.728;         % refraction index parameters, n=A+B/Lambda^2
B = 1.342e-2;      % parameter B in microns^2
LambdaMin = 0.38;  % min wavelength (microns)
LambdaMax = 0.78;  % max wavelength (microns)
dLambda = 0.01;    % wavelegth interval (microns);
dAlpha = pi/128;   % incidence angle interval (rad)

% Find refraction angles
Lambda = (LambdaMin:dLambda:LambdaMax)';

% Plot refraction angle as a function of incidence angle
figure(11)
Alpha = (0:dAlpha:pi/2)';
Lambda = [LambdaMin,LambdaMax];
Beta = PrismRefraction_JMS(Phi,A,B,Lambda,Alpha);
plot(Alpha/pi,Beta/pi,'LineWidth',1.5)
legend(['\lambda=',num2str(Lambda(1)),'\mum'], ...
       ['\lambda=',num2str(Lambda(2)),'\mum']);
xlabel('\alpha/\pi')
ylabel('\beta/\pi')
grid on
set(gca,'FontSize',15)

% Plot refraction angle as a function of wavelength
figure(12)
Alpha = (0:pi/8:pi/2)';
nAlpha = numel(Alpha);
Lambda = (LambdaMin:dLambda:LambdaMax)';
Beta = PrismRefraction_JMS(Phi,A,B,Lambda,Alpha);
plot(Lambda,Beta/pi,'LineWidth',1.5)
legend([repmat('\alpha=',nAlpha,1),num2str(Alpha/pi),repmat('\pi',nAlpha,1)], ...
       'Location','NorthWest');
xlabel('\lambda (\mum)')
ylabel('\beta/\pi')
grid on
set(gca,'FontSize',15)

% Find incidence angles such that Alpha=Beta
for iLambda = 1:numel(Lambda)
  myfunc = @(Alpha)(PrismRefraction_JMS(Phi,A,B,Lambda(iLambda),Alpha)-Alpha);
  alpha(iLambda) = bisectionRoot_JMS( myfunc, 0, pi/2, 1e-6 );
end

% Plot incidence angles such that Alpha=Beta
figure(13)
plot(Lambda,alpha/pi,'LineWidth',1.5)
xlabel('\lambda (\mum)')
ylabel('\alpha/\pi')
grid on
set(gca,'FontSize',15)
