Guida dielettrica a 5 strati generica: studio completo. Programma principale

Next: Studio completo: function Up: Introduzione Previous: Curve di dispersione: esecuzione
Intro Gen: Introduzione Generale Home: Home page

CORSO DI LAUREA IN INGEGNERIA ELETTRONICA

PROGRAMMI IN MATLAB

DI CAMPI ELETTROMAGNETICI

ANNO ACCADEMICO 1997 / 1998

Guida dielettrica a 5 strati generica:

Studio completo e andamenti del campo elettrico nello spazio

Professore : Marco De Sario

Eseguito da Marinelli Vito

PROGRAMMA PRINCIPALE

guida.m

% GUIDA DIELETTRICA A 5 STRATI : ANDAMENTO DEL CAMPO ELETTRICO

clc;

clear all;

gcf;

for figura=1:gcf

close;

end

clear global;

format long e;

whitebg;

disp('Introduzione delle dimensioni dei vari strati :');

disp('Introduci lo spessore (in metri): ');

c2=input('del primo cuore : ');

g=input('del gap : ');

c4=input('del secondo cuore : ');

disp('Introduci l''indice di rifrazione : ');

n(1)=input('del superstrato : ');

n(2)=input('del primo cuore : ');

n(3)=input('del gap : ');

n(4)=input('del secondo cuore : ');

n(5)=input('del sottostrato : ');

l=input('Introduci la lunghezza d''onda (in metri): ');

if n(2)>n(4) % i modi che calcoliamo devono essere

k=n(2); % guidati almeno nel secondo cuore

n(2)=n(4);

n(4)=k;

k=c2;

c2=c4;

c4=k;

k=n(1);

n(1)=n(5);

n(5)=k;

end

k0=2*pi/l;

e0=8.85418781762e-12;

y0=pi*4e-7;

d2=g/2+c2;

d4=g/2+c4;

errore=eps;

j=1;

passo=-(n(4)-n(3))/1000;

% calcolo degli zeri della matrice di dispersione ------------------

for nef=n(4):passo:n(3)-passo

for k=1:5

h(k)=k0*sqrt(abs(n(k)^2-nef^2));

end

pol2=mdisp(h,d2,d4,g,nef,n);

if nef~=n(4) & nef~=n(4)+passo

if nef>n(2) | nef-passo<n(2)

if pol2*pol1<=0

if pol2==0

zero(j)=nef;

j=j+1;

else

% metodo di bisezione ------------------------------

indice=(log((2*nef-passo)/errore))/log(2);

pol11=pol1;

pol21=pol2;

ak=nef;

bk=nef-passo;

indice=ceil(indice);

ciclo=1;

while ciclo<=indice

medio=(ak+bk)/2;

for k=1:5

hmedio(k)=k0*sqrt(abs(n(k)^2-medio^2));

end

determ=mdisp(hmedio,d2,d4,g,medio,n);

if determ==0

ciclo=indice;

else

if determ*pol11<0

ak=medio;

pol21=determ;

else

bk=medio;

pol11=determ;

end

ciclo=ciclo+1;

end

%--------------------------------------------------------------------------------

zero(j)=medio;

j=j+1;

end

pol1=pol2;

end

%------------------------------------------------------------------------------------------------

disp('I possibili valori dell''indice efficace sono :');

nef=zero

lunghezza=length(zero);

i=1;

str='Per l'' indice di rifrazione efficace pari a ';

while i<=lunghezza

clear A b B h coeff Ey x

for k=1:5

h(k)=k0*sqrt(abs(n(k)^2-nef(i)^2));

end

% costruzione del sistema con il quale verranno calcolati i vari coefficienti

errore=0;

A(1,1)=exp(-h(1)*d2);

A(2,1)=h(1)*exp(-h(1)*d2);

A(3,4)=-exp(h(3)*g/2);

A(3,5)=-exp(-h(3)*g/2);

if nef(i)>n(2)

A(1,2)=-exp(h(2)*d2);

A(1,3)=-exp(-h(2)*d2);

A(2,2)=h(2)*exp(h(2)*d2);

A(2,3)=-h(2)*exp(-h(2)*d2);

A(3,2)=exp(h(2)*g/2);

A(3,3)=exp(-h(2)*g/2);

A(4,2)=-h(2)*exp(h(2)*g/2);

A(4,3)=h(2)*exp(-h(2)*g/2);

else

A(1,2)=-cos(h(2)*d2);

A(1,3)=-sin(h(2)*d2);

A(2,2)=-h(2)*sin(h(2)*d2);

A(2,3)=h(2)*cos(h(2)*d2);

A(3,2)=cos(h(2)*g/2);

A(3,3)=sin(h(2)*g/2);

A(4,2)=h(2)*sin(h(2)*g/2);

A(4,3)=-h(2)*cos(h(2)*g/2);

end

A(4,4)=h(3)*exp(h(3)*g/2);

A(4,5)=-h(3)*exp(-h(3)*g/2);

A(5,4)=exp(-h(3)*g/2);

