%% Scalar Difference Frequency Generation
% Simulates DFG for monochromatic, CW, plane waves.
%
% Input:
%   TolAbs  Absolute tolerence used in RungeKutta45
%   TolRel  Relative tolerence used in RungeKutta45
%   Pp1     Signal seed power
%   Pp2     Idler seed power
%   Pp3     Pump power
%   zsteps  Number of steps that shall be returned from the RungeKutta45
%   zlength Length of crystal
%   QPM     Order of QPM structure
%   QPMtype Type of QPM simulation; 'eff' (effeective) or 'domain'
%
% Output:
%   output0     z axis
%   output1     Signal power;
%   output2     Idler power;
%   output3     Pump power
%   output4     Total power
%   output5     Signal photon number
%   output6     Idler photon number
%   output7     Pump photon number
%   output8     Total photon number


function [output0, output1, output2, output3, output4, output5, output6,...
                    output7, output8] = ScalarDFG(TolAbs, TolRel,Pp1,...
                    Pp2, Pp3, zsteps,zlength, QPM, QPMtype)

%% Independent Parameters
Idler.Lambda  = 1.570e-6;           %wavelength of idler in [m]
Pump.Lambda   = 1.064e-6;           %wavelength of pump in [m]
Param.Length  = zlength;            %Length of nonlinear crystal in [m]
Param.Area    = pi*(40e-6)^2;       %Area of the nonlinear crystal in [m^2]
Param.Temp    = 88;                 %Temperature in degree
Param.QPM     = QPM;                %QPM: 0=NON, 1=1st order, 2=2th order,...
Param.QPMtype = QPMtype;            %Type of QPM: domain og eff.
Param.zsteps  = zsteps;             %Number of z-steps


%% Boundary conditions
Signal.PeakPower    = Pp1;          %power in watt - signal
Idler.PeakPower     = Pp2;          %power in watt - idler
Pump.PeakPower      = Pp3;          %power in watt - pump


%% 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;       %eff. nonlinear suceptibility in [m/V]

%% Dependent Parameters
Idler.Omega     = 2*pi*const.c / Idler.Lambda;
Pump.Omega      = 2*pi*const.c / Pump.Lambda;
Signal.Omega    = Pump.Omega - Idler.Omega;
Signal.Lambda   = 2*pi*const.c / Signal.Omega;

Signal.Omega / Pump.Omega;
Pump.Lambda / Signal.Lambda;

%% Refractive index from Sellmeier equations
[n , kVector] = sellmeier([Signal.Lambda Idler.Lambda Pump.Lambda],Param.Temp);
% k vector
Signal.kVector  = kVector(1);
Idler.kVector   = kVector(2);
Pump.kVector    = kVector(3);
% refractive index
Signal.n        = n(1);
Idler.n         = n(2);
Pump.n          = n(3);

deltak  = Pump.kVector - Signal.kVector - Idler.kVector;
lcoh    = 2 * pi / deltak;
if( strcmp(Param.QPMtype , 'eff') )
   const.chi  = 2/(pi*Param.QPM)*const.chi;
   lcoh       = 1e10;
   deltak     = 0;
end


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

%calculate amplitudes
Signal.A    = sqrt(Signal.PeakPower / (const.hbar*Signal.Omega)/Param.Area);
Idler.A     = sqrt(Idler.PeakPower  / (const.hbar*Idler.Omega)/Param.Area);
Pump.A      = 1i*sqrt(Pump.PeakPower   / (const.hbar*Pump.Omega)/Param.Area);

%% Solve set of ODEs
    % Sets the global and local maximum fault tolerance. 
options = odeset('AbsTol',TolAbs,'RelTol',TolRel);
    % Solve ODE: ode45(solve for , interval , init.values , options).

[z, A]  = ode45(@(z,A) ode45equation(z,A,g,deltak,Param.QPM,lcoh), ...
    linspace(0,Param.Length,Param.zsteps), [Signal.A, Idler.A, Pump.A], options);

%% Find number of photons instead of amplitudes
Signal.N    = A(:,1).*conj(A(:,1)) ;
Idler.N     = A(:,2).*conj(A(:,2)) ;
Pump.N      = A(:,3).*conj(A(:,3)) ;
Total.N     = Signal.N + Idler.N + Pump.N;

%% Find electric fields
Signal.E    = A(:,1) * sqrt( 2 * eta * const.hbar * Signal.Omega) .* exp(1i .* Signal.kVector .* z);
Idler.E     = A(:,1) * sqrt( 2 * eta * const.hbar * Idler.Omega ) .* exp(1i .* Idler.kVector .* z);
Pump.E     = A(:,1) * sqrt( 2 * eta * const.hbar * Pump.Omega ) .* exp(1i .* Pump.kVector .* z);

%% Find powers instead of amplitudes
Signal.P    = Signal.N .* Param.Area .* const.hbar .* Signal.Omega;
Idler.P     = Idler.N .* Param.Area .* const.hbar .* Idler.Omega;
Pump.P      = Pump.N .* Param.Area .* const.hbar .* Pump.Omega;
Total.P     = Signal.P + Idler.P + Pump.P;

%% Z axis
Axis.zmili  = z*1e3;

output0 = Axis.zmili;
output1 = Signal.P;
output2 = Idler.P;
output3 = Pump.P;
output4 = Total.P;
output5 = Signal.N;
output6 = Idler.N;
output7 = Pump.N;
output8 = Total.N;

end