% WEBDEMO .m-file (Octave), Institute of Telecommunications, University of Stuttgart, Germany
% AUTHOR:   Stephan ten Brink
% EXAMPLE:  path_loss(10,10000, 900, 30, 1, 0, 0.44, 4, 1, 200, 2000, 6, 1,1, [1 0 0 0], 0, 0, "path_loss")
% PARAMETERS:
%   dmin,                % in m
%   dmax,                % in m
%   fc,                  % in MHz
%   h_b,                 % in m, antenna height, base station
%   h_m,                 % in m, antenna height, mobile terminal
%   OH_scenario,         % NOT USED; Okumura-Hata model, scenario (0=small/medium-sized cities, 1=metropolitan, 2=suburban, 3=rural)
%   MK_alpha,            % MotleyKeenan indoor model: in dB/m
%   BP_exp,              % Breakpoint model, path loss exponent
%   logx_en,             % log on x axis (effectively double-log)
%   nb_points,           % nb of points in plot
%                        %
%   nb_points_measured,  % for emulating measured data; nb of points
%   sigma_S,             % standard deviation, shadow fading
%   shadow_en,           % shadow fading enable
%   rayleigh_en,         % Rayleigh fading enable
%                        %
%   plot_en,             % enable plots 0,1 [free space, Okumura-Hata, MotleyKeenan indoor, Breakpoint model]
%   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


% further examples
% - PL: path_loss(10,10000, 900, 30, 1, 0, 0.44, 4, 1, 200, 2000, 6, 1,1, [1 0 0 0], 0, 0, "path_loss")
% - OH: path_loss(1000,10000, 900, 30, 1, 0, 0.44, 4, 0, 200, 0,0,0,0, [0 1 0 0], 0, 0, "path_loss")
% - OH: path_loss(100,10000, 900, 30, 0.5, 0, 0.44, 4, 0, 200, 0,0,0,0, [0 1 0 0], 0, 0, "path_loss")
% - MK: path_loss(1,100, 2400, 30, 1, 0, 0.44, 4, 1, 200, 0,0,0,0, [0 0 1 0], 0, 0, "path_loss")
% - BP: path_loss(1,10000, 900, 10, 1.5, 0, 0.44, 4, 1, 200, 0,0,0,0, [0 0 0 1], 1, 0, "path_loss")
% - BP: path_loss(20,20000, 900, 10, 1.5, 0, 0.44, 4, 1, 20000, 0,0,0,0, [0 0 0 1], 0, 0, "path_loss_BP")

% notes, can be included in future
% - breakpoint model; CROSS-CHECK OTHER REFERENCES (normalization)
% - OH, graphically mask for d<1km
% - show minor grid lines, when logx enabled

% extensions (separate demo)
% - include link-budget considerations
% - superimpose shadow fading; also plot 2dim view (for fun)

% BEGIN, main function (control)

