% WEBDEMO coherent_optical_transmission_advanced.m-file, Institute of Telecommunications, University of Stuttgart, Germany
% www.inue.uni-stuttgart.de
% author: Michael Bernhard, michibernhard@gmx.de
% date: 15.01.2014
%
% Written for GNU Octave 
%
% EXAMPLE:  
%
%
% PARAMETERS:
%  fiber_length,    % in km
%  symbol_rate,     % in GHz  
%  chromatic_dispersion, % in ps/km/nm
%  impulse_shaper,  % 6: rectangular; 4: raised cosine, alpha = 0.5; 5: raised cosine, alpha = 1  
%  pulse_mode,      % nrz, rz
%  modulation_order,% qam_4, qam_16, qam_64, ook or bpsk
%  eye_diagram_en,  % if 0: off (means constellation diagram), 1: on (means eye diagram)
%  jitter_en,       % if 0: no jitter, 1: jitter
%  equalizer_en,    % if 0: no equalizer, 1: equalizer on
%  snr,             % SNR: additional noise
%  oversampling_equalizer, % if 0: 1 Sample/Symbol, 1: oversampling mode for eqalizer
%  max_number_of_taps, % maximum number of taps for equalizer
%  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 coherent_optical_transmission_time_domain_eq(
  fiber_length,     % in km
  symbol_rate,      % in GHz  
  chromatic_dispersion, % in ps/km/nm
  impulse_shaper,   % 6: rectangular;  4: raised cosine, alpha = 0.5; 5: raised cosine, alpha = 1  
  pulse_mode,       % nrz, rz
  modulation_order, % qam_4, qam_16, qam_64, ook or bpsk
  eye_diagram_en,   % if 0: off (means constellation diagram), 1: on (means eye diagram)
  jitter_en,        % if 0: no jitter, 1: jitter
  equalizer_en,     % if 0: no equalizer, 1: equalizer on 
  snr,              % SNR: additional noise
  oversampling_equalizer, % if 0: 1 Sample/Symbol, 1: oversampling mode for eqalizer
  max_number_of_taps, % maximum number of taps for equalizer
  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
  %snr = 15;
  inset_en = 0;
  if eye_diagram_en == 1;
    inset_en = 1;
  end

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


  set(0, "defaultlinelinewidth", 2);

  intensity_plot = 0;

  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

  hold off;
  % 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/km/nm in s/m/m

 
  % END check input parameters for consistency
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  hold off;
  t1 = 1/symbol_rate;
  M = 200;

  tx_data = zeros(1,M);
  scaling = 1;
  for k=1:M
    if strcmp(modulation_order, 'bpsk')
      scaling = 1;
      tx_data(1,k) = (round(rand)*2 - 1) / scaling;      
    else
      if strcmp(modulation_order, 'ook')
        scaling = 1/sqrt(2);
        tx_data(1,k) = (round(rand)) / scaling;
      else
        if strcmp(modulation_order, 'qam_4')
          scaling = sqrt(2);
          tx_data(1,k) = ((round(rand) + 1j* round(rand)) * 2 - 1 - 1j) / scaling;
        else
          if strcmp(modulation_order, 'qam_16')
            scaling = 3 / sqrt(10);
            tx_data(1,k) = ((round(rand*3) / 3 + 1j* round(rand*3) / 3) * 2 - 1 - 1j) / scaling; 
          else % qam 64
            scaling = 7 / sqrt(42);
            tx_data(1,k) = ((round(rand*7) / 7 + 1j* round(rand*7) / 7) * 2 - 1 - 1j) / scaling;
          end 
        end
      end
    end

  end

  oversampling_factor = 64;
  if impulse_shaper == 6
    oversampling_factor = 256;
  end
 
 
  nb_points = oversampling_factor * size(tx_data,2); % total number of points for plot
 
  % create tx signal
  [tx_signal_time tx_signal_data] = create_tx_signal(tx_data, impulse_shaper, pulse_mode, t1, nb_points);
  
  Es = 1;%mean(abs(tx_data).*abs(tx_data));

  % compute rx signal
  [rx_signal_time rx_signal_data] = compute_rx_signal(tx_signal_time, tx_signal_data, chromatic_dispersion, fiber_length, symbol_rate, oversampling_factor, snr, Es);
 
  % helper variable
  step = size(rx_signal_data,2)/M;
    
  % optical - electrical conversion
  rx_electrical_signal_data = [];
  if intensity_plot == 1
    rx_electrical_signal_data = abs(rx_signal_data).^2;
    disp('intensity')
  else
    rx_electrical_signal_data = rx_signal_data;
  end
 
  % downsample electical data (ADC):
  rx_data_sampled = dac(rx_electrical_signal_data, jitter_en, step, pulse_mode, 0);

  equalizer_before_adc = oversampling_equalizer;

  % equalizer
  offset = 0;
  number_of_used_taps= 0; %(time_shift))
  necessary_taps = 0;
  step = 1;
  if equalizer_before_adc == 1
    [eq_data number_of_used_taps necessary_taps] = equalizer_time_domain(rx_electrical_signal_data, chromatic_dispersion, fiber_length, symbol_rate, oversampling_factor, max_number_of_taps);

    %% eq_data_sampled = dac(eq_data, jitter_en, step, pulse_mode, offset);
    step = size(rx_signal_data,2)/M;
  else
    [eq_data_sampled number_of_used_taps necessary_taps] = equalizer_time_domain(rx_data_sampled, chromatic_dispersion, fiber_length, symbol_rate, 1, max_number_of_taps);
    eq_data = rx_electrical_signal_data;%eq_data_sampled;
    step = size(rx_signal_data,2)/M;
  end

  %equalizer
  if equalizer_en == 1
    index = 0;
  if number_of_used_taps > 0 
    index = step + ceil(number_of_used_taps/2)+1;%mod(number_of_used_taps, step);%number_of_used_taps; %ceil(number_of_used_taps / step) * step - number_of_used_taps;   
    %index = ceil((number_of_used_taps) / step) * step;% - number_of_used_taps;
  end
  
  rx_signal_data = eq_data(1,index+1:end);
  
  %size of eq_data is smaller than the origianl data size: calculate M:
  if equalizer_before_adc == 1
    M = floor(size(rx_signal_data,2) / step);
  else
    M = floor(size(rx_signal_data,2) / step); 
  end 
 end

 
  % optical - electrical conversion
  rx_electrical_signal_data = [];
  if intensity_plot == 1
    rx_electrical_signal_data = abs(rx_signal_data).^2;
  else
    rx_electrical_signal_data = rx_signal_data;
  end
 

  xlim_a =rx_signal_time(1+step/2)*1e9;
  xlim_b =rx_signal_time(step*2+step/2+1)*1e9;
 

  inset_x = [];
  inset_y = [];
  rx_data_all = [];
  
  for m=3:M-3
    a = 1+step*(m-1)+step/2;
    b = a+step*2-1; 
    index = 1 + step*m + step/2;
    index_base = 1 + step*1 + step/2;
    if strcmp(pulse_mode, 'nrz')
      a = 1+step*(m-1);
      b = a+step*2-1;
      index = 1 + step*m;
    end
    
    % sample signal:
    % jitter:
    jitter =  0;
    if jitter_en == 1
      % timing jitter Delta tau = 0.3 T:
      jitter = round(rand*step*0.3 - step*0.15)
    end
    rx_data = rx_electrical_signal_data(index+jitter);
    rx_data_all = [rx_data_all rx_data];
    
    if eye_diagram_en == 0
      plot(real(rx_data), imag(rx_data), 'b.');
      if equalizer_en == 1 && equalizer_before_adc == 0
          hold on;
          plot(real(eq_data_sampled), imag(eq_data_sampled), 'r.')
      end
    else
      subplot(2,1,1);
      plot(rx_signal_time(1+step/2:step*2+step/2)*1e9, real(rx_electrical_signal_data(a:b)), "linewidth", const_lw, "color", "m");
      hold on;
      plot(rx_signal_time(index_base+jitter)*1e9, real(rx_data), 'bo');
      subplot(2,1,2);
      plot(rx_signal_time(1+step/2:step*2+step/2)*1e9, imag(rx_electrical_signal_data(a:b)), "linewidth", const_lw, "color", "m");
      hold on;
      plot(rx_signal_time(index_base+jitter)*1e9, imag(rx_data), 'bo');
      
      if inset_en
        inset_x = [inset_x real(rx_data)];
        inset_y = [inset_y imag(rx_data)];
             
      end
    end
    
    hold on;
  end
 
  if eye_diagram_en == 0
    axis equal
    if strcmp(modulation_order, 'ook')
      xlim([-0.5 1.7]);
      ylim([-1 1]);
    else
      xlim([-1.2 1.2]);
      ylim([-1.2 1.2]);
    end
    xlabel('I');
    ylabel('Q');
    if equalizer_en == 1
      title(['Nr. of used taps: ' num2str(number_of_used_taps) ' of ' num2str(necessary_taps)]);
    end
  else
    if strcmp(modulation_order, 'ook')
      ylim_a = -1.4;
      ylim_b = -1,4
    else
      ylim_a = -1.4;
      ylim_b = 1.4;    
    end
    subplot(2,1,1);
    xlabel('time t/ns');
    ylabel('I');
    if equalizer_en == 1
      title(['Nr. of used taps: ' num2str(number_of_used_taps) ' of ' num2str(necessary_taps)]);
    end
    subplot(2,1,2);  
    xlabel('time t/ns');
    ylabel('Q');
  end
    
  if inset_en       
    inset_pos = [0.75 0.65 0.2 0.3];
    sub_ax = axes('position', inset_pos);        
    plot(inset_x, inset_y, 'b.');
    if equalizer_en == 1 && equalizer_before_adc == 0
      hold on;
      plot(real(eq_data_sampled), imag(eq_data_sampled), 'r.')
    end
    axis equal
    set(sub_ax, 'color', 'w');        
    set(sub_ax, 'XTick', [-1 0 1]);
    set(sub_ax, 'YTick', [-1 0 1]);             
  end
  

  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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DAC
