function demo_pathloss( 
	env_sel,		% environment selection 
	fc,			% 2.4 or 5.0 GHz 
	d_min,			% min distance in m 
	d_max,			% max distance in m 
	nb_points,		% nb of points in plot 
	plot_sel,		% plot selection ((1)none, (2)rand_dist, (3)pdm)
	filename		% output file for webdemo 
)
% create figure and set style 
fig = figure(1); 
set(fig, "visible", "on"); 
landscape = "-S930,350"; 
portrait  = "-S640,480"; 
output_format = portrait; 
set(0, "defaultlinelinewidth", 2); 
set(0, "defaultaxesfontsize", 10); 
set(0, "defaultaxesfontname", "Helvetica"); 
set(0, "defaulttextfontsize", 8); 
set(0, "defaulttextfontname", "Helvetica"); 
const_lw = 2; 
const_msize = 2; 

% set (environment) parameters 
d_step = (d_max-d_min)/nb_points; 
d = (d_min:d_step:d_max);

[d_bp, sigma_S, tap_delay, tap_power, cl_AoA_rAS_AoD_tAS, cl_power, cl_9, cl_1113, cl_15, cl_1012, cl_13, cl_14, rtp] = gep(env_sel);

% calc y-vectors 
y_bpl = calc_y_bpl(fc, d_bp, d_min, d_max, nb_points); 
y_s = calc_y_s(d_bp, sigma_S, d_min, d_max, nb_points); 
y_bpl_s  = y_s + y_bpl; 

% plot 
hold on
xlabel("distance in m"); 
ylabel("path loss in dB"); 
xlim([d_min, d_max]); 
if (plot_sel == 2)
	semilogx(d, y_bpl, "linewidth", const_lw);
	semilogx(d, y_bpl_s, "+","linewidth", const_lw,"color","r", "markersize", const_msize);
	legend("path loss", "shadow fading", "location", "northeast");
elseif (plot_sel == 3)
	clf
	hold on
	xlabel("distance in m"); 
	ylabel("path loss in dB"); 
	set(0, "defaultaxesfontsize", 6); 
	colormap(jet);
	imagesc([d_min d_max], [min(y_bpl) max(y_bpl)], calc_pd_matrix(y_bpl, d_min, d_max, d_bp, sigma_S, nb_points));
	plot(d, y_bpl, "linewidth", const_lw, "color", "m");   
	h = colorbar;
	ylabel(h, 'probability density of shadow fading')
else
	semilogx(d, y_bpl, "linewidth", const_lw, "color", "m");
	legend("path loss", "location", "northeast");	
endif
 
% Create an image file 
print(filename, "-dpng",output_format); 
endfunction

% Other functions
function y_bpl_vec = calc_y_bpl(fc, d_bp, d_min, d_max, nb_points) 
	% y_bpl = calc_y_bpl(fc, d_bp, d_min, d_max, nb_points); 
	lambda = 3*10^8/(fc*10^9); 
	d_step = (d_max-d_min)/nb_points; 
	d = d_min:d_step:d_max; 
	y_bpl_vec = zeros(1,length(d)); 
	k=1; 
	for x = d 
		if (x <= d_bp) 
			y_bpl_vec(k) = 10*log10((lambda/4/pi/x)^2); 
		else 
			y_bpl_vec(k) = 10*log10((lambda/4/pi/d_bp)^2) - 35*log10(x/d_bp); 
		endif 
		k = k + 1; 
	endfor 
endfunction

function y_s_vec = calc_y_s(d_bp, sigma_S, d_min, d_max, nb_points) 
	% y_s = calc_y_s(d_bp, sigma_S, d_min, d_max, nb_points); 
	d_step = (d_max-d_min)/nb_points; 
	d = d_min:d_step:d_max;     
	y_s_vec = zeros(1,length(d)); 
	k = 1; 
	for x = d 
		if(x <= d_bp) 
			y_s_vec(k) = randn*3; 
		else 
			y_s_vec(k) = randn*sigma_S; 
		endif 
		k = k + 1; 
	endfor 
endfunction

function pd_matrix = calc_pd_matrix(y_bpl, d_min, d_max, d_bp, sigma_S, nb_points) 
	% pdm = calc_pd_matrix(y_bpl, d_min, d_max, d_bp, sigma_S, nb_points); 
	x_step = (d_max - d_min) / nb_points; 
	y_step = abs(max(y_bpl) - min(y_bpl)) / nb_points;
	pd_matrix = zeros (nb_points);
	x_in = 1; 
	for x = d_min:x_step:d_max;
		if (x <= d_bp) 
			sigma = 3; 
		else 
			sigma = sigma_S; 
		endif 
		mu = y_bpl(x_in); 
		y_in = 1; 
		for y = min(y_bpl):y_step:max(y_bpl)
			pd_matrix(y_in, x_in) = (1/sqrt(2*pi*sigma^2))*exp(-(y - mu)^2/(2*sigma^2)); 
			y_in = y_in + 1; 
		endfor 
		x_in = x_in + 1; 
	endfor;     
endfunction
