%% Oving 11 i bolgefysikk, program for aa regne ut
%% og plotte interferensmonstret for rektangulaer
%% og sirkulaer aapning.
%%
%% Bolgelengden er i alle tilfelle 633 nm (Bruker SI-enheter)
lambda = 633e-9;
%% Rektangulaer aapning.
%%
%% Vi har disse lengdeparametrene: 
%% Sidekanter a og b, samt avstand z0 fra aapning til obs.skjerm.
a = 0.04e-3;
b = 0.02e-3;
z0 = 1.00;
%% Lager vektorer med verdier for x0 og y0, som representerer
%% et punkt paa observasjonsskjermen. Bruker ogsaa her m.

% Velg området -2 < x < 2 og -2 < z < 2 for plotting, med skrittlengde
% 2/50 langs både x og z, dvs ca 100 x 100 punkter i alt
%[x,z] = meshgrid(-2:2/50:2,-2:2/50:2);
% Ve er det eksakte potensialet fra dipolen
%Ve = (1./(sqrt(x.^2+(z-1/2).^2))-1./(sqrt(x.^2+(z+1/2).^2)));
% Va er dipolpotensialet til ledende orden i a/r når r >> a
%Va = z.*(x.^2.+z.^2).^(-3/2);
% Relavvik er relativ differanse mellom Ve og Va i prosent
%Relavvik = 100*abs((Ve-Va)./(Ve));
% Kommandoen figure ber om en ny figur
%figure;
% Kommandoen mesh(x,y,z) tegner opp et 3D-plott av z som funksjon
% av x og y
%mesh(x,z,Ve);

[x0,y0] = meshgrid(-0.02:0.001001:0.02,-0.02:0.001001:0.02);
%% Lager saa matrise med R-verdier. R = sqrt(x0^2+y0^2+z0^2)
%% R = sqrt(x0.^2 + y0.^2 + z0^2);
%% Nei, vi setter simpelthen R = z0
R = z0;
%% Deretter lager vi matrise I_rekt med funksjonen
%% [sin(alfa)/alfa]^2*[sin(beta)/beta]^2.
%% Her er alfa = pi*x0*a/lambda*R
%% Her er beta = pi*y0*b/lambda*R
alfa = x0*(pi*a/(lambda*R));
beta = y0*(pi*b/(lambda*R));
I_rekt = (sinc(alfa).*sinc(beta)).^2;
I_rekt_log = log(I_rekt);
%% Deretter plotter vi denne intensitetsfordelingen som funksjon av 
%% x0 og y0:
surfc(x0,y0,I_rekt);
axis([-0.02 0.02 -0.02 0.02 -2 2.2])
figure;
mesh(x0,y0,I_rekt_log);
title('Rekt aapning, a=0.04 mm, b = 0.02 mm');
xlabel('x0 (m)');
ylabel('y0 (m)');
zlabel('Normert intensitet');
