%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% .m script of web demo "Channel Characteristics of Transmission Line" 
% Written for GNU Octave
% 
% INSTITUTE OF TELECOMMUNICATIONS  
% University of Stuttgart 
% www.inue-uni-stuttgart.de 
% author: Jing Lin   
% date: 04.07.2015   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This function creats animationNumberOfPictures frames and prints the values
% of the parameters that have been passed as well as an increasing frame number.
% Make sure that the function prototype matches the prototype defined
% in the XML config file. (Parameter names are arbitrary but should be the same
% for better readability) 

function [] = trans_fun_pdp(logscale, r2, alpha_db, v_ph, l1, l2, l3, f_stop, iter_no, filename)

b = 0;
% overall transfer function
if r2 > 1
	r2 = 1
elseif r2 < -1
	r2 = -1
end
% alpha_db, v_ph, l1, l2, l3, f_stop are larger than zero, invert in case of negative input
if alpha_db < 0
   alpha_db = -alpha_db
   b = 1
end  
if v_ph < 0
   v_ph = -v_ph
   b = 1
elseif v_ph > 3 * 1e5
   v_ph = 3 * 1e5
end
if l1 < 0
   l1 = -l1
   b = 1
end
if l2 < 0
   l2 = -l2
   b = 1
end
if l3 < 0
   l3 = -l3
   b = 1
end
if f_stop < 0
	f_stop = - f_stop
	b = 1
end
% iteration number is from 2 to 12
if iter_no > 12
	iter_no = 12
elseif iter_no < 2
   iter_no = 2
end
v_ph = v_ph.*1e3; % change phase velocity from km/s in m/s
N_Frames = 50;
f_start = 0;
f_stop = f_stop .* 1e6; % frequency from MHz to Hz
f = linspace(f_start, f_stop, N_Frames); % frequency
alpha_db = alpha_db./1e3; % attenuetion from db/km to db/m
alpha_np = log(10)/20*alpha_db;
r1 = 0; % reflection factor of the source
r3 = 0; % reflection factor of the sink
r4 = -1/3; % reflection factor at junction
gamma = alpha_np + 1j*2*pi*f/v_ph; % transmission coefficient with attenuation and phase coefficent
h = zeros(iter_no,N_Frames);
h_0 = (1-r1)/2 * exp(-gamma*l1);
h(1, :) = h_0 .* (2/3*exp(-gamma*l3)); % first iteration
h_sum = h(1,:); % overall transfer function

% k is 1
t = zeros (1,iter_no);
P_t = zeros (1,iter_no);
t(1) = (l1+l3)/v_ph;
t_period = (2*l2)/v_ph;
P_t(1) = abs(h(1,1))^2;
pdp_acc = t(1)*P_t(1);
h_acc = P_t(1);
ds_acc = (t(1)^2)*P_t(1);
% k is from 2 to iteration number
for k = 2 : iter_no
    h(k, :) = h(1, :) .* (2/3*(r2+eps)^(k-1)*(-1/3)^(k-2)*exp((k-1)*(-gamma*2*l2)));
    h_sum = h_sum + h(k,:);
    t(k) = t(1) + (k-1)*t_period;
    P_t(k) = abs(h(k,1))^2;
    pdp_acc = pdp_acc + t(k)*P_t(k);
    h_acc = h_acc + P_t(k);
    ds_acc = ds_acc + (t(k)^2)*P_t(k);
end
t_d = pdp_acc/h_acc;
t_rms = sqrt(ds_acc/h_acc - (t_d^2));
ratio = min((abs(h_sum(1,:))).^2)/max((abs(h_sum(1,:))).^2);
if ratio >= 1/2
	B_coh = inf;
else
   B_coh = 1/(6.5*t_rms)/1e6;
end

clf
cla;

uni_colors = [238/255 28/255 35/255]; % rot
uni_colors = [uni_colors; 237/255 0 140/255]; % pink
uni_colors = [uni_colors; 140/255 198/255 62/255]; % apfelgruen
uni_colors = [uni_colors; 1 221/255 0]; % gelb
uni_colors = [uni_colors; 0 126/255 198/255]; % mittelblau
uni_colors = [uni_colors; 244/255 111/255 33/255]; % orange
uni_colors = [uni_colors; 19/255 31/255 31/255]; % anthrazit
uni_colors = [uni_colors; 0/255 170/255 173/255]; % turquoise
uni_colors = [uni_colors; 128/255 40/255 145/255]; % lila


