%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% .m script of web demo "Amplified fiber link" 
% Written for GNU Octave
% 
% INSTITUTE OF TELECOMMUNICATIONS  
% University of Stuttgart 
% www.inue-uni-stuttgart.de 
% author: David Rörich
% date: 04.02.2014
%
% EXAMPLES:
% Simple OSNR plot and dispersion plot.
% amplified_fiber_link([1 0 0], 1000, 50, 5000, 0, 0.2, 5, 1550e-9, 80, 50e9, 17, 0.053, -100, 0.3, 0.5, 1, 1, 0, 1, 'both', 0, 'example1.svg');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BEGIN main function
function amplified_fiber_link(...
    plot_en,...                     % plot enable vector [evolution, span_length]
    total_length,...                % in km
    span_length,...                 % in km
    sampling_length,...             % in m
    input_power,...                 % in dBm
    attenuation_coefficient,...     % in dB/km
    noise_figure,...                % in dB
    center_wavelength,...           % in m
    no_wdm_channels,...             % integer number
    channel_spacing,...             % in Hz
    dispersion_coefficient,...      % in ps / (km*nm)
    dispersion_slope,...            % in ps / (km*nm*nm)
    dcf_dispersion_coefficient,...  % in s / (m*m)
    dcf_dispersion_slope,...        % in ps / (km*nm*nm)
    dcf_attenuation_coefficient,... % in dB/km
    amplification,...               % boolean
    dispersion_compensation,...     % boolean  
    webdemo_en,...                  % if 0, no file output but gnuplot display
    svg_en,...
    subplot_en,...                  % "power"/"dispersion"/"both"
    snr_plot_en,...                 % enable ESNR plot
    filename...                     % name of output file
)

% Number of given paramters must exactly match the num of possible parameters
if nargin < nargin('amplified_fiber_link')
    error('Too few arguments.');
end
if nargin > nargin('amplified_fiber_link')
    error('Too many arguments.');
end

% Parameter checks
% TODO
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BEGIN, set up graphical output

% Clear figure objects (useful when executed at octave prompt)
clf;

% Create an invisible figure.
if (webdemo_en)
    fig = figure(1); set(fig, 'visible', 'off');
else
    fig = figure(1); set(fig, 'visible', 'on');
end

landscape = '-S930,350';
portrait  = '-S640,480';

output_format = landscape;     % default

set(0, 'defaultlinelinewidth', 2);

if (webdemo_en)
    set(0, 'defaultaxesfontsize', 8);
    set(0, 'defaultaxesfontname', 'Helvetica');
    set(0, 'defaulttextfontsize', 8);
    set(0, 'defaulttextfontname', 'Helvetica');
    const_lw = 2;
elseif (svg_en)
    set(0, 'defaultaxesfontsize', 6);
    set(0, 'defaultaxesfontname', 'Helvetica');
    set(0, 'defaulttextfontsize', 6);
    set(0, 'defaulttextfontname', 'Helvetica');
    const_lw = 1;
else
    set(0, 'defaultaxesfontsize', 10);
    set(0, 'defaultaxesfontname', 'Helvetica');
    set(0, 'defaulttextfontsize', 8);
    set(0, 'defaulttextfontname', 'Helvetica');
    const_lw = 2;
end

% Set default color order for multiline plots
set(0, 'DefaultAxesColorOrder', unicmap);


% END, set up graphical output
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Unit/dimension conversions
total_length = total_length * 1e3; % km to m
span_length = span_length * 1e3; % km to m
attenuation_coefficient = attenuation_coefficient / 1e3; % dB/km to dB/m
dispersion_coefficient = dispersion_coefficient / 1e6; % ps/(km*nm) to s/(m*m)
dispersion_slope = dispersion_slope  * 1e3; % ps/(km*nm*nm) to s/(m*m*m)
dcf_dispersion_coefficient = dcf_dispersion_coefficient / 1e6; % ps/(km*nm) to s/( m*m)
dcf_dispersion_slope = dcf_dispersion_slope * 1e3; % ps/(km*nm*nm) to s/(m*m*m)
dcf_attenuation_coefficient = dcf_attenuation_coefficient / 1e3; % dB/km to dB/m

% Determine plot enable vector for plot_parameter_evolution()
if(strcmp(subplot_en, 'power'))
    plot_ena = [1 0];
