% WEBDEMO modal_dispersion_animation.m-file, Institute of Telecommunications, University of Stuttgart, Germany
% www.inue.uni-stuttgart.de
% author: Michael Bernhard, michibernhard@gmx.de
% date: 17.10.2013
%
% Written for GNU Octave 
%
% EXAMPLE:  modal_dispersion(0, 50, 1, 1.46, 1.445, 50, 1550, 5, 0, 0, 0, "modal_dispersion")
% PARAMETERS:
%  fiber_length_start,    % in km
%  fiber_length_end, % in km
%  symbol_rate,     % in GHz  
%  n_c,             % Brechzahl Kern
%  n_cl,             % Brechzahl Mantel
%  d,               % diameter in µm
%  lambda,          % in nm, wavelength
%  impulse_shaper,  % 1: rectangular, 4: raised cosine  
%  // random_data,     % 1: create random data
%  svg_en,          % if 1, then also output svg-file
%  webdemo_en,      % if 0, then test-mode in Octave/Matlab commandline, no need for webserver running
%  filename         % output file for webdemo





% BEGIN, main function (control)

function modal_dispersion_animation(
  fiber_length_start, % in km
  fiber_length_end,   % in km
  symbol_rate,      % in MBd  
  n_c,              % Brechzahl Kern
  n_cl,              % Brechzahl Mantel
  d,                % diameter in µm
  lambda,           % in nm, wavelength
  impulse_shaper,   % 1: rectangular, 4: raised cosine  
  svg_en,           % if 1, then also output svg-file
  webdemo_en,       % if 0, then test-mode in Octave/Matlab commandline, no need for webserver running
  filename          % output file for webdemo
)


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

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


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

  output_format = landscape; %portrait;     % 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;
  endif


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



  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN check input parameters for consistency
  fiber_length_start = fiber_length_start * 1e3; % from km to m
  fiber_length_end = fiber_length_end * 1e3; % from km to m
  symbol_rate = symbol_rate * 1e6; % from MBd to Hz
  lambda = lambda * 1E-9; % from nm to m
  a = d/2 * 1E-6; % from µm to m; from diameter to radius

  nb_frames = 21; % number of frames for animation
  fiber_length_list = linspace(fiber_length_start, fiber_length_end, nb_frames);


 
  % END check input parameters for consistency
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % Constants:
  c = 299792458; % speed of light

 % if random_data == 1
 %   N = 20;
 %   tx_data = zeros(1,N);
 %   for k=1:N
 %     tx_data(1,k) = round(rand);
 %   endfor
 % else
    tx_data = [0 1 0 1 1 0];
 % endif

  for index=0:nb_frames-1
    fiber_length = fiber_length_list(index+1);
    new_filename = [filename '_' num2str(index)];

    A_N = sqrt(n_c^2 - n_cl^2); % numerical aperture
    V = 2*pi / lambda * a * A_N; % normalized frequency parameter
    M = floor(4*(V/pi)^2); % number of modes
    gamma_L_crit = asin(A_N);
    gamma_L_crit_grad = gamma_L_crit/pi*180;

    gamma_crit = acos(n_cl/n_c);% asin(n_c * sin(gamma_L_crit));
    gamma_crit_grad = gamma_crit/pi*180;
    if M == 1
      gamma = 0;
    else 
      gamma = linspace(-gamma_crit, gamma_crit, M); % Winkel gleichverteilen
    endif;
    lengths = fiber_length ./ cos(gamma);

    hold off;
    oversampling_factor = 1000;
    R_s = symbol_rate;
    t1 = 1/R_s; % symbol time
    T_d = t1/oversampling_factor;% smallest discret time resolution

    nb_points = floor(t1/T_d)*size(tx_data,2);

    [tx_signal_time tx_signal_data] = create_tx_signal(tx_data, impulse_shaper, t1, nb_points);
  
    % compute rx signal
    %t = floor(lengths / c * n_c / T_d) * T_d; % calculate time delay and raster it using T_d
    t = lengths / c * n_c;
    t_min = min(t);
    t_max = max(t);
    delta_t_mo = t_max - t_min;
    delta_t_mo = n_c * fiber_length / c *(n_c / n_cl - 1);
    rx_signal_time = t_min:T_d: t_max+tx_signal_time(end);

    attenuation = ones(1,M)/M;
    for m = 1:M
      attenuation(1,m) = cos(gamma(1,m))/M;
    endfor

    rx_signal_data = zeros(1, length(rx_signal_time) );
    for m=1:M
      X = floor((t(m) - t_min) / T_d);
      Y = ceil(length(rx_signal_time) - length(tx_signal_data) - X);
      rx_signal_data = rx_signal_data + [zeros(1, X) tx_signal_data zeros(1,Y)]*attenuation(m);
    endfor

    tx_signal_data = [tx_signal_data zeros(1,length(rx_signal_data) -length(tx_signal_data))];
    time_axis = [0:length(rx_signal_data)-1]/oversampling_factor;
  
    % plot intensity of tx signal, time axis in ns
    plot(time_axis, tx_signal_data, "linewidth",  const_lw, "color", "b");
    hold on;
    plot(time_axis, rx_signal_data, "linewidth", const_lw, "color", "m");
    title(['R_s = ' num2str(R_s/1E6) ' MBd, l = ' num2str(fiber_length/1E3) ' km, A_N = ' sprintf("%.3f",A_N) ', \gamma_{a, crit} = ' sprintf("%.2f",gamma_L_crit_grad) '^\circ, \gamma_{c, crit} = ' sprintf("%.2f",gamma_crit_grad) '^\circ, M \approx ' num2str(M) ', \Delta t_{Mo} = ' num2str(delta_t_mo*1E9) ' ns'  ]);
    xlim([0  length(tx_data)+1]);
    ylim([0 1.1]);
    xlabel("time t/T");
    ylabel("intensity");
    legend("tx signal", "rx signal");
    hold off

    % Create an image file. This image is displayed on the webpage. Its size
    % must be 640x480.  
    print(new_filename, "-dpng",output_format);
  
    if (svg_en)
      print(["fig_",new_filename], "-dsvg",output_format);
    endif
  endfor