function rx_data = dac(rx_signal_data, jitter_en, step, pulse_mode, offset)

  rx_data = [];

  M = floor(size(rx_signal_data,2) / step);
  
  for m=3:M-3
    index = 1 + step*m + step/2 + offset;
    if strcmp(pulse_mode, 'nrz')
      index = 1 + step*m + offset;
    end
    
    % sample signal:
    % jitter:
    jitter =  0;
    if jitter_en == 1
      % timing jitter Delta tau = 0.3 T:
      jitter = round(rand*step*0.3 - step*0.15)
    end
    rx_data = [rx_data rx_signal_data(index+jitter)];
  end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


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

  N1 = N;
  k_step = 1;
  if strcmp(pulse_mode, 'nrz') % NRZ
    N1 = N/2;
    k_step = 1;  
  end
  tx_symbol_1 = zeros(1,nb_points/N1);
  tx_symbol_2 = zeros(1,nb_points/N1);
  tx_symbol = zeros(1,nb_points/N1);

  for k=1:k_step:N % symbol index
    for time=0:nb_points/N1-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/N1/2)/(nb_points/N1/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/N1/2)/(nb_points/N1/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;
      case 6
	t_n = (time - nb_points/N1/2)/(nb_points/N1/2);
        alpha = 0.0000000000001;
	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
  if strcmp(pulse_mode, 'nrz') % NRZ
    a = nb_points/N; %nb_points/N1*3/4
    tx_signal_data(1,1+(k-1)*nb_points/N:k*nb_points/N) = [tx_symbol_2(1,a+1:end) + tx_symbol(1,1:a)];   
    tx_symbol_1 = tx_symbol_2;
    tx_symbol_2 = tx_symbol;
  else
    tx_signal_data(1,1+(k-1)*nb_points/N:k*nb_points/N) = tx_symbol(1,:);
  end
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

if strcmp(pulse_mode, 'nrz') % NRZ
  tx_signal_data = tx_signal_data;
end

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




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute RX signal
function [rx_signal_time rx_signal_data] = compute_rx_signal(tx_signal_time, tx_signal_data, d_c, fiber_length, symbol_rate, oversampling_factor, snr, Es);
  N=size(tx_signal_data,2);
  rx_signal_data = zeros(1, N);
  rx_signal_time = tx_signal_time;


  c = 299792458; % speed of light in m/s
  lm0 = 1550e-9; % reference wavelength in m
  D = d_c;%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( -1j * pi * D * lm0^2 / c * L * f.^2);
  H = my_fftshift(H);

  freq = fft(tx_signal_data);
  freq = freq .* H;
  rx_signal_data = ifft(freq);
  N0 =  Es * 10^(-snr/10);
  for k=1:N  
    rx_signal_data(1,k) = rx_signal_data(1,k) + (rand()*2-1)*sqrt(N0/2) + 1j*(rand()*2-1)*sqrt(N0/2);
  end
endfunction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute time response
function time_response_data = compute_time_response(filter_length, d_c, fiber_length, symbol_rate, oversampling_factor)
  snr = 1000;
  tx_time = 0; % not used
  tx_data = zeros(1, filter_length);
  tx_data(1,1) = 1;
  Es = mean(abs(tx_data).*abs(tx_data));
  [time_response_time, time_response_data] = compute_rx_signal(tx_time, tx_data, d_c, fiber_length, symbol_rate, oversampling_factor, snr, Es);
endfunction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Frequency domain equalizer
function [equalized_data] = equalizer(input_data, d_c, fiber_length, symbol_rate, oversampling_factor)

  N=size(input_data,2);
  equalized_data = zeros(1, N);

  c = 299792458; % speed of light in m/s
  lm0 = 1550e-9; % reference wavelength in m
  D = d_c;%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(- 1j * pi * D * lm0^2 / c * L * f.^2);
  H = my_fftshift(H);
  freq = fft(input_data);
  freq = freq .* H;
  equalized_data = ifft(freq);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [equalized_data time_shift necessary_taps] = equalizer_time_domain(input_data, d_c, fiber_length, symbol_rate, oversampling_factor, max_number_of_taps)
c = 299792458; % speed of light in m/s
lm0 = 1550e-9; % reference wavelength in m
D = d_c;%17e-6; % chromatic dispersion coefficient in s/m^2 at lm0
L = fiber_length; %100e3; % fiber length in m
T = 1 / symbol_rate / oversampling_factor ;

N = 2* floor(abs(D)*lm0*lm0*L/(2*c*T*T)) + 1;
necessary_taps = N;
if N > max_number_of_taps
 N = max_number_of_taps
end


disp(['N: ' num2str(N)]);
a = zeros(1,N);


for k = 1:N
    m = k-floor(N/2)-1;
    a(1,k) = sqrt(1j*c*T*T/(D*lm0*lm0*L))*exp(-1j*pi*c*T*T/(D*lm0*lm0*L)*m*m);
end

equalized_data = filter(a, 1, input_data);
%equalized_data = equalized_data(N:end);
time_shift = N;   

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




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ds_pos = fig2ds_pos(ah, pos)
% FIG2DS_POS Convert a position from normalized figure to data space units.
%
% IN:
%   ah  Handle of axes that contains the dataspace
%   pos Four element position vector [left bottom width height] in normalized
%       figure units. (Values between 0 and 1.0).
% OUT:
%  ds_pos Four element position vector in data space units of the given axes.

    % Position that is to be converted
    x = pos(1);
    y = pos(2);
    width = pos(3);
    height = pos(4);
    
    % Limits of data space of given axes
    ax_xlim = get(ah, 'XLim');
    ax_ylim = get(ah, 'YLim');
    dsx_max = ax_xlim(2);
    dsx_min = ax_xlim(1);
    dsy_max = ax_ylim(2);
    dsy_min = ax_ylim(1);
    
    % Position of given axes inside the figure
    ax_pos = get(ah, 'Position');
    fig_x = ax_pos(1);
    fig_y = ax_pos(2);
    fig_width = ax_pos(3);
    fig_height = ax_pos(4);
    
    % Conversion
    ds_pos(1) = (dsx_max - dsx_min)/fig_width * (x - fig_x) + dsx_min;
    ds_pos(2) = (dsy_max - dsy_min)/fig_height * (y - fig_y) + dsy_min;
    ds_pos(3) = (dsx_max - dsx_min)/fig_width * width;
    ds_pos(4) = (dsy_max - dsy_min)/fig_height * height;
end