elseif(strcmp(subplot_en, 'dispersion'))
    plot_ena = [0 1];
else
    % Default
    plot_ena = [1 1];
end

% Plot parameter evolution
if plot_en(1)
    simple_plot = 1;
    snr_plot_en = 0;
    plot_parameter_evolution(
        total_length,...                
        span_length,...                 
        sampling_length,...             
        input_power,...                 
        attenuation_coefficient,...     
        noise_figure,...                
        center_wavelength,...           
        no_wdm_channels,...             
        channel_spacing,...             
        dispersion_coefficient,...      
        dispersion_slope,...            
        dcf_dispersion_coefficient,...  
        dcf_dispersion_slope,...        
        dcf_attenuation_coefficient,... 
        amplification,...               
        dispersion_compensation,
        plot_ena,
        snr_plot_en,
        simple_plot );

    % Set aspect ratio of image file
    output_format = portrait;
end

% Plot shot/ASE noise, electr. signal power, OSNR, dispersion
if plot_en(2)
    simple_plot = 0;
    plot_parameter_evolution(
        total_length,...                
        span_length,...                 
        sampling_length,...             
        input_power,...                 
        attenuation_coefficient,...     
        noise_figure,...                
        center_wavelength,...           
        no_wdm_channels,...             
        channel_spacing,...             
        dispersion_coefficient,...      
        dispersion_slope,...            
        dcf_dispersion_coefficient,...  
        dcf_dispersion_slope,...        
        dcf_attenuation_coefficient,... 
        amplification,...               
        dispersion_compensation,
        plot_ena,
        snr_plot_en,
        simple_plot );

    % Set aspect ratio of image file
    output_format = portrait;
end

% Create image file that will be displayed on the webpage. Its size
% must be 640x480 for "normal" and 1280x960 for wide demo slides.
if(webdemo_en)
    print(filename, '-dpng',output_format);
end
if (svg_en)
    print(['fig_',filename], '-dsvg',output_format);
end

end
% END main function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BEGIN Function definitions
function plot_parameter_evolution(...
    total_length,...                % in m
    span_length,...                 % in m
    sampling_length,...             % in m
    input_power,...                 % in dBm
    attenuation_coefficient,...     % in dB/m
    noise_figure,...                % in dB
    center_wavelength,...           % in m
    no_wdm_channels,...             % integer number
    channel_spacing,...             % in Hz
    dispersion_coefficient,...      % in s / (m*m)
    dispersion_slope,...            % in s / (m*m*m)
    dcf_dispersion_coefficient,...  % in s / (m*m)
    dcf_dispersion_slope,...        % in s / (m*m*m)
    dcf_attenuation_coefficient,... % in dB/m
    amplification,...               % boolean
    dispersion_compensation,...     % boolean
    plot_ena,...                    % Enable subplots [power_plot dispersion_plot ]
    snr_plot_ena,...                % Enable (E)SNR plot
    simple_ena                      % Enable simple OSNR plot (opt. power and OSNR)
)

% default value of simple_ena
if nargin < 19
    simple_ena = 0;
end

% Call to fiberParameterEvolution()
[l, P, P_s, P_a, osnr_dB, P_t, D_acc, lm_min, lm_max, snr_dB, G, Q, l_DCF] = ...
fiberParameterEvolution(...
    total_length,...
    span_length,...
    sampling_length,...
    input_power,...
    attenuation_coefficient,...
    noise_figure,...
    center_wavelength,...
    no_wdm_channels,...
    channel_spacing,...
    dispersion_coefficient,...
    dispersion_slope,...
    dcf_dispersion_coefficient,...
    dcf_dispersion_slope,...
    dcf_attenuation_coefficient,...
    amplification,...
    dispersion_compensation );

% Dimension conversion
l_km = l / 1e3;    
    
