fdbpm3d_free_space_volume.m

计算有限差分与传播光束法结合的程序

% 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 un volumen del haz propagado tal y como se veria en 3-D

close all;                          % Cierra Ventanas
clear all;                          % Limpia Variables
clc;                                % Limpia Pantalla
tic                                 % Comienza Timer

% ---------- Declaracion de Variables -----------------
x1 = -50e-6;                                % Coordenada Inicial
x2 = 50e-6;                                 % Coordenada Final
num_samples = 64;                           % 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 el pulso -------- 
[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

% -------- Generamos la reticula para volumen --------
[xx,yy,zz] = meshgrid (x1:dx:x2-dx,x1:dx:x2-dx,1:1:1000);
volumen = zeros(size(zz));
clear xx yy zz;

% ---------- 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

% ------------- Liberacion de Memoria --------------
clear W0 dz k0 lambda x;
% --------- 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) = G*modo(num_samples,lc) + (1 - 2*G)*modo(ir,lc) + G*modo(ir+1,lc);
                    d(1) = eps;
                else
                    %d(num_samples) = G*modo(ir-1,lc) + (1 - 2*G)*modo(ir,lc) + G*modo(1,lc);
                    d(num_samples) = eps;
                end
            end
        end
        modo(lc,:) = (matrix\d.')';              % Resuelve la l-esima fila
    end
% --------- Captura de las "rebanadas" del volumen ---------
    volumen(:,:,m) = abs(modo);
end

% ------------- Liberacion de Memoria --------------
clear matrix d ir lc m;

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

% ------------- Liberacion de Memoria --------------
clear B G dx modo num_samples x1 x2;

% ---------- Graficamos el rayo propagado ----------
hpatch = patch(isosurface(xx,yy,zz,volumen,0.09));
isonormals(xx,yy,zz,volumen,hpatch);
set(hpatch,'FaceColor','red','EdgeColor','none');
set(gca,'Color','black','XColor','white','YColor','white','ZColor','white');
set(gcf,'Renderer','zbuffer'); lighting phong;
daspect([1 1 5000000]);
axis tight;
box off;
camlight(hpatch,129)
view([-0.7312   -0.0000    3.6194   -1.4441;0.3353    0.8710    1.9053   -1.5558;0.5944   -0.4913    3.3775   32.1870;    0         0         0    1.0000]);
toc                             %Termina Timer