close all;
clear;
clc;
%definizione variabili
%grandezze indipendenti dal mezzo
global omega f c T epsilonzero muzero EY eta0
c = 3*10^8;
f = 3*10^8;
omega = 2*pi*f;
T=1/f;
epsilonzero = 8.857*10^-12;
muzero = 4*pi*10^-7;
EY=2;
eta0 = sqrt(muzero/epsilonzero);
%grandezze carattestiche del mezzo e dipendenti dal mezzo
global epsilonr1 epsilonr2 eta1 eta2 lambda1 lambda2 beta1 beta2 v1 v2
epsilonr1 = 5;
epsilonr2 = 1;
eta1 = eta0/sqrt(epsilonr1);
eta2 = eta0/sqrt(epsilonr2);
v1 = c/sqrt(epsilonr1);
v2 = c/sqrt(epsilonr2);
lambda1 = v1/f;
lambda2 = v2/f;
beta1 = (2*pi)/lambda1;
beta2 = (2*pi)/lambda2;
%inizializzazione degli angoli
global thetai thetacritico thetat
thetai = 0.3;
thetacritico=asin(sqrt((epsilonr2)/(epsilonr1)));
%il seguente if assegna a thetat il valore pi/2 se l'angolo d'incidenza è maggiore o uguale all'angolo critico
%invece assegna a thetat il valore ricavato dalla legge di Snell per la rifrazione
if(thetai < thetacritico)
thetat = asin((sqrt(epsilonr1/epsilonr2))*sin(thetai));
else
thetat=(pi/2);
end
%impedenze dei mezzi e coefficienti di riflessione e trasmissione
global Z1te Z2te ROte TAUte
Z1te = eta1/(cos(thetai));
Z2te = eta2/(cos(thetat));
ROte = (Z2te-Z1te)/(Z2te+Z1te);
if(thetai >= thetacritico)
TAUte = 0;
else
TAUte = (2*Z2te)/(Z2te+Z1te);
end
%y0i è l'ordinata all'origine all'istante 0 delle rette che rappresentano l'intersezione tra fronti d'onda e il piano XZ
%passo è quanto si incrementa ad ogni frame dell'animazione
global y0i passo passo2
y0i = -3-(sin(thetai)+cos(thetai)*cot(thetai))*23*(lambda1/10);
passo = (sin(thetai)+cos(thetai)*cot(thetai))*(lambda1/10);
passo2 = (sin(thetai)+cos(thetai)*cot(thetai))*(lambda2/10);
global k %contatore
k=0;
% definizione dei vettori e della matrice che verrano passati alla funzione pcolor
n = 100;
x = linspace(3,-3,n);
z = linspace(-3,3,n);
Ci =zeros(n,n);
if(thetai > 0 && thetai < (pi/2))
while(1)
for i = 1:n
for j = 1:n
%questo if fa in modo che il campo venga disegnato solo per valori della retta sopra l'ordinata all'origine
%dopo 23 frame l'onda collide con l'interfaccia
%ad ogni frame "fotografiamo" l'onda dopo T/10
if(rettai(z(j)) > -3+passo*(k-23))
Ci(i,j) = heaviside(rettai(z(j))-x(i))*Ei(x(i),z(j),k*(T/10))+ (heaviside(rettar(z(j))-x(i))*Er(x(i),z(j),k*(T/10)));
else
if(rettat(z(j))>x(i))
Ci(i,j) = Et(x(i),z(j),k*(T/10));
end
end
end
end
%istruzioni per disegnare l'animazione
si = pcolor(z,x,Ci);
si.FaceColor= 'interp';
si.EdgeColor = 'none';
colormap jet
colorbar
xlabel('Asse z')
ylabel('Asse x')
title('Animazione incidenza obliqua polarizzazione TE');
line([0 0], [-3;3],'color', 'black','LineWidth',2); %retta verticale che simula l'interfaccia
F(k+1) = getframe;
%incremento del contatore
k = k+1;
end
end
%corpo istruzioni che regola il caso particolare thetai=0
if(thetai == 0)
while(1)
for i = 1:n
for j = 1:n
Ci(i,j) =heaviside(-(z(j))-23*(lambda1/10)+k*(lambda1/10))*Ei(x(i),z(j),k*(T/10))+ heaviside(k-23)*heaviside((z(j)-23*(lambda1/10)+k*(lambda1/10)))*Er(x(i),z(j),k*(T/10));
if(z(j)>0)
Ci(i,j) = heaviside(k-23)*heaviside(-(z(j))+(k-23)*(lambda2/10))*Et(x(i),z(j),k*(T/10));
end
end
end
%istruzioni per disegnare l'animazione
si = pcolor(z,x,Ci);
si.FaceColor= 'interp';
si.EdgeColor = 'none';
colormap jet
colorbar
xlabel('Asse z')
ylabel('Asse x')
title('Animazione incidenza obliqua polarizzazione TE');
line([0 0], [-3;3],'color', 'black','LineWidth',2); %retta verticale che simula l'interfaccia
F(k+1) = getframe;
%incremento del contatore
k = k+1;
end
end
%le funzioni che verranno definite ora vengono utilizzate per modellare
%l'avanzamento del primo fronte d'onda delle onde incidente, riflessa e
%trasmessa mediante delle rette che traslano nel tempo
%retta che trasla che descrive l'intersezione tra il primo fronte dell'onda incidente e il piano XZ
function y = rettai(x)
global thetai k y0i passo
y = -x*cot(thetai)+y0i+passo*k;
end
%retta che trasla che descrive l'intersezione tra il primo fronte dell'onda riflessa e il piano XZ
function y = rettar(x)
global thetai k passo
y = x*cot(thetai)-3+passo*(k-23);
end
%retta che trasla che descrive l'intersezione tra il primo fronte dell'onda trasmessa e il piano XZ
function y = rettat(x)
global thetat k passo
y = -x*cot(thetat)-3+passo*(k-23);
end
%espressione del campo elettrico dell'onda incidente
function A = Ei(x,z,t)
global EY thetai beta1 omega
A = EY*(cos(omega*t-beta1*(x*sin(thetai)+z*cos(thetai))));
end
%espressione del campo elettrico dell'onda riflessa
function A = Er(x,z,t)
global EY thetai beta1 ROte omega
A =ROte* EY*(cos(omega*t-beta1*(x*sin(thetai)-z*cos(thetai))));
end
%espressione del campo elettrico dell'onda trasmessa
function A = Et(x,z,t)
global EY thetat beta2 TAUte omega
A =TAUte* EY*(cos(omega*t-beta2*(x*sin(thetat)+z*cos(thetat))));
end