%% Power budget plot
if(plot_ena(1))
    if(plot_ena(2))
        subplot(2,1,1); % Two plots, create subplot
        plot_both = 1;
    else
        plot_both = 0;
    end

    % Plot electrical SNR including thermal noise?
    if snr_plot_ena
        second_yaxis_data = [osnr_dB; snr_dB];
    else
        second_yaxis_data = osnr_dB;
    end

    if simple_ena
        [ah_vec, hl1_vec, hl2_vec] = plotyy(l_km, Q, l_km, osnr_dB);
    else
        [ah_vec, hl1_vec, hl2_vec] = plotyy(l_km, [P; P_s; P_a], l_km,  second_yaxis_data);
        set(hl1_vec(2), 'Color', unicmap(5)); % different color for shot noise curve
    end
    
    
    % Style of OSNR curve
    set(hl2_vec(1), 'Marker', 'o', 'Markersize', 2, 'Linewidth', 4);
    
    % Style of SNR curve
    if snr_plot_ena
        set(hl2_vec(2), 'Color', get(hl2_vec(1), 'Color')); % Same Color for OSNR and SNR
        set(hl2_vec(2), 'LineStyle', '--'); % Line style SNR
    end
    
    % Left y axis    
    if simple_ena
        set(get(ah_vec(1), 'YLabel'), 'String', 'opt. signal power in dBm');
        ylims_left = get(ah_vec(1), 'Ylim');
        margin_rel = 0.1;
        y_size = diff(ylims_left);
        ylims_new = [ylims_left(1) - margin_rel*y_size,...
                     ylims_left(2) + margin_rel*y_size];
        set(ah_vec(1), 'YLim', ylims_new);
    else
        % thermal noise floor
        set(gcf, 'currentaxes', ah_vec(1));
        line(l_km, ones(1,length(l_km))*P_t, 'LineStyle', '--', 'Color', 0.2*[1 1 1]); 
        
        set(get(ah_vec(1), 'YLabel'), 'String', 'el. power in dBm');    
        set(ah_vec(1), 'YColor', [0 0 0]);
    end

    % Right y axis
    if(~simple_ena && snr_plot_ena)
        ylabel_str = '(E/O)SNR in dB';
    else
        ylabel_str = 'OSNR in dB';
    end
    set(get(ah_vec(2), 'YLabel'), 'String', ylabel_str);
    ylims_right = get(ah_vec(2), 'Ylim');
    margin_rel = 0.1; % ensure some extra margin of right y axis size for visibility of
                      % ellipse/arrow marker
    if simple_ena
        y_low = (osnr_dB(end) - margin_rel*ylims_right(2))/(1-margin_rel);
    else
        y_low = (min([osnr_dB, snr_dB])...
            - margin_rel*ylims_right(2))/(1-margin_rel);
    end
    set(ah_vec(2), 'Ylim', [y_low ylims_right(2)]);
    
    % x axis    
    set(ah_vec, 'XLim', [min(l_km) max(l_km)*1.05]);
    set(ah_vec(1), 'XTickLabel', {}); % remove x-tick label duplicate
    if(plot_both)
        set(ah_vec(2), 'XTickLabel', {});
    else
        set(get(ah_vec(1), 'XLabel'), 'String', 'length in km');
    end
    
    % Grid
    grid on;
            
    % Labelling
    %%%%%%%%%%%
    if(plot_both)
        t_y_pos = 1.1;  % normalized figure units
    else
        t_y_pos = 1.03; % normalized figure units
    end
    t_x_pos = 0;        % normalized figure units

    if ~simple_ena
        th = text(t_x_pos, t_y_pos, 'Legend:', 'Color', [0 0 0], 'FontSize', 7,...
            'Units', 'Normalized', 'HorizontalAlignment', 'left');
    end
    
    t_x_pos = t_x_pos + 0.14;
    if ~simple_ena
        th = text(t_x_pos, t_y_pos, 'signal', 'Color', unicmap(1), 'FontSize', 7,...
            'Units', 'Normalized', 'HorizontalAlignment', 'left');
    end
    
    t_x_pos = t_x_pos + 0.1;
    if ~simple_ena
        th = text(t_x_pos, t_y_pos, 'shot, ', 'Color', unicmap(5), 'FontSize', 7,...
            'Units', 'Normalized', 'HorizontalAlignment', 'left');
    end
    
    t_x_pos = t_x_pos + 0.08;
    if ~simple_ena
        th = text(t_x_pos, t_y_pos, 'ASE noise', 'Color', unicmap(3), 'FontSize', 7,...
            'Units', 'Normalized', 'HorizontalAlignment', 'left');
    end
    
    t_x_pos = t_x_pos + 0.2;
    if(plot_both)
        y_delta = [-0.04, 0.05];
    else
        y_delta = [-0.01, 0.03];
    end
    th = text(t_x_pos, t_y_pos + y_delta(1), sprintf('G=%.1f dB', G), 'Color',...
        0.4*[1 1 1], 'FontSize', 6, 'Units', 'Normalized',...
        'HorizontalAlignment', 'left');
    th = text(t_x_pos, t_y_pos + y_delta(2), sprintf('l_{DCF}=%.1f km', l_DCF/1e3),...
        'Color', 0.4*[1 1 1], 'FontSize', 6, 'Units', 'Normalized',...
        'HorizontalAlignment', 'left');        
    
    t_x_pos = t_x_pos + 0.18;
    th = text(t_x_pos, t_y_pos, ['OSNR at RX: ' sprintf('%.1f', osnr_dB(end)) ' dB'],...
        'Color', unicmap(2), 'FontSize', 7,...
        'Units', 'Normalized', 'HorizontalAlignment', 'left');
    
    l_xlims = get(ah_vec(1), 'XLim');
    P_ylims = get(ah_vec(1), 'YLim');
    if(P_a(end-1) > P_t)
        t_x_pos = (l_xlims(2) - l_xlims(1))*0.7;
        t_y_pos = P_t + 0.02*(P_ylims(2) - P_ylims(1));
        v_align = 'Bottom';
    else
        t_x_pos = (l_xlims(2) - l_xlims(1))*0.025;
        t_y_pos = P_t;
        v_align = 'Top';
    end
    if ~simple_ena
        text(t_x_pos, t_y_pos, 'Thermal noise floor',...
            'FontSize', 6, 'VerticalAlignment', v_align);
    end


