% Beam Propagation Method
% Originally written by Edgar Guevara Codina 
% Adapted by Peter Bermel
% Propagates a Gaussian pulse through a Y-splitter
close all; clear all;       % clear variables
%----------- Variables ------------------------------- 
xa = -100e-6;               % Initial X-Coordinate (in m)
xb = 100e-6;                % Final X-Coordinate (in m)
num_samples = 1000;         % Number of samples
H = xb - xa;                % Spatial domain
del = H/num_samples;        % X-spacing
x = linspace (xa, xb-del, num_samples);     % Linear Space for X
A = 1;                      % Pulse amplitude
x0 = 0.1e-6;                % Displacement
W0 = 0.8e-6;                % Spatial width of pulse
nmax = 3.492;               % Max refractive index (SiGe)
nmin = 3.481;               % Min refractive index (Si)
lambda = 0.633e-6;          % Wavelength (in m)
k0=(2*pi/lambda);           % Wavevector
deltaz = 1e-6;              % Propagation interval along z
% ---------- Generate Gaussian pulse ---------------------------------
mode = A*exp (-((x+x0)/W0).^2);  % Gaussian pulse
dftmode = fix(fft(mode));   % DFT of Gaussian pulse
%figure(1);
%subplot(3,1,1); plot(x,mode,'r.');grid on; xlim([xa xb]);
%title(sprintf('Original Gaussian pulse %de^{-(x+%g/%g)^2}',A,x0,W0));
%subplot(3,1,2); plot(abs(dftmode),'g.');grid on; title('Magnitude');
%subplot(3,1,3); plot(angle(dftmode),'b.');grid on; title('phase');
% --------------- Calculate average refractive index ------------------
nbar = (nmax + nmin)/2;
% ---------- Generate waveguide profile -------------------------------
zz = imread('ybranch.bmp','BMP');       %Upload image with the profile
%zz = imread('machzender.bmp','BMP');
zz=+zz;                                 %Convert to real numbers
zz=zz.';                                %Transpose
zzdims=size(zz);
zmax=zzdims(1);
% ---------- Graphic display of pulse --------------------------------
[xx,yy] = meshgrid ([xa:del:xb-del],[1:1:zmax]);

for j = 1:1:num_samples,                %Convert to refractive indices
    for k = 1:1:zmax,
        if(zz(k,j)==1)
            n(k,j)=nmin;
        else
            n(k,j)=nmax;
        end
    end
end
%surf(xx,yy,n);             %Display waveguide profile
%view (0,90);
%shading interp;
%xlim([xa xb]);
%plot(x,n,'b.'); title('Refractive index profile'); xlim([xa xb]);
% ---------------- Calculate Phase 1 (normalized) -------------------
kx1 = linspace(0,num_samples/2 + 1,num_samples/2);
kx2 = linspace(-num_samples/2,-1,num_samples/2);
kx = (2*pi/H)*[kx1 kx2];
%figure(3);
%plot(kx,'r.'); title('Wavevector Kx');
phase1 = exp((i*deltaz*kx.^2)./(nbar*k0 + sqrt(max(0,nbar^2*k0*2 - kx.^2))));
% --------------- Calculate signal attenuation -------------------------
for j = 1:1:num_samples,
    if (j >= 370) && (j <= 630)
        od(j) = 0;
    else
        od(j) = 7500;
    end
end
%figure(4);
%plot(x,od,'g.'); title('Signal attenuation'); xlim([xa xb]);

% --------------- Main program loop ------------------------
%figure(5);
%hold on;
for k = 1:1:zmax,
        %Calculate the phase shift contributed by the lens effect
        phase2 = exp(-(od + i*(n(k,:) - nbar)*k0)*deltaz);
        % Convolution via FFT
        mode = ifft((fft(mode).*phase1)).*phase2;
        zz(k,:) = abs(mode);
 end
%surf(xx,yy,zz);
%contour(xx,yy,zz);
contourf(zz);
colormap hot;
colorbar;
%shading flat;
%grid on;
%view (0,90);
%xlim([xa/4 xb/4]);
xlim([350 650]);
title('Y-Splitter: envelope amplitude');
%figure(6);
%contour(xx,yy,zz,30);
%xlim([xa/4 xb/4]);