A(5,5)=exp(h(3)*g/2);

A(5,6)=-cos(h(4)*g/2);

A(5,7)=sin(h(4)*g/2);

A(6,4)=-h(3)*exp(-h(3)*g/2);

A(6,5)=h(3)*exp(h(3)*g/2);

A(6,6)=h(4)*sin(h(4)*g/2);

A(6,7)=h(4)*cos(h(4)*g/2);

A(7,6)=cos(h(4)*d4);

A(7,7)=-sin(h(4)*d4);

A(8,6)=-h(4)*sin(h(4)*d4);

A(8,7)=-h(4)*cos(h(4)*d4);

b(7)=exp(-h(5)*d4);

b(8)=-h(5)*exp(-h(5)*d4);

if rank(A)<7

errore=1;

else

i2=1;

while i2<=8

k2=1;

for j2=1:8

if j2~=i2

B(k2,:)=A(j2,:);

c(k2)=b(j2);

k2=k2+1;

end

if det(B)~=0

i2=8;

A=B;

b=c;

end

i2=i2+1;

end

%--------------------------------------------------------------------------------------------

neff=num2str(nef(i));

if errore==1

disp([str neff]);

disp('il problema ammette infinite soluzioni.');

disp('Sicuramente i dati immessi sono errati.Ricomincia da capo');

i=lunghezza;

else

x=risoluz(A,b);

disp([str neff]);

disp('se normalizziamo i campi in ampiezza ponendo quindi B5=1,');

disp('si hanno i seguenti coefficienti incogniti :');

A1=x(1)

B2=x(2)

A2=x(3)

B3=x(4)

A3=x(5)

B4=x(6)

A4=x(7)

B5=1

% calcolo del coefficiente B5 -------------------------------

for k=1:5

coeff(k)=(1/k0)*sqrt(e0/(y0*abs(1-(n(k)/nef(i))^2)));

end

I1=(coeff(5)/4)*exp(-2*h(5)*d4);

I21=(coeff(4)/4)*((B4^2)-(A4^2))*sin(h(4)*c4)*cos(h(4)*(c4+g));

I22=-(coeff(4)/2)*A4*B4*sin(h(4)*c4)*sin(h(4)*(c4+g));

I23=((A4^2)+(B4^2))*c4*coeff(4)*h(4)/4;

I3=coeff(3)*((1/4)*((A3^2)+(B3^2))*(exp(h(3)*g)-exp(-h(3)*g))+h(3)*A3*B3*g);

if nef(i)>n(2)

I41=(coeff(2)/4)*(A2^2)*(exp(-h(2)*g)-exp(-2*h(2)*d2));

I42=(coeff(2)/4)*(B2^2)*(exp(2*h(2)*d2)-exp(h(2)*g));

I43=A2*B2*c2*coeff(2)*h(2);

else

I41=(coeff(2)/4)*((B2^2)-(A2^2))*sin(h(2)*c2)*cos(h(2)*(c2+g));

I42=(coeff(2)/2)*A2*B2*sin(h(2)*c2)*sin(h(2)*(c2+g));

I43=((A2^2)+(B2^2))*c2*coeff(2)*h(2)/4;

end

I5=coeff(1)*A1*A1*exp(-2*h(1)*d2)/4;

Potenza=I1+I21+I22+I23+I3+I41+I42+I43+I5;

B5=1/sqrt(Potenza);

%-----------------------------------------------------------------------------------------

disp('Con i quali la potenza dell''onda é pari a : ');

Potenza

disp('Adesso normalizzando i campi in potenza utilizzando tali coefficienti');

strn='si ottiene un valore di B5=';

b5=num2str(B5);

disp([strn b5]);

disp('Per cui gli altri coefficienti adesso assumono il seguente valore : ');

A1=A1*B5

B2=B2*B5

A2=A2*B5

B3=B3*B5

A3=A3*B5

B4=B4*B5

A4=A4*B5

% Andamento grafico dei vari modi per z=0 -------------------

passo=2*(d2+d4)/2000;

figure;

hold off;

punti=1;

punto1=0;

for x=-2*d4:passo:-d4

punto1=punto1+1;

y(punto1)=x;

xx(punti)=x;

Ey(punti)=B5*exp(h(5)*x);

punti=punti+1;

end

plot(Ey,y,'k');

clear y

y(1)=x;

xx(punti)=x;

Ey(punti)=Ey(punti-1);

punto2=1;

for x=-d4:passo:-g/2

punto2=punto2+1;

y(punto2)=x;

xx(punti)=x;

Ey(punti)=B4*cos(h(4)*x)+A4*sin(h(4)*x);

punti=punti+1;

end

hold on;

z=Ey(punti-punto2:punti-1);

plot(z,y,'r');

clear y z

y(1)=x;

xx(punti)=x;

Ey(punti)=Ey(punti-1);

punto3=1;

for x=-g/2:passo:g/2

punto3=punto3+1;

y(punto3)=x;

xx(punti)=x;

Ey(punti)=B3*exp(h(3)*x)+A3*exp(-h(3)*x);

