function [dxdy,Attenuation_Signal,Attenuation_Pump] = SRS_power_penalty(paramAttenuationNorm, paramFiberType, paramMinMax, paramReach, paramSignalPower, paramPumpPower, paramSignalWavelength, paramPumpWavelength, paramApproximatedFlag)

% -------------------------------------------------------------------------
%                                   PARAMETERS
% -------------------------------------------------------------------------

% paramAttenuationNorm - fiber attenuation normalisation factor (0 ...
% INF), which will be used to multiply the G652 fibre attenuation for the
% given wavelength to scale it up/down ... 

% paramFiberType - fibre type for G652 cable - A ... D types allowed

% paramMinMax - flag indicating whether max or min values for the fibre
% attenuation are to be used

% paramReach - the length of the fiber in km

% paramSignalPower - video overlay signal power (dBm)

% paramPumpPower - digital signal power (dBm)

% paramSignalWavelength - video overlay wavelength (nm)

% paramApproximatedFlag - digital signal wavelength (nm)

% paramApproximatedFlag - flag indicating whether simplified (true) or
% detailed model are to be used ...

% -------------------------------------------------------------------------
%                                PARAMETER TESTS
% -------------------------------------------------------------------------

% check if all the parameters apart from paramReach have only one value
if (max(size(paramAttenuationNorm)) ~= 1 | max(size(paramSignalPower)) ~= 1 | max(size(paramPumpPower)) ~= 1 | max(size(paramSignalWavelength)) ~= 1 | max(size(paramPumpWavelength)) ~= 1 | max(size(paramApproximatedFlag)) ~= 1)
    % display a warning notice
    disp('All parameters apart from paramReach must be single values !!!');
    return;
end

if (max(size(paramFiberType)) ~= 1 | (strcmpi(paramFiberType,'A') ~= 1 & strcmpi(paramFiberType,'B') ~= 1 & strcmpi(paramFiberType,'C') ~= 1 & strcmpi(paramFiberType,'D') ~= 1))
    disp('Unknown fibre type - A ...D types allowed !!!');
    return;
end;

if (max(size(paramMinMax)) ~= 3 | (strcmpi(paramMinMax,'MAX') ~= 1 & strcmpi(paramMinMax,'MIN') ~= 1))
    disp('MIN/MAX flags only allowed for paramMinMax !!!');
    return;
end;    

% parameters were tested, can continue with the calculations

% -------------------------------------------------------------------------
%                           INITIAL CALCULATIONS
% -------------------------------------------------------------------------

% recalculate the power levels from dBm into mW which is necessary for
% further manipulation
SignalPower = 10.^(paramSignalPower/10);
PumpPower = 10.^(paramPumpPower/10);

% calculate the fibre attenuation  at the signal and pump wavelengths in
% accordance with the normalization factor and the min/max setting 

[min_val, max_val] = AttenuationG652(paramSignalWavelength,paramFiberType);
Attenuation_Signal = min_val;
if (strcmpi(paramMinMax,'MAX') == 1)
    Attenuation_Signal = max_val;
end; 

[min_val, max_val] = AttenuationG652(paramPumpWavelength,paramFiberType);
Attenuation_Pump = min_val;
if (strcmpi(paramMinMax,'MAX') == 1)
    Attenuation_Pump = max_val;
end;

% normalize the fibre attenuation for signal wavelength to paramAttenuationNorm
% and at the pump wavelength proportionally ...
if (paramAttenuationNorm ~= 0.0)
    Attenuation_Pump = Attenuation_Pump * (paramAttenuationNorm/Attenuation_Signal);
    Attenuation_Signal = paramAttenuationNorm;
end; 

% -------------------------------------------------------------------------
%                           FINAL CALCULATIONS
% -------------------------------------------------------------------------

if (paramApproximatedFlag == true)
    % if the approximated version is to be used, use the simplification as
    % presented in the following presentation (page 19 and 20):
    % http://www.ieee802.org/3/av/public/2007_01/3av_0701_ten_2.pdf    
    dxdy = zeros(2,max(size(paramReach)));
    dxdy(1,:) = paramReach;
    dxdy(2,:) = SignalPower .* length_effective(Attenuation_Signal,paramReach) .* Kappa(paramSignalWavelength,paramPumpWavelength) .* Cr(Frequency_delta(paramSignalWavelength,paramPumpWavelength));
    % completed successfully 