function path_loss(
  dmin,                % in m
  dmax,                % in m
  fc,                  % in MHz
  h_b,                 % in m, antenna height, base station
  h_m,                 % in m, antenna height, mobile terminal
  OH_scenario,         % NOT USED; Okumura-Hata model, scenario (0=small/medium-sized cities, 1=metropolitan, 2=suburban, 3=rural)
  MK_alpha,            % MotleyKeenan indoor model: in dB/m
  BP_exp,              % Breakpoint model, path loss exponent
  logx_en,             % log on x axis (effectively double-log)
  nb_points,           % nb of points in plot
                       %
  nb_points_measured,  % for emulating measured data; nb of points
  sigma_S,             % standard deviation, shadow fading
  shadow_en,           % shadow fading enable
  rayleigh_en,         % Rayleigh fading enable
                       %
  plot_en,             % enable plots 0,1 [free space, Okumura-Hata, MotleyKeenan indoor, Breakpoint model]
  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;     % 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;
    const_msize = 4;
  elseif (svg_en)
    set(0, "defaultaxesfontsize", 6);
    set(0, "defaultaxesfontname", "Helvetica");
    set(0, "defaulttextfontsize", 6);
    set(0, "defaulttextfontname", "Helvetica");
    const_lw = 1;
    const_msize = 3;
  else
    set(0, "defaultaxesfontsize", 10);
    set(0, "defaultaxesfontname", "Helvetica");
    set(0, "defaulttextfontsize", 8);
    set(0, "defaulttextfontname", "Helvetica");
    const_lw = 2;
    const_msize = 2;
  endif


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



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

  % only one plot at a time

  if (sum(plot_en)!=1)
    plot_en = zeros(1,100); plot_en(1) = 1;
    warning("setting  'plot_en = [1 0 ... 0]'");
  endif

  if (nb_points>40000)
    nb_points = 40000;
    warning("setting  'nb_points = 40000'");
  end

  if (nb_points_measured>10000)
    nb_points_measured = 10000;
    warning("setting  'nb_points_measured = 10000'");
  end

  % END
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


  d_step = (dmax-dmin)/nb_points;

  d = (dmin:d_step:dmax);

  y_free = compute_PL_free_dB(fc,d);

  clf % other resets/clear commands needed?
  hold on


  xlabel("distance in m");
  ylabel("path loss in dB");
  xlim([dmin, dmax]);


  % free-space path loss
  if (plot_en(1))

    % models measurements according to lognormal/shadow fading and small-scale fading
    N_rnd = nb_points_measured;
    lambda = 3e8/(fc*1e6);
    x_rnd = rand(1,N_rnd)*(dmax-dmin)+dmin;
    if (shadow_en)
      y_rnd_S  = randn(1,N_rnd)*sigma_S;
    else
      y_rnd_S  = 0;    
    end
    y_PL_rnd_S  = y_rnd_S + 10*log10((lambda/4/pi./x_rnd).^2);
    %
    if (rayleigh_en)
      y_PL_rnd_S_R2 = y_PL_rnd_S + 10*log10((randn(1,N_rnd)/sqrt(2)).^2+(randn(1,N_rnd)/sqrt(2)).^2);
    end

    if (logx_en)
      if (shadow_en)
        semilogx(x_rnd, y_PL_rnd_S, "+","linewidth", const_lw,"color","r", "markersize", const_msize);
      end
      if (rayleigh_en)
        semilogx(x_rnd, y_PL_rnd_S_R2, ".","linewidth",const_lw,"color","g");
      end
      semilogx(d, y_free, "linewidth", const_lw, "color", "m");
    else % lin. scale
      if (shadow_en)
        plot(x_rnd, y_PL_rnd_S, "+","linewidth", const_lw,"color","r", "markersize", const_msize);
      end
      if (rayleigh_en)
        plot(x_rnd, y_PL_rnd_S_R2, ".","linewidth",const_lw,"color","g");
      end
      plot(d, y_free, "linewidth",  const_lw, "color", "m");
    endif
    
    if ((shadow_en) & (rayleigh_en)) 
      legend("Shadow fading","Rayleigh fading", "path loss", "location", "northeast");
    elseif (shadow_en) 
      legend("Shadow fading", "path loss", "location", "northeast");
    elseif (rayleigh_en) 
      legend("Rayleigh fading", "path loss", "location", "northeast");
    else
      legend("path loss", "location", "northeast");
    endif
      
    output_format = portrait;
    
  % Okumura-Hata model
  elseif (plot_en(2))

    if (logx_en)
      semilogx(d, y_free, "linewidth", const_lw, "color", "m");
      y_OH = compute_PL_OH_dB(fc,d, h_b, h_m, 3);
      semilogx(d, y_OH, "linewidth", const_lw, "color", "r");
      y_OH = compute_PL_OH_dB(fc,d, h_b, h_m, 2);
      semilogx(d, y_OH, "linewidth", const_lw, "color", "g");
      y_OH = compute_PL_OH_dB(fc,d, h_b, h_m, 1);
      semilogx(d, y_OH, "linewidth", const_lw , "color", "b");
      y_OH = compute_PL_OH_dB(fc,d, h_b, h_m, 0);
      semilogx(d, y_OH, "linewidth", const_lw, "color", "c");
    else
      plot(d, y_free, "linewidth", const_lw, "color", "m");
      y_OH = compute_PL_OH_dB(fc,d, h_b, h_m, 3);
      plot(d, y_OH, "linewidth", const_lw, "color", "r");
      y_OH = compute_PL_OH_dB(fc,d, h_b, h_m, 2);
      plot(d, y_OH, "linewidth", const_lw, "color", "g");
      y_OH = compute_PL_OH_dB(fc,d, h_b, h_m, 1);
      plot(d, y_OH, "linewidth", const_lw, "color", "b");
      y_OH = compute_PL_OH_dB(fc,d, h_b, h_m, 0);
      plot(d, y_OH, "linewidth", const_lw, "color", "c");
    endif

    legend("free-space","rural","suburban", "metropolitan", "small/medium-sized cities",  "location","northeast");
    output_format = portrait;

  % Motely Keenan model
  elseif (plot_en(3))
    y_MK = compute_PL_MK_dB(fc,d, MK_alpha);
    if (logx_en)
      semilogx(d, y_free, "linewidth", const_lw, "color", "m");
      semilogx(d, y_MK, "linewidth", const_lw, "color", "r");
    else
      plot(d, y_free, "linewidth", const_lw, "color", "m");
      plot(d, y_MK, "linewidth", const_lw, "color", "r");
    endif


    legend("free-space","Motley-Keenan model");
    output_format = portrait;

  % Breakpoint model
  elseif (plot_en(4))        % breakpoint (two-path) model, d^{-4}-law
    lambda = 3e8/(fc*1e6);
    d_break = 4*h_b*h_m/lambda;

%text(2e3, -40, ["d_{break}=", mat2str(d_break),"m"], "fontsize", 12, "color", "r");
%ylim([-140 0]);

    y_BP = compute_PL_BP_dB(fc,d, h_b, h_m, BP_exp);
%    y_BP_s = compute_PL_BP_simple_dB(fc, d, h_b, h_m, BP_exp);
    if (logx_en)

      % also plot ref. curve for d^{-4}
      y_pow4=zeros(1,length(d));
      k=1;
      for x=d
        y_pow4(k)=10*log10((h_b*h_m/x^2)^2);
        k=k+1;
      endfor
      %

      % for y-axis scaling: compute minimum value
      pl_min = min(y_pow4(end), y_free(end));
      ylim([1.1*pl_min 0]);

      semilogx(d, y_pow4, "linewidth", const_lw, "color", "b");

      semilogx(d, y_free, "linewidth", const_lw, "color", "m");
      semilogx(d, y_BP, "linewidth", const_lw, "color", "r");

    else
      plot(d, y_free, "linewidth", const_lw, "color", "m");
      plot(d, y_BP, "linewidth", const_lw, "color", "r");
    endif

    legend("d^{-4} reference", "free-space", ["Breakpoint model ", "(d_{break}=", num2str(d_break),"m)"]);
    output_format = portrait;

  endif

  hold off


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


endfunction

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






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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y_vec = compute_PL_free_dB(fc,d_vec)
  lambda = 3e8/(fc*1e6);
  len = length(d_vec);
  y_vec = zeros(1,len);

  k=1;
  for x=d_vec
    y_vec(k) = 10*log10((lambda/4/pi/x)^2);
    k=k+1;
  endfor
endfunction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Breakpoint model
function y_vec = compute_PL_BP_simple_dB(fc, d_vec, h_b, h_m, BP_exp)
  lambda = 3e8/(fc*1e6);
  len = length(d_vec);
  y_vec = zeros(1,len);

  % convention:
  P_T_at1m_dB = 10*log10((lambda/4/pi/1)^2);

  d_break = 4*h_b*h_m/lambda;
  BP_exp=2;    % check path loss exponent after breakpoint; ok

  k=1;
  for x=d_vec
  if (x<d_break)
    y_vec(k) = 10*log10((lambda/4/pi/x)^2);
    k=k+1;
  else
    y_vec(k) = P_T_at1m_dB - 2*10*log10(d_break) - BP_exp*10*log10(x/d_break);
    k=k+1;
  endif
  endfor
endfunction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Breakpoint model 2

function y_vec = compute_PL_BP_dB(fc, d_vec, h_b, h_m, BP_exp)
  lambda = 3e8/(fc*1e6);
  len = length(d_vec);
  y_vec = zeros(1,len);


  k=1;
  for x=d_vec
    y_vec(k) = 10*log10((lambda/4/pi/x)^2) + 10*log10(4* (sin(pi*2*h_b*h_m/x/lambda))^2);
    k=k+1;
  endfor
endfunction

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



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Motley-Keenan model
function y_vec = compute_PL_MK_dB(fc, d_vec, MK_alpha)
  lambda = 3e8/(fc*1e6);
  len = length(d_vec);
  y_vec = zeros(1,len);

  k=1;
  for x=d_vec
    y_vec(k) = 10*log10((lambda/4/pi/x)^2) - MK_alpha*x;
    k=k+1;
  endfor
endfunction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Okumura-Hata model
function y_vec = compute_PL_OH_dB(fc,d_vec, h_b, h_m, OH_scenario);


  % check input parameters for consistency
  if ((h_b<30) || (h_b>200))
    h_b=30;            % reset to default
    warning("setting 'h_b=30'");
  endif

  if ((h_m<1) || (h_m>10))
    h_m=3;            % reset to default
    warning("setting 'h_m=3'");
  endif

  if ((fc<150) || (fc>1500))   % note: 1500..2000 not yet included
    fc=900;            % reset to default
    warning("setting 'fc=900'");
  endif





  if (fc<1500)
    if ((OH_scenario==0) || (OH_scenario==2) || (OH_scenario==3))
      a = (1.1*log10(fc)-0.7)*h_m - (1.56*log10(fc)-0.8);
    endif

    if (OH_scenario==0)      %  small and medium-sized cities
      C = 0;
    elseif (OH_scenario==1)  %  metropolitan
      if (fc<300)
        a = 8.29*(log10(1.54*h_m))^2 - 1.1;
      else % fc>=300
        a = 3.2*(log10(11.76*h_m))^2 - 4.97;
      endif
      C = 0;
    elseif (OH_scenario==2)  %  suburban
      C = -2*(log10(fc/28))^2 - 5.4;
    else                     %  rural
      C = -4.78*(log10(fc))^2 + 18.33*log10(fc) - 40.94;
    endif

    A = 69.55+26.16*log10(fc)-13.82*log10(h_b) - a;
    B = 44.9-6.55*log10(h_b);

  else    % 1500MHz<fc<2000MHz
    a = (1.1*log10(fc)-0.7)*h_m - (1.56*log10(fc)-0.8);
    A = 46.3+33.9*log10(fc) - 13.82*log10(h_b) - a;
    B = 44.9-6.55*log10(h_b); % unchanged
    if (OH_scenario==1)
      C=3;
    else
      C=0;
    endif

  endif

  len = length(d_vec);
  y_vec = zeros(1,len);

  k=1;
  for x=d_vec
    y_vec(k) = -(A + B*log10(x/1e3) + C);
    k=k+1;
  endfor

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



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


