Simulation of SFG and DFG for a Gaussian beam profile

General: omega_3 > omega_2 > omega_1

SFG: omega_1 + omega_2 = omega_3 Fundamental + Idler = signal

DFG: omega_3 - omega_2 = omega_1 Fundamental - Idler = signal

Initialize
Settings
check for workspace
Simulation
Create a complete Gaussian basis
Boundary conditions
Using effective 1st order QPM structure
Propagate electric field
Calculating intensity to plot
Save
Plots figures

Initialize

clear all;
close all;
clc;

Settings

setting.Plot         	= true; %either true/false or vector [1 2]
setting.SaveFigure   	= true;
setting.Movie        	= false;
setting.Depletion       = true;

check for workspace

Checks if there is workspace files in the working directory if so, then choose to stop, ignore, or overwrite the workspace file

if(FindNewestWorkspace('Workspace-2','mat',false) == true)
    choice = questdlg('Do you wish to simulate anyway?', ...
    'There is already a workspace file present.', ...
    'No','Yes - do not save','Yes - override','No');
    % Handle response
    switch choice
        case 'No'
            disp('Not simulating - loading existing workspace');
            setting.runSimulation = false;
            setting.SaveWorkspace = false;
            file = FindNewestWorkspace('Workspace-2','mat',true);
            disp(['Loading file: ' file]);
            load( file );
        case 'Yes - do not save'
            disp('Simulating even though there is an existing workspace');
            setting.runSimulation = true;
            setting.SaveWorkspace = false;
        case 'Yes - override'
            disp('Simulating and deleting old workspace');
            delete('Workspace-2*.mat')
            setting.runSimulation = true;
            setting.SaveWorkspace = true;
    end
else
    setting.runSimulation = true;
    setting.SaveWorkspace = false;
end

Simulation

ticAll = tic;
if(setting.runSimulation == true)

    progressbar

Create a complete Gaussian basis

    % Define parameters
    Param.rsteps        = 5003;
    Param.zsteps        = 1000;
    Param.Lmax          = 10;
    Param.z0            = 20e-3;

    % Independent Parameters
    Pump.lambda         = 1.064e-6;         %wavelength of pump in [m]
    Idler.lambda        = 1.570e-6;         %wavelength of idler in [m]
    omega0pumpmicro     = 1000e-6;          %Beam wais of pump [m]
    Temp                = 88;               %Temperature in degree
    Lqpm                = 30.49e-6;         %Periode of QPM structure

Boundary conditions

    Pump.I0             = 1500e9;           %power in watt
    Signal.I0           = 1e6;              %power in watt
    Idler.I0            = 1e6;              %power in watt

    % Constants
    Const.c             = 299792458;        %Speed of light in [m/s]
    Const.hbar          = 1.05457173e-34;   %Const.hbar in [m2*kg/s]
    Const.epsilon0      = 8.854187817e-12;  %in [F/m]
    Const.mu0           = pi*4e-7;          %in [V*s/A*m]
    Const.chi           = 22e-12;           %effective nonlinear suceptibility in [m/V]

    % Dependent Parameters
    Pump.omega          = 2*pi*Const.c/Pump.lambda;
    Idler.omega         = 2*pi*Const.c/Idler.lambda;
    Signal.omega        = Pump.omega - Idler.omega;
    Signal.lambda       = 2*pi*Const.c/Signal.omega;
    omega0              = [1 1/sqrt(2) 1/(1+sqrt(2))]...
                        .* omega0pumpmicro ./ sqrt(2); %beam waist in [m]

    % Refractive index from Sellmeier equations
    [temp.n , temp.kVector ] = sellmeier([Pump.lambda, Idler.lambda, Signal.lambda],Temp);
    Pump.kVector        = temp.kVector(1);
    Idler.kVector       = temp.kVector(2);
    Signal.kVector      = temp.kVector(3);
    Pump.n              = temp.n(1);
    Idler.n             = temp.n(2);
    Signal.n            = temp.n(3);

    deltak              = Pump.kVector-Idler.kVector-Signal.kVector;
    lcoh                = 2 * pi / deltak;
    rayleigh            = pi.*[omega0(1)^2 / Pump.lambda , omega0(2)^2 / Signal.lambda];

