%% Pulsed DFG simulation
% Simulation of DFG for plane wave pulses
% using the Split-step Fourier method.
% This script varies the pump power and runs the DFG simulation
% for each pump power.

clear all;
close all;
clc;
tic;

%% Simulation paramenters
powersteps            = 10;
power                 = linspace(1,150e3,powersteps);
setting.runSimulation = true;
setting.plot          = true;

%% Indepenent parameters --------------------------------------------------
Param.t_FWHM      = 7e-12 * sqrt(2);       %FWHM of (E) pulse in time domain
Param.t_steps     = 2^11;                 %Steps of t-axis
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     = 20e-6;                %Peak power [W] of signal
Param.P.peak2     = 0;                    %Peak power [W] of idler
Param.Temp        = 88;                   %Param.Temperature in degree
Param.chi_eff     = 22e-12;               %m/V

axis.t                 = linspace(-Param.t_0,Param.t_0,Param.t_steps);
axis.z                 = linspace(0,Param.z_0,Param.z_steps );

%% Simulation
if(setting.runSimulation == true)
    Idler.I.depleted       = zeros(powersteps,Param.z_steps,Param.t_steps);
    Signal.I.depleted      = zeros(powersteps,Param.z_steps,Param.t_steps);
    Pump.I.depleted        = zeros(powersteps,Param.z_steps,Param.t_steps);
    
    Idler.I.undepleted     = zeros(powersteps,Param.z_steps,Param.t_steps);
    Signal.I.undepleted    = zeros(powersteps,Param.z_steps,Param.t_steps);
    Pump.I.undepleted      = zeros(powersteps,Param.z_steps,Param.t_steps);
    
    Signal.Energy.depleted = zeros(powersteps,Param.z_steps);
    Idler.Energy.depleted  = zeros(powersteps,Param.z_steps);
    Pump.Energy.depleted   = zeros(powersteps,Param.z_steps);
    
    Signal.Energy.undepleted = zeros(powersteps,Param.z_steps);
    Idler.Energy.undepleted  = zeros(powersteps,Param.z_steps);
    Pump.Energy.undepleted   = zeros(powersteps,Param.z_steps);
    
    Axis.t_ns   = zeros(powersteps,Param.t_steps);
    Axis.z_mm   = zeros(powersteps,Param.z_steps);
    
    % Initialize progressbar
    progressbar;

    for nn=1:powersteps

        [AX1, AX2, A, B, C, D, E, F] = SplitStepDFG(power(nn),Param.t_steps,Param.t_FWHM,true,false);
        
        Axis.t_ns(nn,:) = AX1;
        Axis.z_mm(nn,:) = AX2;

        Signal.I.depleted(nn,:,:)      = A;
        Idler.I.depleted(nn,:,:)       = B;
        Pump.I.depleted(nn,:,:)        = C;
        
        Signal.Energy.depleted(nn,:)   = D;
        Idler.Energy.depleted(nn,:)    = E;
        Pump.Energy.depleted(nn,:)     = F;
        
        % Update progressbar
        progressbar(nn/(1*powersteps));
        
        [AX1, AX2, A, B, C, D, E, F] = SplitStepDFG(power(nn),Param.t_steps,Param.t_FWHM,false,false);
        
        Axis.t_ns(nn,:) = AX1;
        Axis.z_mm(nn,:) = AX2;

        Signal.I.undepleted(nn,:,:)      = A;
        Idler.I.undepleted(nn,:,:)       = B;
        Pump.I.undepleted(nn,:,:)        = C;
        
        Signal.Energy.undepleted(nn,:)   = D;
        Idler.Energy.undepleted(nn,:)    = E;
        Pump.Energy.undepleted(nn,:)     = F;

        % Update progressbar
        progressbar(nn/powersteps);
        disp(['I am ' sprintf('%0.1f',nn/powersteps*100) '% done']);

    end

    disp(['I am done. It took ' sprintf('%0.2f',toc) ' seconds']);

    %save(['Workspace-',datestr(now, 'yyyy-mm-dd--HH-MM-SS')]);
end

%% Plot

if(setting.plot == true)
    
    RUNpowerPlotWaterfall;
    RUNpowerPlotPropagation;
    
end