end % if(plot_ena(1))...

%% Dispersion budget plot
if(plot_ena(2))
    if(plot_ena(1))
        subplot(2,1,2); % Two plots, create subplot
    end
    lh_vec = plot(l_km, D_acc/1e9*1e12);
    set(lh_vec(1), 'Color', unicmap(9)); % min wavelength: violet
    set(lh_vec(2), 'Color', unicmap(3)); % center wavelength: green
    set(lh_vec(3), 'Color', unicmap(2)); % max wavelength: red
    xlim([min(l_km) max(l_km)*1.05]);
    title('');
    ylabel('Dispersion in ps/nm');
    xlabel('length in km');
    yticks = get(gca, 'YTick');

    % Reduce the gap between the two plots
    if(plot_ena(1))
        lower_pos = get(gca, 'Position');
        upper_pos = get(ah_vec(1), 'Position');
        gap_size = upper_pos(2) - lower_pos(4) - lower_pos(2);
        epsi = gap_size/8; % don't close gap completely
        enlargement = gap_size/2 - epsi;
        % Change position of dispersion plot
        set(gca, 'Position', [lower_pos(1:3) lower_pos(4) + enlargement]);
        % Change position of power/OSNR plot
        set(ah_vec(1), 'Position', [upper_pos(1), ...
            upper_pos(2) - enlargement, upper_pos(3),...
            upper_pos(4) + enlargement]);
    end

    grid on;
    
    legend(sprintf('%d nm',round(lm_min/1e-9)),...
            sprintf('%d nm',round(center_wavelength/1e-9)),...
            sprintf('%d nm',round(lm_max/1e-9)),...
            'Location', 'Northwest');
end

% Marker to show (E/O)SNR at RX
% Plotting of marker is done at last because rescaling of axes would also 
% rescale and distort the marker.
if(plot_ena(1))
    ellipse_a_px = 10;
    ellipse_b_px = ellipse_a_px;
    arrow_size_px = 5;
    ellipse_color = [0 0 0];
    draw_ellipse(l_km, osnr_dB, ah_vec(2), ellipse_a_px, ellipse_b_px,...
        arrow_size_px, ellipse_color);
    if(~simple_ena && snr_plot_ena)    
        draw_ellipse(l_km, snr_dB, ah_vec(2), ellipse_a_px, ellipse_b_px,...
            arrow_size_px, ellipse_color);
    end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ l, P, P_s, P_a, osnr_dB, P_t, D_acc, lm_min, lm_max, snr_dB, G, Q, l_DCF ] = ...
