% WEBDEMO modal_dispersion_spectrum.m-file, Institute of Telecommunications, University of Stuttgart, Germany
% www.inue.uni-stuttgart.de
% author: Michael Bernhard, michibernhard@gmx.de
% date: 17.10.2013
% EXAMPLE:  modal_dispersion_spectrum(10, 1, 1.46, 1.445, 50, 1550, 100, 0, 0, "modal_dispersion")
% PARAMETERS:
%  fiber_length,     	% in km
%  symbol_rate,      	% in GHz  
%  n_K,			% Brechzahl Kern
%  n_M,			% Brechzahl Mantel
%  d,			% diameter in µm
%  lambda,		% in nm, wavelength
%  attenuation, % in dB/km
%  nb_points,   	% number of points for plots
%  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_spectrum(
  fiber_length,    	% in km
  symbol_rate,    	% in GHz  
  n_K,		        % Brechzahl Kern
  n_M ,			% Brechzahl Mantel
  d,     		% diameter in µm
  lambda,		% in nm, wavelength
  attenuation, % in dB/km
  nb_points,            % number of points
  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; 


  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
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


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


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN check input parameters for consistency

  %fiber_length = fiber_length * 1e3; % from km to m
  symbol_rate = symbol_rate * 1e6; % from MHz to Hz
  lambda = lambda * 1E-9; % from nm to m
  %attenuation = attenuation / 1000; % in dB/m
  a = d / 2 * 1E-6; % from µm to m and diameter to radius

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



  % Constants:
  c = 299792458; % speed of light
  hold off;
  impulse_shaper = 5;
  oversampling_factor = 100;
  f_S = 1E6;
  t1 = 1/f_S; % symbol time
  T_d = t1/oversampling_factor; % smallest discret time resolution

  A_N = sqrt(n_K^2 - n_M^2); % numerical aperture
  V = 2*pi / lambda * a * A_N; % normalized frequency parameter
  M = floor(4*(V/pi)^2); % number of modes
  gamma_crit = acos(n_M/n_K);
  if M==1
    gamma = 0;
  else
    if M > 1000
      M = 1000;
    endif
    gamma = linspace(-gamma_crit, gamma_crit, M); % equally distributed angles
  endif
  lengths = fiber_length ./ cos(gamma);

  
  t = lengths / c / n_K; % calculate time delay
  t_min = min(t);
  t_max = max(t);

  discrete_time = linspace(t_min, t_max, nb_points);
  h = zeros(1,nb_points); % impulse response h(t)
  for k=1:ceil(M/2)
    if M > 2
      index = floor((nb_points - 1) / (t_max - t_min) * (t(k) - t_min) + 1);    
    else
      index = 1;
    endif
    if gamma(k) == 0
      h(1, index) = 1/M.* exp(-attenuation .* lengths(k));
    else
      h(1,index) = 2/M * cos(gamma(k)).* exp(-attenuation .* lengths(k));
    endif

  endfor

  sfh = subplot(1,2,1);
  set(sfh,'outerposition',[0 0 0.5 1]);
  stem(discrete_time*1e6, h, "linewidth", const_lw, "color", "b");
  xlabel('time t/ms');
  ylabel('h(t)');
  title('');


  R_s = symbol_rate;
  f = linspace(-R_s, R_s, nb_points+1);
 
  X=zeros(1,nb_points+1);
  for k=1:M
    X = X + exp(1j*2*pi.*f*t(k)*R_s/R_s) * cos(gamma(k)) .* exp(-attenuation .* lengths(k));
  endfor
  X = X / M;

  sfh = subplot(1,2,2)
  
  plot_spectrum(X,f, R_s, const_lw, [0 0 1], [0 0 0]);

  hold off

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

end

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






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



function plot_spectrum(X, f, R_s, const_lw, color_magnitude, color_phase)
  frequency_scaling = 1E6; % unit in MHz instead of Hz
  magnitude = 20*log10(abs(X));
  min_magnitude = min(min(magnitude))-1;
  max_magnitude = max(max(magnitude))+1;

  N = size(X,2);

  phase = angle(X);
  phase = unwrap(phase);
  phase = phase/pi;

  magenta = [[51,0,51]; [102,0,102]; [153,0,153]; [204,0,204]; [255,0,255]; [255,51,255]; [255,102,255]; [255,153,255]; [255,204,255]]/255;
  rectangle("Position",[-R_s/frequency_scaling,min_magnitude+0.01, R_s/2/frequency_scaling,max_magnitude-min_magnitude],"FaceColor",magenta(7,:),"EdgeColor",magenta(7,:)); % add +0.01 to min_magnitude, because of a bug of octave
  hold on;
  rectangle("Position",[R_s/2/frequency_scaling,min_magnitude+0.01, R_s/2/frequency_scaling,max_magnitude-min_magnitude],"FaceColor",magenta(7,:),"EdgeColor",magenta(7,:));
 
  % plot magnitude and phase into the same plot
  [AX,H1, H2] = plotyy(f/frequency_scaling, magnitude,f/frequency_scaling,phase);

  set(get(AX(1),'Ylabel'),'String','|H(f)|^2 in dB')
  set(get(AX(2),'Ylabel'),'String','phase / \pi')

  xlabel('frequency f/MHz')
 
  set(H1,'LineWidth',const_lw,'Color', color_magnitude);
  set(H2,'LineWidth',const_lw,'Color', color_phase);
  set(AX(1),'YColor', color_magnitude);
  set(AX(2),'YColor', color_phase);
  set(AX(1), 'Position', [0.56 0.15 0.3 0.79]);
  set(AX(2), 'Position', [0.56 0.15 0.3 0.79]);
  ylim([min_magnitude max_magnitude]);
  xlim([-R_s/frequency_scaling R_s/frequency_scaling]);

endfunction



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


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

