19
Relatório: Fotónica 1 1. Resolução numérica de equações diferenciais ordinárias Na simulação numérica da modulação directa da corrente de injecção em lasers semicondutores há que resolver as equações das taxas. Considera-se, aqui, o caso dos lasers semicondutores monomodais. Estas equações regulam a evolução temporal, na zona activa do laser semicondutor, das populações totais de electrões N (envolvidos nos processos de geração e recombinação) e de fotões S , tendo-se equações das taxas st sp a p t dS R R R dt dN R R dt = Γ + = em que 1 Γ< representa o factor de confinamento óptico (o modo de oscilação laser não fica exactamente confinado à zona activa). As várias taxas intervenientes têm unidades ( ) 1 s e são: (i) a taxa de bombeamento p R ; (ii) a taxa total de recombinação (de pares electrão- lacuna) t R ; (iii) a taxa efectiva líquida de emissão estimulada st R ; (iv) a taxa contributiva de emissão espontânea sp R; (v) a taxa efectiva de aniquilação de fotões a R (relacionada com o tempo de vida médio p τ dos fotões na cavidade laser). Shun Lien Chuang, Physics of Optoelectronic Devices (New York: Wiley, 1995), Chapter 11 (pp. 487-507). Tem-se

Mat Lab Drog

  • Upload
    ad0001

  • View
    3

  • Download
    2

Embed Size (px)

DESCRIPTION

MATLAB

Citation preview

Page 1: Mat Lab Drog

Relatório: Fotónica 1

1. Resolução numérica de equações diferenciais ordinárias Na simulação numérica da modulação directa da corrente de injecção em lasers

semicondutores há que resolver as equações das taxas. Considera-se, aqui, o caso dos lasers

semicondutores monomodais. Estas equações regulam a evolução temporal, na zona activa do

laser semicondutor, das populações totais de electrões N (envolvidos nos processos de

geração e recombinação) e de fotões S , tendo-se

equações das taxasst sp a

p t

d S R R Rdtd N R Rdt

′= Γ + −→

= −

em que 1Γ < representa o factor de confinamento óptico (o modo de oscilação laser não fica

exactamente confinado à zona activa). As várias taxas intervenientes têm unidades ( )1s− e

são: (i) a taxa de bombeamento pR ; (ii) a taxa total de recombinação (de pares electrão-

lacuna) tR ; (iii) a taxa efectiva líquida de emissão estimulada stR ; (iv) a taxa contributiva de

emissão espontânea spR′ ; (v) a taxa efectiva de aniquilação de fotões aR (relacionada com o

tempo de vida médio pτ dos fotões na cavidade laser).

Shun Lien Chuang, Physics of Optoelectronic Devices (New York: Wiley, 1995),

Chapter 11 (pp. 487-507).

Tem-se

Page 2: Mat Lab Drog

2 Carlos R. Paiva

, , , ,p t st st sp sp sp ap

I NR R R R R G S R R Rq

βτ

′= = + = = =

em que I representa a corrente de injecção (sendo q o valor da carga do electrão,

considerada aqui como positiva), G a taxa elementar líquida de emissão estimulada e spβ o

coeficiente de emissão espontânea contributiva (i.e., a pequena parcela da taxa de emissão

espontânea total spR que contribui para a emissão estimulada). A taxa R , que integra a taxa

total de recombinação tR , representa a taxa de aniquilação de electrões e é dada por

c

NRτ

= .

Em geral cτ , que é o tempo de vida médio dos electrões, depende de N . Para simplificar a

análise considera-se, frequentemente, que cτ é uma constante. A taxa spR′ relaciona-se com

G através da expressão

sp sp sp spR R n Gβ′ = =

em que spn é o chamado coeficiente de inversão da população. As equações das taxas

assumem, portanto, a forma explícita

equações das taxassp

p

c

d S SG S n Gdt

d N I NG Sdt q

τ

τ

= Γ + −

→= − −

e, para serem resolvidas, necessitam que se adopte um modelo que relacione G com N e S .

Um modelo possível (não linear) é

( )1

N tG N NG

Sε−

=+

Page 3: Mat Lab Drog

Relatório: Fotónica 3

onde NG , tN e ε são constantes que caracterizam este modelo específico. A resolução