fiberParameterEvolution(...
    total_length,...                % in m
    span_length,...                 % in m
    sampling_length,...             % in m
    input_power,...                 % in dBm
    attenuation_coefficient,...     % in dB/m
    noise_figure,...                % in dB
    center_wavelength,...           % in m
    no_wdm_channels,...             % integer number
    channel_spacing,...             % in Hz
    dispersion_coefficient,...      % in s / (m*m)
    dispersion_slope,...            % in s / (m*m*m)
    dcf_dispersion_coefficient,...  % in s / (m*m)
    dcf_dispersion_slope,...        % in s / (m*m*m)
    dcf_attenuation_coefficient,... % in dB/m
    amplification,...               % boolean
    dispersion_compensation)        % boolean
%FIBERPARAMETEREVOLUTION Calculate noise power and dispersion on an amplified fiber link.
%   
%  Calculate the evolution of signal power, shot noise power, ASE noise
%  power, OSNR and accumulated dispersion of a signal transmitted over an
%  amplified WDM fiber link with dispersion compensation. The link is
%  divided into spans. Each span consists of a standard single mode fiber
%  (SMF), dispersion compensating fiber (DCF) and an amplifier. DCF and
%  amplifier are optional. The DCF is inserted before amplification and the
%  amplifier exactly compensates for the power loss caused by SMF and DCF.
%  Similarly the length of the DCF is chosen is such a way that it exactly
%  compensates for 
%
%  IN
%   total_length - Length of fiber link in m.
%   span_length - Length of SMF of one span in m.
%   sampling_length - All output parameters are calculated at equidistant
%                     points on the link. This is the distance between two
%                     points.
%   input_power - Optical launch power in dBm.
%   attenuation_coefficient - Attenuation coefficient of SMF in dB/m.
%   noise_figure - Optical amplifier noise figure in dB.
%   center_wavelength - Wavelength of the center WDM channel in m.
%   no_wdm_channels - Number of WDM channels.
%   channel_spacing - Bandwidth of one WDM channel in Hz.
%   dispersion_coefficient - Dispersion coefficient of the SMF at 1550nm 
%                            wavelength in s / (m*m).
%   dispersion_slope - Dispersion slope of the SMF at 1550nm wavelength in 
%                      s / (m*m*m).
%   dcf_dispersion_coefficient - Dispersion coefficient of the DCF at
%                                1550nm wavelength in s / (m*m).
%   dcf_dispersion_slope - Dispersion slope of the DCF at 1550nm wavelength 
%                          in s / (m*m*m).
%   dcf_attenuation_coefficient - Attenuation coefficient of the DCF in
%                                 dB/m.
%   amplification - Set to 1/0 to switch amplifciation on/off.
%   dispersion_compensation - Set to 1/0 to switch dispersion compensation
%                             on/off.
%
%  OUT
%   l - Vector of lengthes in m gives the distances at which the other
%       output parameters have been evaluated.
%   P - Power of current that is caused by the optical signal in an ideal
%       photodetector. Unit: dBm. This vector has the same length as vector 
%       l.
%   P_s - Shot noise power in an ideal photodetector in an electrical 
%         bandwidth of 12.5/2 GHz. Unit: dBm. This vector has the same 
%         length as vector l.
%   P_a - Power of current that is caused by the ASE noise in an ideal
%         photodetector in a bandwidth of 12.5/2 GHz. Unit: dBm. This 
%         vector has the same length as vector l.
%   osnr_dB - The signal to noise ratio calculated by the above signal and
%             noise powers. This vector has the same length as l.
%   P_t - Thermal noise power in the electrical receiver in a bandwidth of
%         of 12.5/2 GHz. Unit: dBm. This is a scalar value.
%   D_acc - Row vector with three rows containing the accumulated
%           dispersion on the WDM channels with lowest, middle and highest
%           wavelength, respectively. Unit: dBm. This vector has the same 
%           length as l.
%   lm_min - Wavelength of channel with lowest wavelength. Unit: m.
%   lm_max - Wavelength of channel with biggest wavelength. Unit: m.
%   snr_dB - Electrical signal-to-noise ratio after photodetection including
%            thermal noise in a bandwidht of 12.5/2 GHz. Unit: dB.
%   G - Amplifier gain in dB.
%   Q - Optical signal power in dBm. Vector has same length as vector l.
%   l_DCF - Length of dispersion compensating fiber. Unit: m.

% Convention:
% ----------
% In the following, all optical powers are denoted Q
%                all electrical powers are denoted P

