I'm trying to simulate in Matlab a 2D gaussian beam propagation in free propagation region using multi-step Wide-Angle BPM (Pade(2, 2)) with Crank-Nicholson scheme. According to book Pade polynomials coefficients are given by:
ξ1=34−iKΔz(1−α)2
ξ2=116−iKΔz(1−α)4
χ1=34+iKΔzα2
χ2=116+iKΔzα4
where
α is backward/forward component ratio and
α=0.5.
Roots of the numerator and denominator polynomials can be obtained:
α1=12[ξ1−√ξ21−4ξ2]
α2=12[ξ1+√ξ21−4ξ2]
β1=12[χ1−√χ21−4χ2]
β2=12[χ1+√χ21−4χ2]
The propagation is split into two steps:
um+12=(1+α1P)(1+β1P)um
um+1=(1+α2P)(1+β2P)um+12
which yields finite-difference equation
um+12j+β1K2um+12j−1−2um+12j+um+12j+1Δx2+β1K2(k2−K2)um+12j=umj+α1K2umj−1−2umj+umj+1Δx2+α1K2(k2−K2)umj
where k=k0n is position-dependent wavenumber and K=k0n0 is the reference wavenumber. This equation forms a tridiagonal system, from where its coefficient aj,bj,cj,rj are given by:
aj=β1K2Δx2
bj=1−2β1K2Δx2+β1K2(k2−K2)
cj=β1K2Δx2
and
rj is the right side of equation.
This is tridiagonal system I need to solve. To do this I use Thomas Method Algorithm. This method says:
1. Set β=b1,u1=r1β
2. Evaluate for j=2 to n:
γj=cj−1β
β=bj−ajγj
uj=rj−ajuj−1β
- Find for j=1 to n-1
k=n−j
uk=uk−γk+1uk+1
Now I have problem, because I can't really figure out how to correctly implement Thomas algorithm. Since in this situation aj,bj,cj=const only rj changes but I don't understand 1st point in this method. This is my Matlab code:
function u_m_total = WideAnglePade2(k0, step_z, step_x, lateral_x, ...
input_field, z_steps, n)
%% COEFFICIENTS
alpha = 0.5;
xi_1 = 0.75 - (1j * k0 * step_z * (1 - alpha))/2;
xi_2 = 1/16 - (1j * k0 * step_z * (1 - alpha))/4;
chi_1 = 0.75 + (1j * k0 * step_z * alpha)/2;
chi_2 = 1/16 + (1j * k0 * step_z * alpha)/4;
alfa1 = 0.5 * (xi_1 - sqrt(xi_1^2 - 4 * xi_2));
alfa2 = 0.5 * (xi_1 + sqrt(xi_1^2 - 4 * xi_2));
beta1 = 0.5 * (chi_1 - sqrt(chi_1^2 - 4 * chi_2));
beta2 = 0.5 * (chi_1 + sqrt(chi_1^2 - 4 * chi_2));
% caching coefficients to speed up computation
coeff_a1 = beta1 / (k0^2 * step_x^2);
coeff_b1 = 1 - 2 * beta1 / (k0^2 * step_x^2) + beta1 / k0^2 * ((k0 * n)^2 - k0^2);
coeff_c1 = beta1 / (k0^2 * step_x^2);
coeff_a2 = beta2 / (k0^2 * step_x^2);
coeff_b2 = 1 - 2 * beta2 / (k0^2 * step_x^2) + beta2 / k0^2 * ((k0 * n)^2 - k0^2);
coeff_c2 = beta2 / (k0^2 * step_x^2);
beta_1 = coeff_b1;
beta_2 = coeff_b2;
% creating vectors holding values for m_half_z and m_z
% I have problem in here, I don't know what to assign and where
u_m_half = zeros(z_steps, lateral_x);
u_m_total(1, :) = input_field;
u_m_total(1, 1) = input_field(1) / beta2;
for m=1:z_steps-1
disp(m);
%% STEP 1 PROPAGATION ΔZ/2
gamma = coeff_c1 / beta_1;
beta = coeff_b1 - coeff_a1 * gamma;
for j=2:lateral_x-1
r = u_m_total(m, j) + alfa1/k0^2 * (u_m_total(m, j-1) - ...
2*u_m_total(m,j) + u_m_total(m, j+1))/step_x^2 + alfa1/k0^2 *...
((k0 * n)^2 - k0^2) * u_m_total(m, j);
u_m_half(m, j) = (r - coeff_a1 * u_m_total(m, j-1))/beta;
end
for j=1:lateral_x-1
k=lateral_x-j;
u_m_half(m, k) = u_m_half(m, k) - gamma * u_m_half(m, k+1);
end
%% STEP 2 PROPAGATION ΔZ
gamma = coeff_c2 / beta_2;
beta = coeff_b2 - coeff_a2 * gamma;
for j=2:lateral_x-1
r = u_m_half(m, j) + alfa2/k0^2 * (u_m_half(m, j-1) - ...
2*u_m_half(m,j) + u_m_half(m, j+1))/step_x^2 + alfa2/k0^2 ...
* ((k0 * n)^2 - k0^2) * u_m_half(m, j);
u_m_total(m+1, j) = (r - coeff_a2 * u_m_half(m, j-1))/beta;
end
for j=1:lateral_x-1
k=lateral_x-j;
u_m_total(m+1, k) = u_m_total(m+1, k) - gamma * u_m_total(m+1, k+1);
end
end
end
But instead of simulating dispersion and broadening of beam, its amplitude goes towards +Inf.
This post imported from StackExchange Physics at 2017-01-11 10:31 (UTC), posted by SE-user Colonder