Using effective 1st order QPM structure

    deltak = 0;
    Const.chi = 2/pi * Const.chi;

    %Nonlinear parameters
    Pump.eta    = sqrt(Const.mu0 / Const.epsilon0) / Pump.n;
    Idler.eta   = sqrt(Const.mu0 / Const.epsilon0) / Idler.n;
    Signal.eta  = sqrt(Const.mu0 / Const.epsilon0) / Signal.n;
    eta         = nthroot(Pump.eta * Idler.eta * Signal.eta, 3);
    g           = Const.epsilon0 * Const.chi .* sqrt(1/2 * Const.hbar * eta^3 * Pump.omega * Idler.omega * Signal.omega);

    %Find the radial vector and the propagation vector
    z       = linspace(-Param.z0,Param.z0,Param.zsteps);
    zmili   = z.*1e3;
    r       = linspace(0,10*max(omega0),Param.rsteps);
    deltaz  = z(2)-z(1);
    deltar  = r(2)-r(1);

    disp(['Time spend defining: ' sprintf('%0.2f',toc(ticAll))]);

Time spend defining: 0.12

Propagate electric field

    ticPropagating  = tic;

    % Initialize matrices
    Pump.E          = zeros(Param.rsteps    ,Param.zsteps);
    Pump.C          = zeros(Param.Lmax      ,Param.zsteps);
    Idler.E         = zeros(Param.rsteps    ,Param.zsteps);
    Idler.C         = zeros(Param.Lmax      ,Param.zsteps);
    Signal.E        = zeros(Param.rsteps    ,Param.zsteps);
    Signal.C        = zeros(Param.Lmax      ,Param.zsteps);

    Pump.E(:,1)     = sqrt(2*Pump.I0 /(Const.c*Const.epsilon0*Pump.n))...
                .* GaussianMode(1,r, z(1),Pump.lambda,omega0(1),Pump.n) ./ max(abs(GaussianMode(1,r, z(1),Pump.lambda,omega0(1),Pump.n)));
    Idler.E(:,1)    = sqrt(2*Idler.I0       /(Const.c*Const.epsilon0*Idler.n))...
                .* GaussianMode(1,r, z(1),Idler.lambda,omega0(2),Idler.n) ./ max(abs(GaussianMode(1,r, z(1),Idler.lambda,omega0(2),Idler.n)));
    Signal.E(:,1)   = 1i.*sqrt(2*Signal.I0  /(Const.c*Const.epsilon0*Signal.n))...
                .* GaussianMode(1,r, z(1),Signal.lambda,omega0(3),Signal.n) ./ max(abs(GaussianMode(1,r, z(1),Signal.lambda,omega0(3),Signal.n)));

    PowerChange = zeros(1,Param.zsteps-1);
    for countk=1:Param.zsteps-1

        Pump.C(:,countk)     = GaussianBasis(Param.Lmax,r,z(countk),Pump.lambda,omega0(1),Pump.n,deltar)    * Pump.E(:,countk);
        Idler.C(:,countk)    = GaussianBasis(Param.Lmax,r,z(countk),Idler.lambda,omega0(2),Idler.n,deltar)  * Idler.E(:,countk);
        Signal.C(:,countk)   = GaussianBasis(Param.Lmax,r,z(countk),Signal.lambda,omega0(3),Signal.n,deltar)* Signal.E(:,countk);

        Pump.Eprime          = 0;
        Idler.Eprime         = 0;
        Signal.Eprime        = 0;

        for countn=1:Param.Lmax
            Pump.Eprime      = Pump.Eprime   + Pump.C(countn,countk)   * transpose(GaussianMode(countn,r, z(countk+1),Pump.lambda,omega0(1),Pump.n));
            Idler.Eprime     = Idler.Eprime  + Idler.C(countn,countk)  * transpose(GaussianMode(countn,r, z(countk+1),Idler.lambda,omega0(2),Idler.n));
            Signal.Eprime    = Signal.Eprime + Signal.C(countn,countk) * transpose(GaussianMode(countn,r, z(countk+1),Signal.lambda,omega0(3),Signal.n));
        end

        Pump.A               = Pump.Eprime   ./ sqrt(2*Pump.eta*Const.hbar*Pump.omega)  ;%  .* exp(-1i*Pump.kVector*z(countk+1));
        Idler.A              = Idler.Eprime  ./ sqrt(2*Idler.eta*Const.hbar*Idler.omega) ;%  .* exp(-1i*Idler.kVector*z(countk+1));
        Signal.A             = Signal.Eprime ./ sqrt(2*Signal.eta*Const.hbar*Signal.omega);%  .* exp(-1i*Signal.kVector*z(countk+1));

        if( setting.Depletion == true )
            Pump.deltaA      = 1i.*g.*deltaz.*(Signal.A.* Idler.A).* exp(-1i * deltak*z(countk+1));
        else
            Pump.deltaA      = 0;
        end
        Idler.deltaA         = 1i.*g.*deltaz .* ( Pump.A  .* conj( Signal.A) )    .* exp(1i * deltak*z(countk+1));
        Signal.deltaA        = 1i.*g.*deltaz .* ( Pump.A  .* conj( Idler.A) )     .* exp(1i * deltak*z(countk+1));

        Pump.A               = Pump.A   + Pump.deltaA;
        Idler.A              = Idler.A  + Idler.deltaA;
        Signal.A             = Signal.A + Signal.deltaA;

        PowerChange(countk)  = Const.hbar .* Pump.omega .* 2*pi*deltar*sum(transpose(r).*(conj(Pump.deltaA).*Pump.deltaA)) - Const.hbar .* Signal.omega .* 2*pi*deltar*sum(transpose(r).*(conj(Signal.deltaA).*Signal.deltaA));

        Pump.E(:,countk+1)   = sqrt(2*Pump.eta*Const.hbar*Pump.omega)  .* Pump.A   ;% .* exp(1i*Pump.kVector*z(countk+1));
        Idler.E(:,countk+1)  = sqrt(2*Idler.eta*Const.hbar*Idler.omega) .* Idler.A         ;% .* exp(1i*Signal.kVector*z(countk+1));
        Signal.E(:,countk+1) = sqrt(2*Signal.eta*Const.hbar*Signal.omega).* Signal.A        ;% .* exp(1i*Signal.kVector*z(countk+1));

        progressbar(countk/(Param.zsteps-1));
    end

    Pump.C(:,Param.zsteps)   = GaussianBasis(Param.Lmax,r,z(Param.zsteps),Pump.lambda,omega0(1),Pump.n,deltar)      * Pump.E(:,Param.zsteps);
    Idler.C(:,Param.zsteps)  = GaussianBasis(Param.Lmax,r,z(Param.zsteps),Idler.lambda,omega0(2),Idler.n,deltar)    * Idler.E(:,Param.zsteps);
    Signal.C(:,Param.zsteps) = GaussianBasis(Param.Lmax,r,z(Param.zsteps),Signal.lambda,omega0(2),Signal.n,deltar)  * Signal.E(:,Param.zsteps);
    clearvars Pump.A Idler.A Signal.A Pump.Eprime Idler.Eprime Signal.Eprime;