% constants
c = 299792458;       % vacuum speed of light in m/s
h = 6.62606957e-34;  % planck constant in Js
q = 1.602176565e-19; % electron charge in C (Coulomb)
k_B = 1.3806488e-23; % Boltzmann's constant in J/K

% conversions
nu = c / center_wavelength; % center frequency in 1/s
input_power_lin = 10^(input_power/10) * 1e-3; % in W

% Fixed parameters
quantum_efficiency = 1; % Quantum efficiency of photodetector
responsivity = quantum_efficiency * q / h / nu;
electrical_bandwidth = 12.5e9 / 2; % in 1/s
temperature = 290; % receiver temperature in K

no_spans = ceil(total_length / span_length);
no_length_samples = floor(span_length/sampling_length);
max_wavelength = c / (nu - ceil(no_wdm_channels / 2 - 1) * channel_spacing);
min_wavelength = c / (nu + ceil(no_wdm_channels / 2) * channel_spacing);
% dispersion coefficient of SMF at min/center/max wavelength
max_dispersion = dispersion_coefficient + ...
    dispersion_slope * (max_wavelength - center_wavelength);
min_dispersion = dispersion_coefficient + ...
    dispersion_slope * (min_wavelength - center_wavelength);
dispersion_coeffs = [min_dispersion; dispersion_coefficient; max_dispersion];
% dispersion coefficient of DCF at min/center/max wavelength
max_dcf_dispersion = dcf_dispersion_coefficient + ...
    dcf_dispersion_slope * (max_wavelength - center_wavelength);
min_dcf_dispersion = dcf_dispersion_coefficient + ...
    dcf_dispersion_slope * (min_wavelength - center_wavelength);
dcf_dispersion_coeffs = [min_dcf_dispersion; dcf_dispersion_coefficient; max_dcf_dispersion];

% Initialization
l = 0; % length in m
Q = input_power; % signal power in dBm (opt.)
P_s = calcShotNoise(Q, center_wavelength, electrical_bandwidth,...
    quantum_efficiency); % shot noise power in dBm (el.)
P_a = -inf; % ASE noise power in dBm (el.)
D_acc = [0; 0; 0]; % accumulated dispersion (1st row: min wavel., 2nd: center, 3rd: max)

for spi = 1:no_spans
    for li = 1:no_length_samples
       % Attenuation on fiber
       l(end+1) = l(end) + sampling_length;
       Q(end+1) = Q(end) - sampling_length * attenuation_coefficient;       
       P_a(end+1) = P_a(end) - 2 * sampling_length * attenuation_coefficient;
       Q_tot = 10*log10(10^(Q(end)/10) + sqrt(10^(P_a(end)/10))/responsivity);
       P_s(end+1) = calcShotNoise(Q_tot, center_wavelength,...
           electrical_bandwidth, quantum_efficiency);
       D_acc(:, end+1) = D_acc(:, end) + sampling_length .* dispersion_coeffs; 
    end
    % Remaining fiber length
    l_diff = span_length - no_length_samples * sampling_length;
    if (l_diff ~= 0)
        l(end+1) = l(end) + l_diff;
        Q(end+1) = Q(end) - l_diff * attenuation_coefficient;
        P_a(end+1) = P_a(end) - 2*l_diff * attenuation_coefficient;
        Q_tot = 10*log10(10^(Q(end)/10) + sqrt(10^(P_a(end)/10))/responsivity);
        P_s(end+1) = calcShotNoise(Q_tot, center_wavelength,...
           electrical_bandwidth, quantum_efficiency);
        D_acc(:, end+1) = D_acc(:, end) + sampling_length .* dispersion_coeffs; 
    end
    
    % Amplification and dispersion are considered as "lumped" -> no distance
    % is added
    l(end+1) = l(end);
    
    % Dispersion compensation
    if(dispersion_compensation)
        dcf_length = - dispersion_coefficient / dcf_dispersion_coefficient...
            * span_length;
        dcf_attenuation = dcf_attenuation_coefficient * dcf_length;
        D_acc(:, end+1) = D_acc(:, end) + dcf_length * dcf_dispersion_coeffs;        
    else
        dcf_length = 0;
        dcf_attenuation = 0;
        D_acc(:, end+1) = D_acc(:, end);
    end
    Q(end+1) = Q(end) - dcf_attenuation;
    P_a(end+1) = P_a(end) - 2*dcf_attenuation;
    Q_tot = 10*log10(10^(Q(end)/10) + sqrt(10^(P_a(end)/10))/responsivity);
    P_s(end+1) = calcShotNoise(Q_tot, center_wavelength,...
        electrical_bandwidth, quantum_efficiency);
    
    % Amplification (compensates for attenuation of one span)
    l(end+1) = l(end);
    if(amplification)
        gain = span_length * attenuation_coefficient + dcf_attenuation; % in dB
        P_amp = noiseFigureToAseNoisePower(noise_figure, gain, Q(end),...
            quantum_efficiency, center_wavelength, electrical_bandwidth);
        P_a(end+1) = 10*log10(10^((P_a(end) + 2*gain)/10) + 10^(P_amp/10));
        Q(end+1) = Q(end) + gain; % amplification
    else
        gain = 0;
        Q(end+1) = Q(end);
        P_a(end+1) = P_a(end);
    end    
    Q_tot = 10*log10(10^(Q(end)/10) + sqrt(10^(P_a(end)/10))/responsivity);
    P_s(end+1) = calcShotNoise(Q_tot, center_wavelength,...
        electrical_bandwidth, quantum_efficiency);
    D_acc(:,end+1) = D_acc(:,end); % amp adds no dispersion