numérica das equações das taxas enquadra-se, portanto, no quadro mais geral da resolução

numérica de equações diferenciais ordinárias.

Apresenta-se, a seguir, um caso concreto de um sistema de duas equações diferenciais

ordinárias que pertence à mesma classe das equações das taxas. Considere-se, com efeito, que

( ) ( )

( ) ( )

( ) ( )( ) ( )

1 111 2

1 2

2 222 1

1 2

11

11

B p zd p A p zd z p z p z

B p zd p A p zd z p z p z

⎧= −⎡ ⎤⎪ ⎣ ⎦ + +⎪

⎨⎪ = −⎡ ⎤⎣ ⎦⎪ + +⎩

onde 1 1 2 2, , ,A B A B ∈ e 1 2, :p p → são funções conhecidas. Note-se que, neste

problema, a coordenada espacial z desempenha o papel do tempo t nas equações das taxas.

Para a resolução deste novo sistema é necessário conhecer as condições iniciais

( )( )

1 0

2 0

00

0p a

zp b

=⎧⎪= → ⎨ =⎪⎩.

Apresenta-se, então, um programa MATLAB que permite resolver este novo problema.

1

02

01

2

1101constantes12escolhidas

1

AaAbB

B

= −==

→==

=

[ ]intervalo de pesquisa 0,10z→ ∈

Programa MATLAB: «diferen.m».

Nota: O programa «diferen» chama, durante a sua execução, a subrotina «derv» através do

programa «ode45».

Page 4: Mat Lab Drog

4 Carlos R. Paiva

Programa «diferen.m»

% DIFEREN Resolução do seguinte sistema de equações diferenciais

%

% d B1 * p1(z)

% ---- p1(z) = ----------------------- * [A1 * p2(z) - 1]

% d z 1 + p1(z) + p2(z)

%

% d B2 * p2(z)

% ---- p2(z) = -------------------- * [A2 * p1(z) - 1]

% d z 1 + p1(z) + p2(z)

%

clear all % Limpa todas as variáveis da memória

close all % Fecha todas as janelas de figuras

A1 = -1; A2 = 1; B1 = 2; B2 = 1; % Constantes do sistema de equações

z=linspace(0,10,1000); % Vector distância onde z = 0 até z = 10 com 1000 pontos

pini=[10 1]; % Valores iniciais das variaveis p1(z) = 10 e p2(z) = 1 no ponto z = 0

OPTIONS=odeset('AbsTol',1e-9,'RelTol',1e-6); % Alteração do valor da tolerância do

% método de RUNGE-KUTTA

% (função ODE45.m)

% Método de RUNGE-KUTTA. 'derv' é a função onde está descrito o sistema diferencial

[z,p]=ode45('derv',z,pini,OPTIONS,A1,A2,B1,B2);

plot(z,p(:,1),'b-',z,p(:,2),'r-'); % Gráfico de p1(z) e p2(z)

title('p1(z) - Azul p2(z) - Vermelho');

xlabel('Distância');

ylabel('Amplitude');

Page 5: Mat Lab Drog

Relatório: Fotónica 5

Função «derv.m»

% DERV Função executada por ode45.m com a descrição do sistema de equações diferenciais

% z - posição onde calcular o valor das derivadas de p1(z) e p2(z)

% p - vector com os valorer p1(z) e p2(z) na posição z

% flag - flags usadas pela função ode45.m

% A1, A2, B1, B2 - constantes do sistema diferencial

%

function dp=derv(z,p,flag,A1,A2,B1,B2)

dp(1)=B1*p(1)/(1+p(1)+p(2))*(A1*p(2)-1); % Equação dp1(z)/dz

dp(2)=B2*p(2)/(1+p(1)+p(2))*(A2*p(1)-1); % Equação dp2(z)/dz

dp=dp'; % Passagem do vector das derivadas de linha para coluna

Figura 1 – Resultado obtido através do programa DIFEREN.

William E. Boyce and Richard C. DiPrima, Elementary Differential Equations and

Boundary Value Problems (Hoboken, New Jersey: Wiley, 8th ed., 2005).

Page 6: Mat Lab Drog

6 Carlos R. Paiva

Apresentam-se, de seguida, alguns gráficos correspondentes à simulação numérica das

