% Beam Propagation Method 30 Marzo 2007
% Edgar Guevara Codina 
% Dispositivos Optoelectronicos 
% Propagacion de un pulso gaussiano por el modulador Mach-Zender
close all; clear all;       % cierra ventanas y borra variables
%----------- Declaracion de Variables ------------------------------- 
xa = -100e-6;               % Coordenada Inicial
xb = 100e-6;                % Coordenada Final
num_samples = 1000;         % Numero de muestras (par)
H = xb - xa;                % Tama駉 del dominio espacial
del = H/num_samples;        % Espaciado de las muestras
x = linspace (xa, xb-del, num_samples);     % Dominio espacial
A = 1;                      % Amplitud del pulso
x0 = 0.1e-6;                  % Desplazamiento
W0 = 1.2e-6;                  % Radio de la cintura del pulso
nmax = 3.492;               % Indice de refraccion maximo 3.492
nmin = 3.481;               % Indice de refraccion minimo 3.481
lambda = 1e-6;           % Longitud de onda
k0=(2*pi/lambda);           % Numero de onda
deltaz = 1e-6;              % Intervalo de Propagacion
zmax = 1800;                 % Distancia Maxima de Propagacion (um)
% ---------- Reticula para graficar pulso propagado -------------------
[xx,yy] = meshgrid ([xa:del:xb-del],[1:1:zmax]);
% ---------- Generacion del pulso gaussiano ---------------------------
modo = A*exp (-((x+x0)/W0).^2);  % Pulso Gaussiano
dftmodo = fix(fft(modo));   % Transformada Discreta de Fourier de dicho pulso
%figure(1);
% Grafica el pulso original
%subplot(3,1,1); plot(x,modo,'r.');grid on; xlim([xa xb]);
%title(sprintf('Pulso Gaussiano Original %de^{-(x+%g/%g)^2}',A,x0,W0));
% Grafica la magnitud del pulso original
%subplot(3,1,2); plot(abs(dftmodo),'g.');grid on; title('Magnitud');
% Grafica la fase del pulso original
%subplot(3,1,3); plot(angle(dftmodo),'b.');grid on; title('Fase');
% --------------- Obtencion del indice promedio -----------------------
nbar = (nmax + nmin)/2;
% ---------- Generacion del perfil de la guia de onda -----------------
zz = imread('machzender.bmp','BMP');
zz=+zz;
zz=zz.';
%figure(2);
for j = 1:1:num_samples,
    for k = 1:1:zmax,
        if(zz(k,j)==1)
            n(k,j)=nmin;
        else
            n(k,j)=nmax;
        end
    end
end
%Variacion del indice para un desfase de 180