有限差分光束传输法

% Finite Difference Beam Propagation Method     18 Mayo 2007
% Edgar Guevara Codina 
% Dispositivos Optoelectronicos 
% Genera un pulso gaussiano en 3D y lo propaga 1500 um a lo largo del eje z
% Utiliza Diferencias Finitas para resolver la ecuacion parabolica en 3D
% Mediante el metodo Implicito de Direccion Alternante (ADI)
% Genera una animacion y la guarda en formato AVI

close all;                          % Cierra Ventanas
clear all;                          % Limpia Variables
clc;                                % Limpia Pantalla

% ---------- Declaracion de Variables -----------------
x1 = -50e-6;                                % Coordenada Inicial
x2 = 50e-6;                                 % Coordenada Final
num_samples = 128;                          % Numero de muestras (potencia de 2)
dx = (x2-x1)/num_samples;                   % Espaciado de las muestras en x
dz = 0.25e-6;                               % Incremento en z
x = linspace (x1, x2-dx, num_samples);      % Dominio espacial (simetrico en x y en y)
W0 = 8e-6;                                  % Radio de la cintura del pulso
lambda = 0.8e-6;                            % Longitud de onda
k0 = 2*pi/lambda;                           % Numero de onda

% -------- Generamos la reticula para graficar -------- 
[xx,yy] = meshgrid ([x1:dx:x2-dx],[x1:dx:x2-dx]);

% ------------ Generacion del pulso -------------------
modo = exp (-(xx/W0).^2-(yy/W0).^2);    % Pulso Gaussiano en 3D

% ---------- Constantes para metodo ADI -----------------
B = j/(2*k0);                       % Constante de difusion
G = B*dz/(dx^2);                    % Parametro de ganancia
d = zeros(1,num_samples);           % Terminos Independientes

matrix = zeros(num_samples);        % Inicializa Matriz

% --------- Generacion de la matriz tridiagonal ---------
for m = 1:1:num_samples,
    if ((m>1) && (m<num_samples))
        matrix(m,m-1) = -G;
        matrix(m,m) = 1 + 2*G;
        matrix(m,m+1) = -G;
    else
        matrix(1,1) = 1 + 2*G;
        matrix(1,2) = -G;
        matrix(num_samples,num_samples-1) = -G;
        matrix(num_samples,num_samples) = 1 + 2*G;
    end
end
matrix=sparse(matrix);  %la convierte a matriz escasamente poblada

% -------------- Posiciona la grafica --------------
scrsz = get(0,'ScreenSize');
figure('Position',[1 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2]);

% --------------- Abre Archivo Video ---------------
%aviobj = avifile('FDBPM_3D_free_space.avi','fps',50,'quality',100);

% --------- Ciclo Principal de Propagacion ---------  

% ------------------ Primer paso -------------------
    for ir = 1:1:num_samples,
        for lc = 1:1:num_samples,
            if ((lc>1) && (lc<num_samples))
                d(lc) = G*modo(ir,lc-1) + (1 - 2*G)*modo(ir,lc) + G*modo(ir,lc+1);
            else
                if (lc == 1)
                    d(1) = eps;
                else
                    d(num_samples) = eps;
                end
            end
        end
        modo(:,ir) = matrix\d.';              % Resuelve la i-esima columna
    end

% ------------------ Segundo paso ------------------
for m = 1:1:1500,
    for lc = 1:1:num_samples,
        for ir = 1:1:num_samples,
            if ((ir>1) && (ir<num_samples))
                d(ir) = G*modo(ir-1,lc) + (1 - 2*G)*modo(ir,lc) + G*modo(ir+1,lc);
            else
                if (ir == 1)
                    d(1) = eps;
                else
                    d(num_samples) = eps;
                end
            end
        end
        modo(lc,:) = (matrix\d.')';              % Resuelve la l-esima fila
    end

% -- Graficamos la magnitud de nuestros pulsos propagados --
    subplot(1,2,1);
    contour(xx,yy,abs(modo),16);
    axis equal;
    xlabel('x');    ylabel('y');
    subplot(1,2,2);
    surf(xx,yy,abs(modo));
    shading interp;
    axis([x1 x2 x1 x2 0 1]);
    title(strcat(sprintf('Pulso Propagado en z = %3.0f ',m),' \mum'));
    xlabel('x');    ylabel('y');    zlabel('Amplitud');
    subplot(1,2,1);
    pause(0.005);
    %if (m==250)
        %print -djpeg100 3DFDBPM_0250.jpg       %Exporta la imagen
    %end
% ------- Captura y a馻de un cuadro al video -------
    %frame = getframe(gca,[-80 -40 scrsz(3)/2 scrsz(4)/2+40]);
    %aviobj = addframe(aviobj,frame);
end

% -------------- Cierra Archivo Video --------------
%aviobj = close(aviobj);

 D141