%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% .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_im(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 50
if iter_no > 50
	iter_no = 50
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_1 = (l1+l3)/v_ph;
%t_period = (2*l2)/v_ph;
%P_t1 = abs(h(1,1))^2;
%pdp_acc = t_1*P_t1;
%h_acc = P_t1;
%ds_acc = (t_1^2)*P_t1;
% 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_tk = abs(h(k,1))^2;
    %pdp_acc = pdp_acc + t_k*P_tk;
    %h_acc = h_acc + P_tk;
    %ds_acc = ds_acc + (t_k^2)*P_tk;
end
%t_d = pdp_acc/h_acc;
%t_rms = sqrt(ds_acc/h_acc - (t_d^2));
%B_coh = 1/(t_rms+eps)/1e6;

fn_counter = 0;
for m = 1:N_Frames
	 clf
	 cla;
    ax1 = subplot(1,2,1);
    plot(max(abs(h(1,2:end)))*exp(-1j*2*pi*f), 'w'); %for fixing x scaling
    set(ax1, 'XLimMode', 'manual');
    set(ax1, 'YLimMode', 'manual');
    set(ax1, 'XLim', [-max(abs(h_sum(1,:)))*1.3, max(abs(h_sum(1,:)))*1.3]);
    set(ax1, 'YLim', [-max(abs(h_sum(1,:)))*1.3, max(abs(h_sum(1,:)))*1.3]);
    hold on; axis equal;
    set(gca, 'position', [0.1300 0.3100 0.4000 0.4000]);
    xlabel('re\{H_{sum}(f)\}');
    ylabel('left \& right: im\{H_{sum}(f)\}');
    title(sprintf('f = %0.1fMHz',f(m)/1e6)); 
    quiver(0,0, real(h_sum(1,m)), imag(h_sum(1,m)),0);% vector for overall transfer function
    a = plot(real(h_sum(1,1:m)),imag(h_sum(1,1:m)),'b');
    set(a,'LineWidth',3)
    plot([real(h_sum(1,m)) sum(abs(h(:,1)))*1.3], [imag(h_sum(1,m)) imag(h_sum(1,m))], 'g');
    plot(real(h_sum(1,m)), imag(h_sum(1,m)),'go','MarkerEdgeColor','k','MarkerFaceColor',[.49 1 .63],'MarkerSize',3);
    if b == 1 	 
	 text(-sum(abs(h(:,1))),sum(abs(h(:,1)))*1.1,'Note: Negative input values are inverted for \alpha_{dB}, v_{ph}, L_i, f_{stop}','color','r');   
	 end    
    hold off;
    
    ax2 = subplot(1,2,2);
    set(ax2, 'XLimMode', 'manual');
    set(ax2, 'YLimMode', 'manual');
    set(ax2, 'XLim', [0, max(f(:)/1e6)]);
    set(ax2, 'YLim', [-max(abs(h_sum(1,:)))*1.3, max(abs(h_sum(1,:)))*1.3]);
    hold on;
    set(gca, 'position', [0.5703 0.3100 0.4000 0.4000]);
    xlabel('f/MHz');
    %ylabel('im\{H_{sum}(f)\}');
    title(sprintf('Number of Frames = %d/%d', m, N_Frames));
    plot(f(:)/1e6, imag(h_sum(1,:)),'r');
    plot(f(1:m)/1e6, imag(h_sum(1,1:m)),'b','LineWidth',2);
    plot([0 max(f(m)/1e6)], [imag(h_sum(1,m)) imag(h_sum(1,m))], 'g');
    plot(max(f(m)/1e6), imag(h_sum(1,m)),'go','MarkerEdgeColor','k','MarkerFaceColor',[.49 1 .63],'MarkerSize',3);
	 %plot([B_coh,B_coh],[-1,1],'m'); 
    %s1 = sprintf ('%0.2f', B_coh);
    %s2 = ['B_{coh} = ' s1 'MHz']; 
    %text(0, abs(max(h_sum(1,:)))*1.15, s2,'color','m');
    %plot(f(1:m)/1e6, imag(h_sum(1,1:m)),'-bo','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor',[.49 1 .63],'MarkerSize',3);
    % Create arrow across two subplots  
    %x1 = 0.33;
	 %y1 = 0.51+abs(h_sum(1,m))/max(abs(h_sum(1,:)))*0.2;
	 %x2 = 0.5703+(m-1)/(N_Frames-1)*0.4;
	 %y2 = 0.51+abs(h_sum(1,m))/max(abs(h_sum(1,:)))*0.2;
	 %annotation('arrow',[x1 y1],[x2 y2]);
    fn = strcat(filename,'_',num2str(fn_counter));
	 %print(fn,'-dpng',"-S640,480");
	 print(fn,'-dgif',"-S640,480", "-r0");
	 fn_counter = fn_counter + 1;
end;
