%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function calculates the Second Harmonic Generation
% of a Gaussian of frequencies with a given phase function.
%
%    phi - an input COLUMN vector, containing the phase function.
%    SHG - the output of the calculation; scalar.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [SHG] = SHG(phi);

%constant for consistency with the Fortran Calculation...
c_fortran = 153.7687;

% Generate the Gaussian and the phase function consistently.
Np = length(phi(:,1));
Nv = 4000;
v = linspace(-300,300,Nv); %Linearly Spaced Vector
G = 40;
Ain = exp(-(v/G).^2); %The Gaussian of Frequencies

%Distribute the phase function according to the desired resolution
step = round((2600-1400)/Np);
step = step + (step==1);
phase = zeros(1,Nv);
k = 1;
for j = 1400:step:2600-step+1,
    phase([j:j+step-1]) = phi(k);
    if (k < Np)
        k = k+1;
    else
        k = Np;
    end
end
    % *** The Core: SHG *** %

%Fourier Transform with phase shift (phi) on the Gaussian *Ain*
E_t = fftshift(ifft(fftshift(exp(i*phase).*Ain)));

%plot(abs(E_t));
%Integrate the result to yield the SHG
SHG = sum(abs(E_t).^4)/c_fortran;


%%%%%%%%%%%%%%%%%%%%%%%%%%% E O F %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%