end


% END, main program (control)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BEGIN, function definitions



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create TX signal
function [tx_signal_time tx_signal_data] = create_tx_signal(tx_data, impulse_shaper, t1, nb_points);
  N=size(tx_data,2);
  tx_signal_data = zeros(1,nb_points);
  tx_signal_time = linspace(0,N*t1,nb_points);
  tx_symbol = zeros(1,nb_points/N);


  for k=1:N % symbol index
    for time=0:nb_points/N-1 % time index
      switch (impulse_shaper)
        case 1 % rectangular
          tx_symbol(1,time+1) = tx_data(1,k);
        case 2 % rectangular with filter (see some lines below)
          tx_symbol(1,time+1) = tx_data(1,k);
        case 3 % rectangular with filter (see some lines below)
          tx_symbol(1,time+1) = tx_data(1,k);
        case 4
          t_n = (time - nb_points/N/2)/(nb_points/N/2);
          alpha = 0.5;
          if (abs(t_n) <= (1-alpha)/2)
            tx_symbol(1, time+1) = 1 * tx_data(1,k);
          else
            if (((1-alpha)/2) < abs(t_n) && abs(t_n) <= (1+alpha)/2)
              tx_symbol(1,time+1) = (cos(pi/2/alpha * (abs(t_n)-(1-alpha)/2) ))^2 * tx_data(1,k);
            else
              tx_symbol(1,time+1) = 0;
            endif
          endif;
        case 5
          t_n = (time - nb_points/N/2)/(nb_points/N/2);
          alpha = 1;
          if (abs(t_n) <= (1-alpha)/2)
            tx_symbol(1, time+1) = 1 * tx_data(1,k);
          else
            if (((1-alpha)/2) < abs(t_n) && abs(t_n) <= (1+alpha)/2)
              tx_symbol(1,time+1) = (cos(pi/2/alpha * (abs(t_n)-(1-alpha)/2) ))^2 * tx_data(1,k);
            else
              tx_symbol(1,time+1) = 0;
            endif
          endif;
        otherwise
      endswitch
    endfor
    tx_signal_data(1,1+(k-1)*nb_points/N:k*nb_points/N) = tx_symbol(1,:);
  endfor

  if impulse_shaper == 2
    freq = fft(tx_signal_data);
    M = nb_points/3;
    freq(1,nb_points/2-M:nb_points/2+M) = 0;
    tx_signal_data = ifft(freq);
  endif
  if impulse_shaper == 3
    freq = fft(tx_signal_data);
    M = nb_points/2-10;
    freq(1,nb_points/2-M:nb_points/2+M) = 0;
    tx_signal_data = ifft(freq);
  endif

endfunction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




% END, function definitions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