punti=punti+1;

end

hold on;

z=Ey(punti-punto3:punti-1);

plot(z,y,'b');

clear y z

y(1)=x;

xx(punti)=x;

Ey(punti)=Ey(punti-1);

punto4=1;

for x=g/2:passo:d2

punto4=punto4+1;

y(punto4)=x;

xx(punti)=x;

if nef(i)>n(2)

Ey(punti)=B2*exp(h(2)*x)+A2*exp(-h(2)*x);

else

Ey(punti)=B2*cos(h(2)*x)+A2*sin(h(2)*x);

end

punti=punti+1;

end

hold on;

z=Ey(punti-punto4:punti-1);

plot(z,y,'m');

clear y z

y(1)=x;

xx(punti)=x;

Ey(punti)=Ey(punti-1);

punto5=1;

for x=d2:passo:2*d2

punto5=punto5+1;

y(punto5)=x;

xx(punti)=x;

Ey(punti)=A1*exp(-h(1)*x);

punti=punti+1;

end

hold on;

z=Ey(punti-punto5:punti-1);

plot(z,y,'k');

clear y z

punti=punti-1;

%-----------------------------------------------------------------------------------------

for ind=1:length(xx)

Eyy(ind,i)=Ey(ind);

end

i=i+1;

xlabel('Ey');

ylabel('x');

grid on;

end

if errore==0

risp=input('Se vuoi disegnare il grafico tridimensionale premi "s" : ','s');

% Costruzione dei grafici tridimensionali -----------------------

if risp=='s' | risp=='S'

clear Ey B i

x=xx;

dimensione=length(xx);

if n(2)==n(4) & n(1)==n(5) & c2==c4 & lunghezza==2

zi=l/(2*(nef(1)-nef(2)));

else

zi=input('Introduci il valore della quota z (in metri): ');

end

figure;

hold off;

i=sqrt(-1);

punto2=punto1+punto2-1;

punto3=punto2+punto3-1;

punto4=punto3+punto4-1;

punto5=dimensione;

if n(2)==n(4) & n(1)==n(5) & c2==c4 & lunghezza==2

disp('L''andamento nello spazio del modulo del campo elettrico é il seguente:');

pass=0;

for z1=0:zi/10:zi

pass=pass+1;

indd=0;

for ind=1:30:dimensione

indd=indd+1;

Ey(indd,pass)=0;

for ii=1:lunghezza

Ey(indd,pass)=Ey(indd,pass)+Eyy(ind,ii)*exp(-i*k0*nef(ii)*z1);

end

Ey(indd,pass)=abs(Ey(indd,pass));

end

z=0:zi/10:zi;

map=[0.35,0.35,0.35;0,0,1;1,0,0;1,.5,0];

colormap(map);

for ind1=1:pass

for ind2=1:punto1

C(ind2,ind1)=100;

end

for ind2=punto1+1:punto2

C(ind2,ind1)=200;

end

for ind2=punto2+1:punto3

C(ind2,ind1)=150;

end

for ind2=punto3+1:punto4

C(ind2,ind1)=300;

end

for ind2=punto4+1:punto5

C(ind2,ind1)=100;

end

C=C(1:30:dimensione,:);

mesh(x(1:30:dimensione),z,Ey',C');

xlabel('x');

ylabel('z');

zlabel('|Ey(x,z)|');

grid on;

figure;

hold off;

end

disp('Alcuni di questi andamenti nello spazio sono messi in risalto in questo grafico :');

clear Ey;

for z1=0:zi/4:zi

for ind=1:dimensione

Ey(ind)=0;

for ii=1:lunghezza

Ey(ind)=Ey(ind)+Eyy(ind,ii)*exp(-i*k0*nef(ii)*z1);

end

Ey(ind)=abs(Ey(ind));

z(ind)=z1;

end

if z1==0 | z1==zi

plot3(x(1:punto1),z(1:punto1),Ey(1:punto1),'k');

hold on;

plot3(x(punto1:punto2),z(punto1:punto2),Ey(punto1:punto2),'r');

hold on;

plot3(x(punto2:punto3),z(punto2:punto3),Ey(punto2:punto3),'b');

hold on;

plot3(x(punto3:punto4),z(punto3:punto4),Ey(punto3:punto4),'r');

hold on;

plot3(x(punto4:punto5),z(punto4:punto5),Ey(punto4:punto5),'k');

hold on;

else if z1==zi/4

plot3(x,z,Ey,'m');

hold on;

end

if z1==2*zi/4

plot3(x,z,Ey,'c');

hold on;

end

if z1==3*zi/4

plot3(x,z,Ey,'g');

hold on

end

clear Ey

end

xlabel('x');

ylabel('z');

zlabel('|Ey(x,z)|');

grid on;

%---------------------------------------------------------------------------------------------

end

clear all

Next: Studio completo: function Up: Introduzione Previous: Curve di dispersione: esecuzione
Intro Gen: Introduzione Generale Home: Home page
Vito Marinelli
13-6-2000

BPath Contatore