equações das taxas num laser semicondutor com modulação directa da corrente de injecção.

Page 7: Mat Lab Drog

Relatório: Fotónica 7

Figura 2 – Alguns resultados do Trabalho T1 – alínea (b).

Page 8: Mat Lab Drog

8 Carlos R. Paiva

2. Resolução numérica de equações modais Em vários pontos do programa de Fotónica é necessário resolver numericamente equações

modais. É o caso dos lasers semicondutores, onde é necessário determinar – usando o método

do índice de refracção efectivo – os valores de effn e de n , ver Trabalho T3. É ainda o caso

das fibras ópticas, onde é necessário resolver a respectiva equação modal – seja de forma

rigorosa (modos HE, EH, TE e TM), seja de forma aproximada (modos LP), ver Trabalhos T4

e T5. Em qualquer dos casos a utilização do programa MATLAB «fzero» é perfeitamente

suficiente. O primeiro exemplo que se vai apresentar refere-se ao Trabalho T2 sobre mecânica

quântica. Neste trabalho considera-se um modelo de Kronig-Penney simplificado para

potenciais periódicos. O objectivo é o de entender o aparecimento de bandas de energia

permitidas e proibidas neste tipo de potenciais – uma base indispensável para a física do

estado sólido. Na Fig. 3 apresenta-se um gráfico de (ver Trabalho T2)

( ) ( ) 2 2 2

0 20

,2

k aq a ma

ξ ζ ζ ξ πζ π= ⎡ ⎤

→ = =⎢ ⎥= ⎣ ⎦

EE

E.

Figura 3 – Resultados da segunda parte do Trabalho T2.

Page 9: Mat Lab Drog

Relatório: Fotónica 9

Nas Figs. 4 e 5 apresentam-se os resultados reslativos ao trabalho T3. Nas Figs. 6 e 7

apresentam-se os resultados relativos aos trabalhos T4 e T5 (respectivamente).

Figura 4 – Resultados da primeira fase do Trabalho T3.

George Lindfield and John Penny, Numerical Methods Using MATLAB (London: Ellis

Horwood, 1995).

Page 10: Mat Lab Drog

10 Carlos R. Paiva

Figura 5 – Resultados da segunda fase do Trabalho T3.

Figura 6 – Resultados do Trabalho T4.

Page 11: Mat Lab Drog

Relatório: Fotónica 11

Figura 7 – Resultados do Trabalho T5.

Amnon Yariv and Pochi Yeh, Photonics: Optical Electronics in Modern

Communications (New Yor: Oxford University Press, 6th ed., 2007), pp. 110-155, pp.

797-811.

Bahaa E. A. Saleh and Marvin Carl Teich, Fundamentals of Photonics (Hoboken, New

Jersey: Wiley, 2nd ed., 2007), pp. 289-364.

Page 12: Mat Lab Drog

12 Carlos R. Paiva

3. Propagação de impulsos em fibras ópticas Para a simulação numérica da propagação de impulsos em fibras ópticas há que recorrer à

FFT (fast Fourier transform). Consideremos apenas o caso das fibras monomodais operadas

em regime linear. Despreza-se, aqui, a influência da dispersão de ordem superior (i.e.,

considera-se que é razoável fazer 3 0β = ). Assim, só o efeito da DVG (dispersão da

velocidade de grupo) é tido em consideração. As perdas são desprezadas de forma a poder

analisar, isoladamente, o efeito da DVG sobre a propagação os impulsos. Nestas

circunstâncias e em termos das variáveis normalizadas ζ e τ

220

1

0

D

z zL

t z

βζ

τ

βττ

= =

−=

a equação de propagação dos impulsos é dada pela equação diferencial

( )2

2 2

1 sgn 02

u ui βζ τ∂ ∂

− =∂ ∂

.

Consideremos o caso específico da zona de dispersão anómala em que

( )2

zona de dispersãosgn 1

anómalaβ→ = −

de forma que a equação de propagação dos impulsos se reduz a

2

2

12

u uiζ τ∂ ∂

= −∂ ∂

.

Esta equação é facilmente resolvida no domínio da frequência (normalizada) ξ

( ) ( ) ( ), , ,u u uζ τ ζ ξ ζ τ= ⎡ ⎤⎣ ⎦F