end

% electrical signal power
Q_lin = 10.^(Q/10)*1e-3; % dBm to Watt
P_lin = Q_lin.^2 * responsivity^2;
P = 10*log10(P_lin/1e-3); % Watt to dBm

% Calculate OSNR
% P_lin = 10.^(P/10)*1e-3;
P_a_lin = 10.^(P_a/10)*1e-3;
P_s_lin = 10.^(P_s/10)*1e-3;
% the commented line holds for _optical_ powers
% osnr = P_lin.^2 * responsivity.^2 ./ (P_a_lin + P_s_lin); 
osnr = P_lin ./ (P_a_lin + P_s_lin);
osnr_dB = 10*log10(osnr);

% Thermal noise of electric receiver
P_t_lin = k_B * temperature * electrical_bandwidth;
P_t = 10*log10(P_t_lin / 1e-3);

% Electrical SNR
snr = P_lin ./ (P_a_lin + P_s_lin + P_t_lin);
snr_dB = 10*log10(snr);

% Short symbols for return values
lm_min = min_wavelength;
lm_max = max_wavelength;
G = gain;
l_DCF = dcf_length;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function ase_noise_power_dBm = noiseFigureToAseNoisePower(...
    noise_figure,...       % in dB
    gain,...               % amp. gain in dB
    input_power,...        % optical input power in dBm
    quantum_efficiency,... % Quantum efficiency of photodetector
    wavelength,...         % carrier wavelength in m
    bandwidth )            % electrical bandwidth in Hz    
%NOISEFIGURETOASENOISEPOWER Convert optical amplifier noise figure to ASE noise power.
%
%   Calculate the power of the current that is caused by ASE noise in an
%   ideal photodetector. Output is in dBm.

if nargin < 6
    bandwidth = 12.5e9 / 2;
end

if nargin < 5
    wavelength = 1550e-9;
end

if nargin < 4
    quantum_efficiency = 1;
end

% constants
c = 299792458;       % vacuum speed of light in m/s
h = 6.62606957e-34;  % planck constant in Js
q = 1.602176565e-19; % electron charge in C (Coulomb)

% conversions
F = 10^(noise_figure/10); % noise factor (linear scale)
nu = c / wavelength; % carrier frequency in Hz
P_in = 10^(input_power/10) * 1e-3; % input power in W
G = 10^(gain/10); % gain, linear scale
responsivity = quantum_efficiency * q / h / nu; % in A/W

P_Ase = (F*G^2 -G)*2*responsivity^2*h*nu*P_in*bandwidth;
ase_noise_power_dBm = 10*log10(P_Ase/1e-3);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function shot_noise_power_dBm = calcShotNoise(...
    optical_signal_power,... % in dBm
    wavelength,...           % in m
    electrical_bandwidth,... % in 1/s
    quantum_efficiency )     % Quantum efficiency of photodetector
%CALCSHOTNOISE Calculate shot noise power for a given optical signal power.
%   

