FDBPM算法matlab程序

有限差分，beam propagation method,matlab

% 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 -------- 
[xxyy] = 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（1num_samples）;           % Terminos Independientes

matrix = zeros（num_samples）;        % Inicializa Matriz

% --------- Generacion de la matriz tridiagonal ---------
for m = 1:1:num_samples
    if （（m>1） && （m        matrix（mm-1） = -G;
        matrix（mm） = 1 + 2*G;
        matrix（mm+1） = -G;
    else
        matrix（11） = 1 + 2*G;
        matrix（12） = -G;
        matrix（num_samplesnum_samples-1） = -G;
        matrix（num_samplesnum_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                d（lc） = G*modo（irlc-1） + （1 - 2*G）*modo（irlc） + G*modo（irlc+1）;
            else
                if （lc == 1）
                    d（1） = eps;
                else
                    d（num_samples） = eps;
               

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件       96287  2007-08-26 19:35  FDBPM3D.jpg
     文件        4452  2007-08-26 19:34  FDBPM3D_free_space.m
     文件        5208  2007-08-26 19:33  FDBPM3D_free_space_slice.m
     文件        4890  2007-08-26 19:26  FDBPM3D_free_space_volume.m
     文件        1320  2014-02-12 14:35  license.txt