Page 13: Mat Lab Drog

Relatório: Fotónica 13

( ) ( ) ( )

( ) ( ) ( )

, , exp

1, , exp2

u u i d

u u i d

ζ ξ ζ τ ξ τ τ

ζ ξ ζ τ ξ τ ξπ

−∞

−∞

=

= −

uma vez que

( )21 ,2

u i uξ ζ ξζ∂

=∂

( ) ( ) 21, 0, exp2

u u iζ ξ ξ ξ ζ⎛ ⎞∴ = ⎜ ⎟⎝ ⎠

.

Para determinar a forma do sinal, após uma determinada distância de propagação ao longo da

fibra óptica, há que voltar para o domínio do tempo – aplicando agora uma IFFT (a inversa de

uma FFT). No programa que a seguir se apresenta, intitulado «propag», apresenta-se o efeito

da DVG sobre um impulso com uma forma inicial

( ) ( ) ( )0 0, sechu uτ τ τ= = .

Sublinhe-se a necessidade de alterar a ordem interna do vector das frequências devido à forma

como a FFT apresenta a função no domínio da frequência. Este mesmo efeito poderia ser

alcançado através da instrução «fftshift» que, em vez de alterar a ordem interna do vector das

frequências, altera a própria disposição da transformada (colocando o valor correspondente à

frequência nula no centro). Em qualquer caso é necessário garantir que a multiplicação pela

função de transferência da fibra óptica se processa, no domínio da frequência, de forma

correcta.

Amnon Yariv and Pochi Yeh, Photonics: Optical Electronics in Modern

Communications (New Yor: Oxford University Press, 6th ed., 2007), pp. 317-322.

Page 14: Mat Lab Drog

14 Carlos R. Paiva

% PROPAG Propagação de impulsos regida pela equação diferencial

%

% d i d ^ 2

% ---- a(z, t) = --- ---------- a(z,t)

% d z 2 d t ^ 2

%

clear all; % Limpa todas as variáveis da memória

ztotal=5; % Distância total considerada

Nz=500; % Número de passos na distância

z=linspace(0,ztotal,Nz); % Vector distância com 'Nz' posições

t0=20; % Limite da janela temporal

Nt=1024; % Número de amostras temporais

t=linspace(-t0,t0,Nt); % Vector temporal com 'Nt' pontos (Janela temporal)

Ts=t(2)-t(1); % Separação entre amostras

Ws=2*pi/Ts; % Largura total da janela espectral

W=Ws*[0:1:Nt/2-1 -Nt/2:1:-1]/Nt; % Vector de frequências [rad/s]

a=sech(t); % Impulso inicial: a(0,t)=sech(t)

A=fft(a); % Impulso no domínio espectral: A(0,w)

for i=2:Nz % Ciclo para cada posição de z

Az=A.*exp(-1i/2*W.^2*z(i)); % Impulso no domínio espectral: A(z,w)

a(i,:)=ifft(Az); % Impulso no domínio temporal: a(z,t)

end

figure(1); % Gráfico 2D do impulso inicial/final

plot(t,a(1,:),'b',t,abs(a(end,:)),'r');

title('entrada - azul saída - vermelho');

xlabel('Tempo');

ylabel('Amplitude');

axis([-t0 t0 0 1]);

[X,Y]=meshgrid(t,z); % Cria uma grelha rectangular de dimensões (Nt,Nz) usada por 'mesh'

figure(2);

Page 15: Mat Lab Drog

Relatório: Fotónica 15

mesh(X,Y,abs(a)); % Gráfico da evolução do impulso

xlabel('Tempo');

ylabel('Distância');

zlabel('Amplitude');

figure(3);

mesh(X,Y,abs(a));

xlabel('Tempo');

ylabel('Distância');

zlabel('Amplitude');

view(30,45); % view (azimute, elevação);

Figura 8 – Primeira figura obtida através do programa PROPAG.

Page 16: Mat Lab Drog

16 Carlos R. Paiva

Figura 9 – Segunda figura obtida através do programa PROPAG.

Figura 10 – Terceira figura obtida através do programa PROPAG.

Page 17: Mat Lab Drog

Relatório: Fotónica 17