Calculating intensity to plot

    rplot               = floor(Param.rsteps/4);
    raxis               = [fliplr(-r(1:rplot)) r(1:rplot)] .* 1e3;
    Pump.Eplot          = [fliplr( transpose(Pump.E(1:rplot,:)) )    transpose(Pump.E(1:rplot,:))];
    Idler.Eplot         = [fliplr( transpose(Idler.E(1:rplot,:)) )   transpose(Idler.E(1:rplot,:))];
    Signal.Eplot        = [fliplr( transpose(Signal.E(1:rplot,:)) )  transpose(Signal.E(1:rplot,:))];
    Pump.Iplot          = Const.c*Const.epsilon0/2 * Pump.n    * abs(Pump.Eplot).^2;
    Idler.Iplot         = Const.c*Const.epsilon0/2 * Idler.n   * abs(Idler.Eplot).^2;
    Signal.Iplot        = Const.c*Const.epsilon0/2 * Signal.n  * abs(Signal.Eplot).^2;

    % Calculate power
    Pump.Power          = zeros(Param.zsteps,1);
    Idler.Power         = zeros(Param.zsteps,1);
    Signal.Power        = zeros(Param.zsteps,1);
    Power               = zeros(Param.zsteps,1);
    for countn=1:Param.zsteps
        Pump.Power(countn)        = Const.c*Const.epsilon0/2 * Pump.n       * 2*pi*deltar*sum(transpose(r).* ( Pump.E(:,countn)   .* conj(Pump.E(:,countn)) ) );
        Idler.Power(countn)       = Const.c*Const.epsilon0/2 * Idler.n      * 2*pi*deltar*sum(transpose(r).* ( Idler.E(:,countn)  .* conj(Idler.E(:,countn))      ) );
        Signal.Power(countn)      = Const.c*Const.epsilon0/2 * Signal.n     * 2*pi*deltar*sum(transpose(r).* ( Signal.E(:,countn) .* conj(Signal.E(:,countn))      ) );
        Power(countn)             = Pump.Power(countn)+Idler.Power(countn)+Signal.Power(countn);
    end

    Pump.ModePower                = zeros(Param.Lmax,Param.zsteps);
    Idler.ModePower               = zeros(Param.Lmax,Param.zsteps);
    Signal.ModePower              = zeros(Param.Lmax,Param.zsteps);
    EaendMode                     = zeros(Param.rsteps,Param.Lmax);

    for countn=1:Param.Lmax
        EaendMode(:,countn)       = Pump.C(countn,end) * transpose(GaussianMode(countn,r, z(end),Pump.lambda,omega0(1),Pump.n));
        Pump.ModePower(countn,:)  = sum( abs(Pump.C(1:countn,:)).^2 ,1);
        Idler.ModePower(countn,:) = sum( abs(Idler.C(1:countn,:)).^2 ,1 );
        Signal.ModePower(countn,:)= sum( abs(Signal.C(1:countn,:)).^2 ,1 );
    end
    IaendMode = Const.c*Const.epsilon0*Pump.n/2 * abs(EaendMode).^2;

    % Calculate FWHM of both fields
    Pump.FWHM                     = zeros(Param.zsteps,3);
    Idler.FWHM                    = zeros(Param.zsteps,3);
    Signal.FWHM                   = zeros(Param.zsteps,3);
    for countn=1:Param.zsteps