else    
    % unapproximated approach uses the complete model with undepleted signal
    % approximation as posted in the same document on page 19    
    % http://www.ieee802.org/3/av/public/2007_01/3av_0701_ten_2.pdf
    
    % make two calculation passes 
    %     1º - calculate the system with the video overlay
    %     2º - calculate the system with no video overlay
    %     then calculate the abs difference between two systems for digital
    %     signal propagation
    [x,y_with] = ode45(@delta_dB,paramReach,[SignalPower PumpPower],[],Attenuation_Signal,Attenuation_Pump,paramSignalWavelength,paramPumpWavelength); 
    [x,y_without] = ode45(@delta_dB,paramReach,[0 PumpPower],[],Attenuation_Signal,Attenuation_Pump,paramSignalWavelength,paramPumpWavelength);     
    result_abs = y_with(:,2) ./ y_without(:,2);
    % the result is given in mW - must be recalculated into dBm
    result_abs = abs(10*log10(result_abs));
    % prepare the output format ... 
    dxdy = zeros(2,max(size(paramReach)));
    dxdy(1,:) = x;
    dxdy(2,:) = result_abs;
    % completed successfully 
end 


% -------------------------------------------------------------------------
%                          INTERNAL FUNCTIONS
% -------------------------------------------------------------------------

% -------------------------------------------------------------------------
% --------------------- subfunction attenuation ---------------------------
% -------------------------------------------------------------------------
function dx = recalculateAttenuation(paramAttenuationNorm)

% paramAttenuationNorm     - the attentuation of the fibre given in dB/km

% description: recalculate the dB/km of attenuation into 1/km units

dx = abs(paramAttenuationNorm.*(-1).*log(10)/10); 

% -------------------------------------------------------------------------
% --------------------- subfunction length_effective ----------------------
% -------------------------------------------------------------------------
function dx = length_effective(paramAttenuation, paramReach)

% paramAttenuationNorm - the attentuation of the fibre given in dB/km
% paramReach       - the length of the fiber in km

% description: calculate the effective length of the fibre using the 1/km
% attenuation figure and the target reach value

% recalculate the attenuation from dB/km to 1/km, where the 
paramAttenuationNorm_recalc = recalculateAttenuation(paramAttenuation);
% calculate the effective length and return it to the calling method
dx = (1 - 1./exp(paramAttenuationNorm_recalc.*paramReach))./paramAttenuationNorm_recalc;

% -------------------------------------------------------------------------
% -------------------------- subfunction Cr -------------------------------
% -------------------------------------------------------------------------
function dx = Cr(param_delta_f)

% param_delta_f - frequency difference between the signal and pump lamda

% description: % calculate the value of the approximated Cr coefficient 
% (polarization averaged Raman gain coefficient) for the given system
% use the approximated equation for Cr based on the following material:
% http://www.ieee802.org/3/av/public/2007_01/3av_0701_ten_2.pdf

dx = 0.000000230181.*param_delta_f^3-0.00000143492.*param_delta_f^2+0.00000747534.*param_delta_f+0.0000913764;

% -------------------------------------------------------------------------
% -------------------------- subfunction Kappa ----------------------------
% -------------------------------------------------------------------------
function dx = Kappa(paramSignalWavelength,paramPumpWavelength)

% paramSignalWavelength    - signal wavelength in nm
% paramPumpWavelength      - pump wavelength in nm

% description: calculate the Kappa parameter using paramSignalWavelength and
% paramPumpWavelength parameters (no units)

dx = 10*log10(exp(1)).*paramSignalWavelength./paramPumpWavelength;

% -------------------------------------------------------------------------
% ---------------------- subfunction Frequency_delta ----------------------
% -------------------------------------------------------------------------
function dx = Frequency_delta (paramSignalWavelength,paramPumpWavelength)

% description: calculate the frequency difference between the signal and pump 
% wavelengths (expressed in THz)

dx = abs(3./paramSignalWavelength-3./paramPumpWavelength)*10^5;

% -------------------------------------------------------------------------
% -------------------------- subfunction delta_dB -------------------------
% -------------------------------------------------------------------------
function dxdy = delta_dB(param_x,param_y,Attenuation_Signal,Attenuation_Pump,paramSignalWavelength,paramPumpWavelength)

% param_y(1) = Ps 
% param_y(2) = Pp

% precalculate some of the values
delta_f = Frequency_delta(paramSignalWavelength,paramPumpWavelength);
Attenuation_Signal_recalculated = recalculateAttenuation(Attenuation_Signal);
Attenuation_Pump_recalculated = recalculateAttenuation(Attenuation_Pump);
% initialize the output vector 
dxdy = zeros(size((param_y),1),1);
% calculate the elements of the Ps and Pp
dxdy(1) = param_y(1) .* param_y(2) .* Cr(delta_f) - abs(Attenuation_Signal_recalculated .* param_y(1));
dxdy(2) = - paramSignalWavelength ./ paramPumpWavelength .* Cr(delta_f) .* param_y(1) .* param_y(2) - abs(Attenuation_Pump_recalculated .* param_y(2));