4. Propagação de solitões em fibras ópticas Para estudar a propagação de impulsos numa fibra óptica em regime não linear usa-se o

SSFM (split-step Fourier method). Em termos das variáveis normalizadas ζ e τ , tem-se

então

( ) ( )2

22 2

equação de propagação de impulsos 1 sgnregime não linear 2 2

u ui u u i uβζ τ∂ ∂ Γ

→ − + = −∂ ∂

que se reduz à chamada equação NLS (não linear de Schrödinger) quando ( )2sgn 1β = −

(zona de dispersão anómala) e 0Γ = (sem perdas). Note-se que os efeitos não lineares e

dispersivos de ordem superior são, aqui, desprezados. O SSFM baseia-se na separação entre

os efeitos não lineares (operador N ) e os efeitos dispersivos (operadores Dτ ou Dξ ):

( ) ( )( )

( )

2

2 22

22

1 sgn2 1, sgn

22

D iu D N u D i

N i u

τ

τ ξ

βτ

ζ τ β ξζ

⎧ ∂=⎪ ∂∂ ⎪= + → → = −⎨∂ Γ⎪ = − +⎪⎩

( ) ( ) ( ), , ,u u uζ τ ζ ξ ζ τ= ⎡ ⎤⎣ ⎦F .

Assim, nesta versão do SSFM, considera-se o processo iterativo de passo h que permite ir do

impulso ( ) ( )0 0,u uτ τ= até à saída ( ),Lu ζ τ , com L DL Lζ = .

( ) ( ) ( )SSFM : , , ,u u h wζ τ ζ τ ζ τ→ + =

( ) ( ) ( )( ) ( )( ) ( ) ( )( ) ( )

, exp ,

, FFT ,

, exp ,

, IFFT ,

v h N u

V v

W h D V

w W

ξ

ζ τ ζ τ

ζ ξ ζ τ

ζ ξ ζ ξ

ζ τ ζ ξ

=

= ⎡ ⎤⎣ ⎦

=

= ⎡ ⎤⎣ ⎦

Page 18: Mat Lab Drog

18 Carlos R. Paiva

% SOL_N programa para representar a evolução de um SOLITÃO de ordem N

% Método: SSFM (Split-Step Fourier Method)

%

clear all;

close all;

N=3; % Ordem N do solitão sob consideração

Ntau=2^11; % Número de pontos da janela temporal

tmax=5; % Limites máximo e mínimo da janela temporal

tau=linspace(-tmax,tmax,Ntau); % Definiçao da janela temporal (em 'tau')

zetaMax=pi/2; % Definiçãoo da distância normalizada 'zeta'

h=0.005; % Passo incremental do SSFM

Nstop=round(zetaMax/h); % Número de iterações do SSFM

zeta=linspace(0,zetaMax,Nstop); % Definição da janela espacial (em 'zeta')

Ts=tau(2)-tau(1); % Discretização temporal

Ws=2*pi/Ts; % Discretização na frequência

xi=Ws*[0:1:Ntau/2-1 -Ntau/2:1:-1]/Ntau; % Vector das frequências normalizadas

u=N*sech(tau); % Definição do impulso inicial

beta2=-1; % Sinal de beta 2 (DVG)

gama=0; % Perdas normalizadas

impulso(1,:)=abs(u); % Módulo do impulso

D=i/2*sign(beta2)*xi.^2; % operador D de dispersão (DVG)

for k=2:Nstop

N=-gama/2+i*abs(u).^2;

v=(u).*exp(h*N);

V=fft(v);

U=V.*exp(h*D);

u=ifft(U);

impulso(k,:)=abs(u);

end

Page 19: Mat Lab Drog

Relatório: Fotónica 19

[XN,YN]=meshgrid(tau,zeta);

mesh(XN,YN,impulso);

axis([-tmax tmax 0 zetaMax 0 max(max(impulso))]);

xlabel('\tau = (t - \beta_1 z) / \tau_0');

ylabel('\zeta = z / L_D');

zlabel('|u| (\zeta, \tau)');

title('Solitão de ordem N = 3')

Figura 11 – Resultado obtido através do programa SOL_N.

Govind P. Agrawal, Nonlinear Fiber Optics (San Diego, California: Academic Press,

4th ed., 2007), pp. 41-45.