% WEBDEMO chromatic_dispersion.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:  chromatic_dispersion(10,5, 17, 1, 0, 0, "chromatic_dispersion")
% PARAMETERS:
%  fiber_length,    % in km
%  symbol_rate,     % in GHz  
%  chromatic_dispersion, % in ps/nm/km   
%  unwrwap_phase,   % if 1 then unwrap phase
%  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 chromatic_dispersion(
  fiber_length,     % in km
  symbol_rate,      % in GHz  
  chromatic_dispersion, % in ps/nm/km  
  unwrap_phase,     % if 1 then unwrap phase
  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 = portrait;  

  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 = fiber_length * 1e3; % from km to m
  symbol_rate = symbol_rate * 1e9; % from GHz to Hz
  chromatic_dispersion = chromatic_dispersion * 1e-6; % from ps/nm/km in s/m/m

 
  % END check input parameters for consistency
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  hold off;

  oversampling_factor = 1;

  N=500;

  c = 299792458; % speed of light in m/s
  lm0 = 1550e-9; % reference wavelength in m
  D = chromatic_dispersion;%17e-6; % chromatic dispersion coefficient in s/m^2 at lm0
  L = fiber_length; %100e3; % fiber length in m
  B = oversampling_factor * symbol_rate; % bandwidth
  f = linspace(-B/2, B/2, N);

  H = exp( j * pi * D * lm0.^2 * L * f.^2 /c);

  frequency_scaling = 1E9*oversampling_factor; % unit in GHz instead of Hz

  magnitude = abs(H).^2;
  phase = angle(H);

  if unwrap_phase == 1
    phase = unwrap(phase) / pi;
  else
    phase = phase / pi;
  endif
 
  % 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')
  set(get(AX(2),'Ylabel'),'String','phase/(\pi)')

  xlabel('frequency f/GHz')

  set(H1,'LineWidth',const_lw,'Color',[1 0 1])
  set(H2,'LineWidth',const_lw,'Color',[0 0 1])
  set(AX(1),'YColor',[1 0 1]);
  set(AX(2),'YColor',[0 0 1]);
  set(AX(1), 'Position', [0.15 0.15 0.7 0.79]);
  set(AX(2), 'Position', [0.15 0.15 0.7 0.79]);

  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 retval = my_fftshift (x, dim)
%http://lists.gnu.org/archive/html/octave-bug-tracker/2011-07/msg00014.html

  retval = 0;

  if (nargin != 1 && nargin != 2)
    print_usage ();
  endif

  if (nargin == 2)
    if (!isscalar (dim))
      error ("fftshift: dimension must be an integer scalar");
    endif
    nd = ndims (x);
    sz = size (x);
    sz2 = ceil (sz(dim) / 2);
    idx = cell ();
    for i = 1:nd
      idx{i} = 1:sz(i);
    endfor
    idx{dim} = [sz2+1:sz(dim), 1:sz2];
    retval = x(idx{:});
  else
    if (isvector (x))
      # Error in Version 3.4.0: x must not be overwritten!
      # x = length (x);  
      # xx = ceil (x/2);
      # retval = x([xx+1:x, 1:xx]);
      # Fix: use xl for the length of x
      xl = length (x);  
      xx = ceil (xl/2);
      retval = x([xx+1:xl, 1:xx]);
    elseif (ismatrix (x))
      nd = ndims (x);
      sz = size (x);
      sz2 = ceil (sz ./ 2);
      idx = cell ();
      for i = 1:nd
        idx{i} = [sz2(i)+1:sz(i), 1:sz2(i)];
      endfor
      retval = x(idx{:});
    else
      error ("fftshift: expecting vector or matrix argument");
    endif
  endif

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

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