%        [Pump.FWHM(countn,1), Pump.FWHM(countn,2), Pump.FWHM(countn,3)]    = fwhm(raxis,Pump.Iplot(countn,:));
%        [Idler.FWHM(countn,1), Idler.FWHM(countn,2), Idler.FWHM(countn,3)]                      = fwhm(raxis,Idler.Iplot(countn,:));
%        [Signal.FWHM(countn,1), Signal.FWHM(countn,2), Signal.FWHM(countn,3)]                   = fwhm(raxis,Signal.Iplot(countn,:));
    end

    disp(['Max power gain/loss in percentage: ' sprintf('%0.2f',max(abs(Power(:)-Power(1)))/Power(1)*100)]);
    disp(['Time spend propagating: ' sprintf('%0.2f',toc(ticPropagating)) '(' sprintf('%0.2f',toc(ticAll)) ')']);

Max power gain/loss in percentage: 0.19
Time spend propagating: 88.73(88.87)

end

Save

if(setting.SaveWorkspace == true)
   ticSave = tic;
   save(['Workspace-',datestr(now, 'yyyy-mm-dd-HH-MM-SS')], '-regexp', '^(?!(setting|Ea|Eb|)$).');
   disp(['Time spend saving: ' sprintf('%0.2f',toc(ticSave)) '(' sprintf('%0.2f',toc(ticAll)) ')']);
end

Plots figures

GaussianBasisPlot

save(['Gaussiandata-',datestr(now, 'yyyy-mm-dd-HH-MM-SS-'),int2str(omega0pumpmicro*1e6),'mumI',int2str(Pump.I0*1e-9),'GW']);

disp(['Time spend TOTAL: ' sprintf('%0.2f',toc(ticAll))]);

Time spend making figures: 79.17(168.10)
Time spend making video: 0.00(168.10)
Time spend TOTAL: 182.33