% constants
c = 299792458;       % vacuum speed of light in m/s
h = 6.62606957e-34;  % planck constant in Js
q = 1.602176565e-19; % electron charge in C (Coulomb)

% conversions
signal_power_lin = 1e-3 * 10^(optical_signal_power/10); % in W
nu = c / wavelength; % frequency in 1/s
responsivity = quantum_efficiency * q / h / nu; % in A/W

shot_noise_power = 2 * responsivity^2 * h * nu * signal_power_lin *...
    electrical_bandwidth;
shot_noise_power_dBm = 10*log10(shot_noise_power/1e-3);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ uni_colors ] = unicmap(n)
%UNICMAP Get a matrix of colors defined by the University of Stuttgart
%   unicmap returns a matrix of 9 colors defined by the University of 
%           Stuttgart.
%   unicmap(n) returns the University color with index n. n must be in
%              [1..9].

% uni_hellgrau = [224/255 224/255 224/255];
% uni_mittelgrau = [128/255 128/255 128/255];
% uni_dunkelgrau = [80/255 80/255 80/255];    

narginchk(0, 1);

uni_colors = [0 126/255 198/255]; % mittelblau
uni_colors = [uni_colors; 238/255 28/255 35/255]; % rot
uni_colors = [uni_colors; 140/255 198/255 62/255]; % apfelgruen
uni_colors = [uni_colors; 1 221/255 0]; % gelb
uni_colors = [uni_colors; 244/255 111/255 33/255]; % orange
uni_colors = [uni_colors; 237/255 0 140/255]; % pink
uni_colors = [uni_colors; 19/255 31/255 31/255]; % anthrazit
uni_colors = [uni_colors; 0/255 170/255 173/255]; % turquoise
uni_colors = [uni_colors; 128/255 40/255 145/255]; % lila

if nargin == 1
    if (~isscalar(n) || ~isnumeric(n))
        error('First argument must be a numeric scalar.')
    end
    if (n < 1 || n > size(uni_colors, 1))
        error(['First argument mus be postive and less than or equal to '...
            num2str(size(uni_colors, 1))]);
    end
    uni_colors = uni_colors(n,:);
end

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function draw_ellipse(x_vec, y_vec, ah, a_px, b_px, ar_size_px, color)

% default color
if nargin < 6    
    color = [0 0 0];
end

% Save current axes handle and Y limits
fh = get(ah, 'Parent');
current_ah = get(fh, 'currentaxes');
current_ylims = get(ah, 'YLim');
set(fh, 'currentaxes', ah);

% Get axes and figure position
apos = get(ah, 'Position');                % in normalized figure units

current_figure_units = get(fh, 'units');
set(fh, 'units', 'pixels');
fpos = get(get(ah, 'Parent'), 'Position'); % in pixels

% Length of x/y axis in data units
xlims = get(ah, 'XLim'); % in km
ylims = get(ah, 'YLim'); % in dB

a_height_px = apos(4)*fpos(4); % axes height in pixels
a_width_px = apos(3)*fpos(3);  % axes width in pixels
w = xlims(2) - xlims(1);       % length of x axis in data units
h = ylims(2) - ylims(1);       % length of y axis in data units

% Draw ellipse
x_0 = x_vec(end);
y_0 = y_vec(end);
a = a_px/a_width_px * w;
b = b_px/a_height_px * h;
t = linspace(0,2*pi,100);
x = a*cos(t) + x_0;
y = b*sin(t) + y_0;
lh = line(x, y);
set(lh, 'Color', color);

% Draw arrow shaft (straight line)
x_start = max(x);
x_end = xlims(2);
y_start = y(x == x_start);
y_start = y_start(1); % remove duplicate values
y_end = y_start;
lh = line([x_start x_end], [y_start y_end]);

% Draw arrowhead
y_ar_start = y_end;
x_ar_start = x_end;
ar_size = ar_size_px/a_height_px * h;                      % in dB
x_ar_end = x_ar_start - ar_size *(apos(4)/apos(3))*(w/h);  % in km
line([x_ar_start x_ar_end], [y_ar_start y_ar_start + ar_size]);
line([x_ar_start x_ar_end], [y_ar_start y_ar_start - ar_size]);

% Restore current axes, Y limits and figure units
set(ah, 'YLim', current_ylims);
set(fh, 'currentaxes', current_ah);
set(fh, 'units', current_figure_units);

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% END Function definitions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%