ax1 = subplot(2,1,1);
set(ax1, 'XLimMode', 'manual');
set(ax1, 'YLimMode', 'manual');
xlabel('t/ns');
title('Power Delay Profile (Up) and |H_{sum}(f)| (Down)'); 
set(gca,'xtick',t(1):t_period:t(k));
set(gca,'xticklabel',t*10^9)
s1 = sprintf ('%d', t_rms*10^9);
s2 = ['\tau_{rms} = ' s1 'ns'];
% linear scalexlabel('t/ns');
if logscale == 0
set(ax1, 'XLim', [0, max(t(:))]);
set(ax1, 'YLim', [0, 1]);
hold on;
rectangle('Position',[t_d-t_rms/2,0,t_rms,1],'FaceColor','m','EdgeColor','m');
text((t_d+t_rms/2)*1.1,0.9,s2,'color','m');
ylabel('P(t)/P_{Total}');
hold all; grid on;
% k = 1
P_1 = stem(t(1), P_t(1)/h_acc);
set(P_1,'color', uni_colors(1,:), 'markersize', 3);
% k is from 2 to iteration number
for k = 2: iter_no
    P_k = stem(t(k), P_t(k)/h_acc, 'color', uni_colors(mod(k,9)+1,:), 'markersize', 3);
    set(P_k,'color', uni_colors(mod(k,9)+1,:), 'markersize', 3);
end
P_td = stem(t_d,1,'color','b','LineWidth',2);
legend(P_td,'t_d'); 

% log scale
else
set(ax1, 'XLim', [0, max(t(:))]);
set(ax1, 'YLim', [min(10*log10(P_t()/h_acc)), abs(min(10*log10(P_t()/h_acc)))/6]);
hold on;
rectangle('Position',[t_d-t_rms/2,min(10*log10(P_t()/h_acc)),t_rms,abs(min(10*log10(P_t()/h_acc)))],'FaceColor','m','EdgeColor','m');
text((t_d+t_rms/2)*1.1,min(10*log10(P_t()/h_acc))*0.5,s2,'color','m');
ylabel('P(t)/P_{Total} in dB');
hold all; grid on;
% k = 1
P_1 = stem(t(1), 10*log10(P_t(1)/h_acc));
set(P_1,'color', uni_colors(1,:), 'markersize', 3);
% k is from 2 to iteration number
for k = 2: iter_no
    P_k = stem(t(k), 10*log10(P_t(k)/h_acc), 'color', uni_colors(mod(k,9)+1,:), 'markersize', 3);
    set(P_k,'color', uni_colors(mod(k,9)+1,:), 'markersize', 3);
end
P_td = stem(t_d,0,'color','b','LineWidth',2);
legend(P_td,'t_d');
end

ax2 = subplot(2,1,2);
set(ax2, 'XLimMode', 'manual');
set(ax2, 'YLimMode', 'manual');
set(ax2, 'XLim', [0, max(f(:)/1e6)]);
set(ax2, 'YLim', [0, max(abs(h_sum(1,:)))*1.2]);
hold on;
xlabel('f/MHz');
ylabel('|H_{sum}(f)|');
s3 = sprintf ('%0.2f', B_coh);
s4 = ['B_{coh} = ' s3 'MHz']; 
if ratio <=1/2
rectangle('Position',[2/t_period/1e6-B_coh/2,0,B_coh,1],'FaceColor','m','EdgeColor','m');
plot([0,f_stop],[sqrt(1/2)*max(abs(h_sum(1,:))),sqrt(1/2)*max(abs(h_sum(1,:)))],'m');
text((2/t_period/1e6+B_coh/2)*1.05,max(abs(h_sum(1,:)))*1.1,s4,'color','m');
else
plot([0,f_stop],[sqrt(1/2)*max(abs(h_sum(1,:))),sqrt(1/2)*max(abs(h_sum(1,:)))],'m');
text(0, max(abs(h_sum(1,:)))*1.1, s4,'color','m');
end
plot(f(:)/1e6, abs(h_sum(1,:)),'b','Linewidth',2);hold on
print(filename,'-dgif',"-S640,480", "-r0");
