% Programa para calcular o Indice efetivo (constantes de propagacao)
% em guias de ondas opticos de tres camadas

clear;
n2=2.0;                 % Indice do nucleo
n3=1.95;				% Indice do substrato
n1=1.8;				    % Indice do superstrato (casca superior)
lambda=0.6328;          % Comprimento de onda (entre com o valor em microns)
ko=2*pi/lambda;         % Numero de onda
tg=1.0;                 % espessura do nucleo  (em microns)

% As variaveis xi e xf (logo abaixo) representam a faixa de variacao do indice 
% efetivo (Beta/ko) onde voce espera encontrar uma raiz.
% Voce pode escolher a faixa de indice efetivo indo de n3+delta ate n2-delta 
% (para este exemplo). Esse delta e apenas para evitar que haja uma divisao por
% zero (delta=1e-6, por exemplo). Uma vez que voce tenha visto a localizacao da raiz,
% e so definir xi e xf em torno desse ponto.

% Experimente utilizar xi=n3+1e-6; e xf=n2-1e-6; Voce vai ver que a raiz estara
% entre 1.985 e 1.99, assim, defina novamente xi=1.985; e xf=1.99 que a raiz
% sera calculada mais precisamente (sem a possibilidade de cair em outro ponto
% qualquer que possa ser um outro modo ou mesmo um ponto de divergencia que nao e
% uma raiz.

% Lembre que este metodo de achar a raiz e apenas ilustrativo, sem compromisso de achar 
% valores exatos para  as mesmas. Se quiser precisao, experimente utilizar o comando "fmin"
% com a mesma faixa de variacao para xi e xf. So que neste caso voce tera que definir sua
% equacao transcendental em uma "function" separada do programa principal.


xi=n3+1e-6;             % valor inicial
xf=n2-1e-6;             % valor final
                        % Veja na Figure 1 a melhor faixa de valores para xi e xf para voce
                        % encontrar a raiz de seu interesse.

stp=(xf-xi)/20000;              % passo (aqui eu defino o tamanho do vetor Beta 
Beta=(xi:stp:xf)*ko;            % Vetor Beta com 20000 pontos
q=sqrt(Beta.^2-(n1*ko)^2);      % Constante de propagacao do meio 1 (seria o seu k1)
h=sqrt((n2*ko)^2-Beta.^2);      % Constante de propagacao do meio 2 (seria o seu k2)
p=sqrt(Beta.^2-(n3*ko)^2);      % Constante de propagacao do meio 3 (seria o seu k3)
f=tan(h*tg)-(p+q)./(h.*(1-p.*q./(h.^2))); % Equacao transcendental para modos TE
figure(1)
plot(Beta/ko,f,'.') 
grid
axis([xi xf -10 10])            % eu escolhi este eixo apenas para visualizar melhor
                                % o ponto por onde passa a raiz. Experimente outro.

% Os comandos a seguir sao uma maneira bem simples (mas nao tao exata) de se encontrar 
% uma raiz

% Inicio da rotina para se encontrar uma raiz (So vale para uma raiz de cada vez)
% Quanto mais proximo o intervalo de variacao xi e xf estiverem da raiz, melhor.

a=min(abs(f));
idx=find(abs(f)==a);
Beta=Beta(idx);
Neff=Beta/ko                    % Aqui esta o indice efetivo que voce tanto procura.

% fim da rotina que encontra a raiz

%axis([n3 n2 -10 10])
Mi_zero=4*pi*1e-7;
Eps_zero=8.854e-12;
co=1/sqrt(Mi_zero*Eps_zero);	        % Velocidade da luz em m/s
co=co/1e-6;			                    % Velocidade da luz em microns/s
omega=2*pi*co/lambda;                   % Frequencia angular

% Agora o que eu quero e plotar o campo. Para isso eu preciso encontrar os valores
% de k1, k2 e k3 (no caso desse programa h, p e q) para o modo em questao. Para isso 
% eu preciso utilizar o valor de Beta que eu encontrei como minha raiz.

q=sqrt(Beta.^2-(n1*ko)^2);
h=sqrt((n2*ko)^2-Beta.^2);
p=sqrt(Beta.^2-(n3*ko)^2);

% Se voce se lembra, os campos foram escritos em funcao de apenas uma constante de integracao.
% No seu caderno escolhemos a constante B, aqui eu escolhi escrever em termos de C.
% O Valor de C foi obtido via integral de normalizacao.

C=2*h*sqrt(omega*Mi_zero/(Beta*(tg+1/q+1/p)*(h^2+q^2)));  

% Agora eu comeco a definir o tamanho de cada uma das tres regioes (camadas) para poder
% plotar o campo

% REGIAO 1 (casca superior)
x1=0:0.01:2;                    % Estou variando de 0 a 3um
y1=C*exp(-q*x1);

% REGIAO 2 (nucleo)

x2=-tg:0.01:0;                  % Estou variando de -tg a 0 (Os eixos aqui estao diferentes do utilizado
                                % em seu caderno
y2=C*(cos(h*x2)-(q/h)*sin(h*x2));

% REGIAO 3 (casca inferior)

x3=-3:0.01:-tg;                % Estou variando de -5um ate -tg
y3=C*(cos(h*tg)+(q/h)*sin(h*tg))*exp(p*(x3+tg));

% Agora e so plotar o campo

figure(2)
plot(x1,y1,x2,y2,x3,y3)


% Agora eu vou plotar as interfaces do guia de onda

eixo=get(gca,'ytick');
line( [-tg -tg], [min(eixo) max(eixo)] )    % Desenha a interface entre os meios 2 e 3
line( [0 0], [min(eixo) max(eixo)] )        % Desenha a interface entre os meios 1 e 2

% Divirtam-se.
% Qualquer duvida, me contactem.
% Ben-Hur Viana Borges
% Sala: 3008
% Fone: 273-9344
% E-mail: benhur@sel.eesc.sc.usp.br