%% 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

%% 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