% Réseau de diffraction en réflection "blazé"
clear; close all
rad = 180/pi;

% une longueur d'onde
lambda = 550e-9
N = 4000
a = 0.01 / N

theta_i = 10/rad;
b  = a*cos(theta_i)  % 0.5e-6

theta = [-1 : 0.001 : 1];

% maximum de réflection dans la direction normale au réseau
% beta = pi*b/lambda * sin(theta);
% alfa = pi*a/lambda * sin(theta);
beta = pi*b/lambda * sin(theta + theta_i);
alfa = pi*a/lambda * (sin(theta) + sin(theta_i));

I = 1/ N^2 * (sinc1(beta)).^2 .* (sin(N*alfa) ./ sin(alfa)).^2;
figure
plot (theta*rad, I, 'r');
xlabel ('\theta (deg)'); grid on
tit = ['\theta_i = 10 deg - \lambda = ' num2str(lambda*1e9) ' nm - a = ' num2str(a*1e6) ' \mum - b = ' num2str(b*1e6) ' \mum'];
title (tit)

% maximum de réflection dans la direction correspondante à l'ordre 2 pour
% lambda (550 nm)

theta_2 = asin (2 * lambda / a - sin(theta_i))
theta_2_deg = theta_2 * rad
gamma = (theta_2 + theta_i) /2
gamma_deg = gamma*rad
beta = pi*b/lambda * sin(theta + theta_i - 2*gamma);
I = 1/ N^2 * (sinc1(beta)).^2 .* (sin(N*alfa) ./ sin(alfa)).^2;
figure
plot (theta*rad, I, 'r');
xlabel ('\theta (deg)'); grid on
tit = ['\theta_i = 10 deg - \lambda = ' num2str(lambda*1e9) ' nm - a = ' num2str(a*1e6) ' \mum - b = ' num2str(b*1e6) ' \mum'];
title (tit)