function demo_multipath( 
	env_sel,			% environment selection 
  	plot_t_sel,			% plot type selection
	domain,				% 1 = time domain 2 = frequency domain
	channel_sel,				% enable 3 Dimensional View
	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;

% channel & environment parameters
[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);
tap_delay = tap_delay./10^-9;
[f_min, f_max, fc] = gcp(channel_sel);
rgbcolormat= [0 0 1;
	      0 1 0;
              0 1 1;
	      1 0 0;
	      1 0 1;
	      0 0 0;];
stringcolormat=['k', 'm', 'c', 'r', 'g', 'b'];
[cz cs] = size(cl_power);
% plot PDP
if (domain == 1) % time domain
	if(plot_t_sel == 1)
		if (tap_power(1) == 0)
			tap_power(tap_power == 0) = NaN; % not drawing circles at 0 values
			cl_power(cl_power == 0) = NaN; % not drawing circles at 0 values
			tap_power(1) = 0;
			cl_power(1,1) = 0;
		else
			tap_power(tap_power == 0) = NaN; % not drawing circles at 0 values
			cl_power(cl_power == 0) = NaN; % not drawing circles at 0 values
		endif
		set(0, "defaultaxesfontsize", 8);
		subplot(2,1,1)
		stem(tap_delay, calc_tap_power_ref(tap_power, rtp), 'Color', [0 0 1], 'MarkerFaceColor', [0 0 1], 'BaseValue', min(tap_power));
		ylim([min(tap_power) 0]);
		xlim([-10 max(tap_delay)+10]);
		xlabel("distinctive clusters");
		ylabel("relative power in dB");
		title("averaged PDP")
		
		subplot(2,1,2)
		hold on;
		for clindex = cz:-1:1
			stem(tap_delay, calc_tap_power_ref(cl_power(clindex,:), rtp), 'Color', rgbcolormat(clindex,:), 'MarkerFaceColor', rgbcolormat(clindex,:), 'BaseValue', min(tap_power));
		endfor
		ylim([min(tap_power) 0]);
		xlim([-10 max(tap_delay)+10]);
		xlabel("delay in ns");
		ylabel("relative power in dB");


	elseif (plot_t_sel == 2) || (plot_t_sel == 3)
		set(0, "defaultaxesfontsize", 6);
		subplot(2,1,1)
		if(plot_t_sel == 2)
			AoX = cl_AoA_rAS_AoD_tAS(:,1);
			AS = cl_AoA_rAS_AoD_tAS(:,2);
			labelxc = char("Cluster Angle of Arrival in degree");
			labelxt = char("Tap Angle of Arrival in degree"); %labelxt = char("Tap(Marker) + Cluster(Bar) Angle of Arrival in degree"); 
		else
			AoX = cl_AoA_rAS_AoD_tAS(:,3);
			AS = cl_AoA_rAS_AoD_tAS(:,4);
			labelxc = char("Cluster Angle of Departure in degree");
			labelxt = char("Tap Angle of Departure in degree"); %labelxt = char("Tap(Marker) + Cluster(Bar) Angle of Departure in degree"); 
		endif
		cmap = colormap(jet(360));
		caxis([0 360]);
		h = colorbar('NorthOutside');
		xlabel(h, labelxc)
		hold on
		cl_power_ref = zeros(size(cl_power));
		for clindex = cz:-1:1
			cl_power_ref(clindex,:) = calc_tap_power_ref(cl_power(clindex,:),rtp);
		endfor
		if (cl_power_ref(1,1) == 0)
			cl_power_ref(cl_power_ref == 0) = -40; 
			cl_power_ref(1,1) = 0;
		else
			cl_power_ref(cl_power_ref == 0) = -40; % Make 0 dB Taps invisible (Stem draws 0 and NaN crashes Octave for some reason)
		endif
		for clindex = cz:-1:1
			stem(tap_delay, cl_power_ref(clindex,:), 'Color', cmap(round(AoX(clindex)),:) , 'MarkerFaceColor', cmap(round(AoX(clindex)),:), 'LineWidth', 5, 'BaseValue', min(tap_power));
		endfor

		xlim([-10 max(tap_delay)+10]);
		ylim([min(tap_power) 0]);
		xlabel(labelxt);
		ylabel("relative power in dB");
		
		subplot(2,1,2)
		xlim([-10 max(tap_delay)+10]);
		ylim([min(tap_power) 0]);
		xlabel("delay in ns");
		ylabel("relative power in dB");
		legend;
		hold on
		for clindex = cz:-1:1
			for tapindex = cs:-1:1
				single_tap_vec = zeros(1,cs);
				single_tap_vec(tapindex) = cl_power_ref(clindex,tapindex);
				if (single_tap_vec(1) == 0) && (cz == 1)
					single_tap_vec(single_tap_vec == 0) = -40;
					single_tap_vec(1) = 0;
				else
					single_tap_vec(single_tap_vec == 0) = -40; % Make 0 dB Taps invisible (Stem draws 0 and NaN crashes Octave for some reason)
				endif
				
				color = cmap(round(gen_rand_l(AoX(clindex), AS(clindex))),:);
				stem(tap_delay, single_tap_vec, 'Color', color, 'MarkerFaceColor', color, 'LineWidth', 5, 'BaseValue', min(tap_power)); % cmap(round(AoX(clindex)),:)
			endfor
		endfor
	elseif(plot_t_sel == 4) || (plot_t_sel == 5)
		set(0, "defaultaxesfontsize", 6);
		subplot(2,1,1)
		if(plot_t_sel == 4)
			AoX = cl_AoA_rAS_AoD_tAS(:,1);
			AS = cl_AoA_rAS_AoD_tAS(:,2);
			labelxc = char("Tap(Circles) + Cluster(Bar) Angle of Arrival in degree");
			labelxt = char("Angle of Arrival in degree"); %labelxt = char("Tap(Marker) + Cluster(Bar) Angle of Arrival in degree"); 
		else
			AoX = cl_AoA_rAS_AoD_tAS(:,3);
			AS = cl_AoA_rAS_AoD_tAS(:,4);
			labelxc = char("Tap(Circles) + Cluster(Bar) Angle of Departure in degree"); 
			labelxt = char("Angle of Departure in degree"); %labelxt = char("Tap(Marker) + Cluster(Bar) Angle of Departure in degree"); 
		endif
		cmap = colormap(jet(360));
		caxis([0 360]);
		h = colorbar('NorthOutside');
		xlabel(h, labelxc)
		hold on
		cl_power_ref = zeros(size(cl_power));
		for clindex = cz:-1:1
			cl_power_ref(clindex,:) = calc_tap_power_ref(cl_power(clindex,:),rtp);
		endfor
		if (cl_power_ref(1,1) == 0)
			cl_power_ref(cl_power_ref == 0) = -40; 
			cl_power_ref(1,1) = 0;
		else
			cl_power_ref(cl_power_ref == 0) = -40; % Make 0 dB Taps invisible (Stem draws 0 and NaN crashes Octave for some reason)
		endif
		xlim([-10 max(tap_delay)+10]);
		ylim([min(tap_power) 0]);
		xlabel("delay in ns");
		ylabel("relative power in dB");
		hold on
		color = [];
		for clindex = cz:-1:1
			for tapindex = cs:-1:1
				single_tap_vec = zeros(1,cs);
				single_tap_vec(tapindex) = cl_power_ref(clindex,tapindex);
				if (single_tap_vec(1) == 0) && (cz == 1)
					single_tap_vec(single_tap_vec == 0) = -40;
					single_tap_vec(1) = 0;
				else
					single_tap_vec(single_tap_vec == 0) = -40; % Make 0 dB Taps invisible (Stem draws 0 and NaN crashes Octave for some reason)
				endif
				color = cmap(round(gen_rand_l(AoX(clindex), AS(clindex))),:);
				stem(tap_delay, single_tap_vec, 'Color', cmap(round(AoX(clindex)),:), 'MarkerFaceColor', color, 'LineWidth', 5, 'BaseValue', min(tap_power)); % cmap(round(AoX(clindex)),:)
			endfor
		endfor
		
		subplot(2,1,2)
		hold on
		for clindex = 1:cz		
			%subplot(2,cz,clindex+cz)
	%		x = linspace(AoX(clindex)-AS(clindex)/2,AoX(clindex)+AS(clindex)/2,100);
			laplace = calc_laplace(AoX(clindex), AS(clindex));
			x = linspace(AoX(clindex) - 180, AoX(clindex) + 180, 100);
			plot(x,laplace, 'Color', cmap(round(AoX(clindex)),:));
%			if(AoX(clindex) > 180)
%				set(gca,'XTick',[AoX(clindex)-180 AoX(clindex)]);
%				set(gca,'XTickLabel',[AoX(clindex)-180 AoX(clindex)]);
%			else
%				set(gca,'XTick',[AoX(clindex) AoX(clindex)+180]);
%				set(gca,'XTickLabel',[AoX(clindex) AoX(clindex)+180]);
%			endif
%			set(gca, 'YTick', []);
%			%ylim([0 max(laplace)]);
%			xlim([AoX(clindex) - 180 AoX(clindex) + 180]);
%			title(['Cl.' num2str(clindex) ' PDF']);
%			if (clindex == 1)
%				ylabel('probability density');
%				xlabel(labelxt);
%			endif			
		endfor
		xlabel(labelxt);
		ylabel('probability density')
		
	elseif ((plot_t_sel == 6)||(plot_t_sel == 7))
		hold on

		cl_power_ref = zeros(size(cl_power));
		for clindex = 1:cz
			cl_power_ref(clindex,:) = calc_tap_power_ref(cl_power(clindex,:),rtp);
		endfor
		cl_power_ref(cl_power_ref == 0) = NaN;
		[cz cs] = size(cl_power);
		if(plot_t_sel == 6)
			AoX = cl_AoA_rAS_AoD_tAS(:,1);
			AS = cl_AoA_rAS_AoD_tAS(:,2); 
		else
			AoX = cl_AoA_rAS_AoD_tAS(:,3);
			AS = cl_AoA_rAS_AoD_tAS(:,4);
		endif	
		x = zeros(cs,1);
		y = zeros(cs,1);
		z = zeros(cs,1);
		for clindex = 1:cz	
			for tapindex = 1:cs
				tap_angle = gen_rand_l(AoX(clindex),AS(clindex));
				x(tapindex) = tap_delay(tapindex)*cosd(tap_angle);
				y(tapindex) = tap_delay(tapindex)*sind(tap_angle);
				z(tapindex) = cl_power_ref(clindex,tapindex);
			endfor	
			stem3(x,y,z,stringcolormat(clindex));
		endfor
		clmean = linspace(0,max(tap_delay),50);
		for clindex = 1:cz
			for ii = clmean
				stem3(cosd(AoX(clindex))*ii,sind(AoX(clindex))*ii,0,stringcolormat(clindex));
%				stem3(cosd(AoX(clindex)+(AS(clindex)/2))*ii,sind(AoX(clindex)+(AS(clindex)/2))*ii,0,'k'); % Angular Spread
%				stem3(cosd(AoX(clindex)-(AS(clindex)/2))*ii,sind(AoX(clindex)-(AS(clindex)/2))*ii,0,'k'); % Angular Spread
			endfor
		endfor
	view(30, -30);
	grid on;
	set(gca,'FontSize',6);
	xlim([-max(tap_delay) max(tap_delay)]);
	ylim([-max(tap_delay) max(tap_delay)]);
	zlabel('relative Power in dB');
	xlabel('delay in ns');
	ylabel('delay in ns');
	endif
else % |H(f)|
	set(0, "defaultaxesfontsize", 8); 

	subplot(2,1,1);
	stem(tap_delay, calc_tap_power_ref(tap_power, rtp), 'filled', 'BaseValue', min(tap_power));
	xlim([-10 max(tap_delay)+10]);
	ylabel("relative Power in dB");
	xlabel("delay in ns");

	total_tap_power = sum(tap_power);
	tap_power = tap_power ./sqrt(total_tap_power); % normalize tap power so that overall gain = 0db;

	subplot(2,1,2);
	H_f = 0;
	f_step = (f_max - f_min) / 1000;
	f = f_min:f_step:f_max;
	for tap_ind = 1:length(tap_delay)		
		H_f = H_f + tap_power(tap_ind)*exp(-j*2*pi*f*tap_delay(tap_ind)*10^-9);
	endfor
	h_f_mag = 20*log10(abs(H_f));
	plot(f, h_f_mag, 'color', 'red');
	xlim([f_min f_max]);
	set(gca,'XTick',[f_min f_min + (f_max-f_min)/4 f_min + 2*(f_max-f_min)/4 f_min + 3*(f_max-f_min)/4 f_max]);
	set(gca,'XTickLabel',[f_min f_min + (f_max-f_min)/4 f_min + 2*(f_max-f_min)/4 f_min + 3*(f_max-f_min)/4 f_max]/10^6);
	xlabel("frequency in MHz");
	ylabel("|H(f)|^2 in dB");
endif





hold off

% Create an image file 
print(filename, "-dpng",output_format); 
endfunction

% other functions
function tap_power_ref = calc_tap_power_ref(tap_power, rtp)
%	tap_power_ref = zeros(1, length(tap_power));
%	for k = 1:length(tap_power)
%		if (tap_power(k) == 0)
%			tap_power_ref(k) == 0;
%		else
%			tap_power_ref(k) = rtp + tap_power(k);
%		endif
%	endfor
%	if (tap_power_ref(1) == 31)
%		tap_power_ref(1) = rtp;
%	endif
tap_power_ref = tap_power;
endfunction

function randlaplace_value = gen_rand_l(AoX, AS)
		randlaplace_value = (laplace_rnd(1)*AS) + AoX;
		while (randlaplace_value > 360)
			randlaplace_value = randlaplace_value - 360;
		endwhile
		while (randlaplace_value < 0)
			randlaplace_value = randlaplace_value + 360;
		endwhile;
		if (randlaplace_value < 0.5)
			randlaplace_value = 1;
		endif
endfunction

function laplace = calc_laplace(AoX, AS, clindex)
		x = linspace(AoX - 180, AoX + 180, 100);
		laplace = (1./(2.*AS)).*exp(-abs(x-AoX)./AS);
endfunction
