Contents
Split-step DFG
Simulation of DFG for plane wave pulses using the Split-step Fourier method.
Pump.Omega > Signal.Omega > Idler.Omega Pump.Omega - Idler.Omega => Signal.Omega
function [AX1, AX2, output_1, output_2, output_3, output_4, output_5, output_6] = SplitStepDFG(peakpower,tsteps,tfwhm,depletion,plotting)
Check for optional parameters
if nargin < 5 plotting = true; end if nargin < 4 depletion = true; end if nargin < 3 tfwhm = 7e-9 * sqrt(2); end if nargin < 2 tsteps = 2^9; end if nargin < 1 peakpower = 10e3; end
Setting
setting.Plot_graphs = plotting;
setting.SaveFigure = true;
setting.Movie = false;
setting.Depletion = depletion;%true;
Indepenent parameters
Param.t_FWHM = tfwhm; %FWHM of (E) pulse in time domain Param.t_steps = tsteps; %Steps of t-axis (2^9) Param.t_0 = 11*Param.t_FWHM; %Span of t; -Param.t_0 : Param.t_0 Param.z_steps = 5000; %Steps of z-axis Param.z_0 = 40e-3; %Span of z; 0 : Param.z_0 Param.A_eff = pi*(50e-6)^2; %Area of beam Param.lambda_3 = 1064e-9; %Wavelength [m] of signal Param.lambda_1 = 1600e-9; %Wavelength [m] of idler Param.P.peak1 = 1e-6; %Peak power [W] of signal Param.P.peak2 = 0; %Peak power [W] of idler Param.P.peak3 = peakpower;%20e3; %Peak power [W] of pump Param.Temp = 88; %Param.Temperature in degree
Constants
const.c = 3e8; %Speed of light in [m/s] const.hbar = 1.05457173e-34; %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_eff = 22e-12; %Effective nonlinear suceptibility in [m/V]
Dependent parameters
t_sigma = Param.t_FWHM / (2*sqrt(2*log(2))); t = linspace(-Param.t_0,Param.t_0,Param.t_steps); %Time axis dt = t(2)-t(1); %Time steps f = linspace(0, 1/dt, Param.t_steps); %Frequency axis f_shift = linspace(-1/(2*dt), 1/(2*dt), Param.t_steps); %Centered frequency axis w = f .*2*pi; %Angular frequency axis w_shift = f_shift .*2*pi; %Angular frequency shifted z = linspace(0,Param.z_0,Param.z_steps); %Crystal axis dz = z(2)-z(1); %Crystal steps % Frequency and wavelength, DFG Idler.Omega = 2*pi*const.c / Param.lambda_1; Pump.Omega = 2*pi*const.c / Param.lambda_3; Signal.Omega = Pump.Omega - Idler.Omega; lambda_2 = 2*pi*const.c / Signal.Omega; % Pulse energy % Sellmeier equation [n,kVector] = sellmeier([Param.lambda_1 lambda_2 Param.lambda_3],Param.Temp); deltak = kVector(3) - kVector(2) - kVector(1); lcoh = 2 * pi / deltak; Idler.v = groupvelocity(Param.lambda_1,Param.Temp); Signal.v = groupvelocity(lambda_2,Param.Temp); Pump.v = groupvelocity(Param.lambda_3,Param.Temp); % Use effective 1st order QPM deltak = 0; const.chi_eff = 2/pi * const.chi_eff;
Initial pulse
Idler.E.t = zeros(Param.z_steps,Param.t_steps); %Allocate time matrix Idler.E.f = zeros(Param.z_steps,Param.t_steps); %Allocate frequency matrix Idler.E.t(1,:) = gaussmf(t,[t_sigma t(ceil(end/2))]); %Defining gaussian pulse Idler.E.t(1,:) = Idler.E.t(1,:) * sqrt( 2 * Param.P.peak1 ... / (const.c * const.epsilon0 * n(1) * Param.A_eff) ); Signal.E.t = zeros(Param.z_steps,Param.t_steps); %Allocate time matrix Signal.E.f = zeros(Param.z_steps,Param.t_steps); %Allocate frequency matrix Signal.E.t(1,:) = gaussmf(t,[t_sigma t(ceil(end/2))]); %Defining gaussian pulse Signal.E.t(1,:) = Signal.E.t(1,:) * sqrt( 2 * Param.P.peak2 ... / (const.c * const.epsilon0 * n(2) * Param.A_eff) ); Pump.E.t = zeros(Param.z_steps,Param.t_steps); %Allocate time matrix Pump.E.f = zeros(Param.z_steps,Param.t_steps); %Allocate frequency matrix Pump.E.t(1,:) = gaussmf(t,[t_sigma t(ceil(end/2))]); %Defining gaussian pulse Pump.E.t(1,:) = Pump.E.t(1,:) * sqrt( 2 * Param.P.peak3 ... / (const.c * const.epsilon0 * n(3) * Param.A_eff) );
Propagate pulses
for nn=1:Param.z_steps-1 %Linear part Idler.E.f(nn,:) = fftshift(FourierT(Idler.E.t(nn,:), dt)); Idler.E.f(nn,:) = Idler.E.f(nn,:) .* exp(1i.*w_shift.*dz * ... (1/Idler.v-1/Signal.v)); %Linear dispersion Idler.E.f(nn,:) = fftshift( Idler.E.f(nn,:) ); Idler.E.t(nn+1,:) = IFourierT(Idler.E.f(nn,:), dt); Signal.E.f(nn,:) = fftshift(FourierT(Signal.E.t(nn,:), dt)); Signal.E.f(nn,:) = Signal.E.f(nn,:) .* exp(1i.*w_shift.*dz * ... (1/Signal.v-1/Signal.v)); %Linear dispersion Signal.E.f(nn,:) = fftshift( Signal.E.f(nn,:) ); Signal.E.t(nn+1,:) = IFourierT(Signal.E.f(nn,:), dt); Pump.E.f(nn,:) = fftshift(FourierT(Pump.E.t(nn,:), dt)); Pump.E.f(nn,:) = Pump.E.f(nn,:) .* exp(1i.*w_shift.*dz * ... (1/Pump.v-1/Signal.v)); %Linear dispersion Pump.E.f(nn,:) = fftshift( Pump.E.f(nn,:) ); Pump.E.t(nn+1,:) = IFourierT(Pump.E.f(nn,:), dt); %Non linear part Idler.E.t(nn+1,:) = Idler.E.t(nn+1,:) + dz.*1i.*Idler.Omega ... .* const.chi_eff ./ (2*n(1)*const.c) .* conj( Signal.E.t(nn,:) )... .* Pump.E.t(nn,:) .* exp(1i * deltak*dz); Signal.E.t(nn+1,:) = Signal.E.t(nn+1,:) + dz.*1i.*Signal.Omega ... .* const.chi_eff ./ (2*n(2)*const.c) .* conj( Idler.E.t(nn,:) ) ... .* Pump.E.t(nn,:) .* exp(1i * deltak*dz); if( setting.Depletion == true ) Pump.E.t(nn+1,:) = Pump.E.t(nn+1,:) + dz.*1i.*Pump.Omega ... .* const.chi_eff ./ ( 2*n(3) * const.c) .* Idler.E.t(nn,:) ... .* Signal.E.t(nn,:) .* exp(-1i * deltak*dz); end end
Calculate Power
Idler.I = (const.c * const.epsilon0 * n(1) /2) * abs(Idler.E.t).^2; Signal.I = (const.c * const.epsilon0 * n(2) /2) * abs(Signal.E.t).^2; Pump.I = (const.c * const.epsilon0 * n(3) /2) * abs(Pump.E.t).^2; Idler.P = Idler.I .* Param.A_eff; %(c * epsilon0 * n(1) * area /2) * sum(abs(Idler.E.t).^2,2); Signal.P = Signal.I .* Param.A_eff; %(c * epsilon0 * n(2) * area /2) * sum(abs(Signal.E.t).^2,2); Pump.P = Pump.I .* Param.A_eff; %(c * epsilon0 * n(3) * area /2) * sum(abs(Pump.E.t).^2,2); Total.P = Idler.P + Signal.P + Pump.P; Idler.Energy = dt.*sum(Idler.P,2); Signal.Energy = dt.*sum(Signal.P,2); Pump.Energy = dt.*sum(Pump.P,2); Total.Energy = dt.*sum(Total.P,2);
Define axis
Axis.t_ns = t .* 1e9; Axis.f_gh_shift = f_shift .* 1e-9; Axis.z_mm = z .* 1e3;
Calculate FWHM
Idler.FWHM = zeros(Param.z_steps,3); Signal.FWHM = zeros(Param.z_steps,3); Pump.FWHM = zeros(Param.z_steps,3); for count_n=1:Param.z_steps Idler.FWHM(count_n,:) = fwhm( Axis.t_ns , Idler.I(count_n,:) ); Signal.FWHM(count_n,:) = fwhm( Axis.t_ns , Signal.I(count_n,:) ); Pump.FWHM(count_n,:) = fwhm( Axis.t_ns , Pump.I(count_n,:) ); end
plot figures
if(setting.Plot_graphs == true) SplitStepDFGplot; end
Warning: Ignoring extra legend entries.
Function output
AX1 = Axis.t_ns; AX2 = Axis.z_mm; output_1 = Signal.I(:,:); output_2 = Idler.I(:,:); output_3 = Pump.I(:,:); output_4 = Signal.Energy(:); output_5 = Idler.Energy(:); output_6 = Pump.Energy(:);
end