%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% .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: 27.06.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(r2, alpha_db, v_ph, l1, l2, l3, iter_no, filename)

b = 0;
% r2 is from -1 to 1
if r2 > 1
	r2 = 1
elseif r2 < -1
	r2 = -1
end
% alpha_db, v_ph, l1, l2, l3 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
% 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;
N_Rot = 2; % number of rotation of H1(f)
f_start = 0;
t_1 = (l1+l3)/v_ph; 
B_Rot = 1/t_1;
f_stop = B_Rot*N_Rot;
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
for k = 2 : iter_no
    h(k, :) = h(1, :) .* (2/3*r2^(k-1)*(-1/3)^(k-2)*exp((k-1)*(-gamma*2*l2)));
    h_sum = h_sum + h(k,:);
end

clf
ax1 = gca;
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

fn_counter = 0;
for m = 1:N_Frames
    plot(max(abs(h(1,2:end)))*exp(-1j*2*pi*f), 'w'); %fix x scaling
    set(ax1, 'XLimMode', 'manual');
    set(ax1, 'YLimMode', 'manual');
    set(ax1, 'XLim', [-sum(abs(h(:,1))), sum(abs(h(:,1)))]);
    set(ax1, 'YLim', [-sum(abs(h(:,1))), sum(abs(h(:,1)))]);
    hold on; axis equal;
    p = quiver(0,0, real(h(1,m)), imag(h(1,m)), 0);% vector for first iteration
    set(p, 'color', uni_colors(1,:));     
    plot(real(h(1,:)),imag(h(1,:)),'color', uni_colors(1,:)); %circle for first iteration
    z = h (1,m); 
    for k = 2 : iter_no
     	  q = quiver(sum(real(h(1:k-1,m))), sum(imag(h(1:k-1,m))), real(h(k,m)), imag(h(k,m)), 0); 
        set(q, 'color', uni_colors(mod(k,9)+1,:)); 
        plot(z + abs(h(k,:)).* exp(-1j*angle(h(k,:))),'color', uni_colors(mod(k,9)+1,:));
        z = z + h(k,m);
    end
    a = plot(real(h_sum(1,1:m)),imag(h_sum(1,1:m)),'b');
    set(a,'LineWidth',3);
    hold off;
    xlabel('re\{H(f)\}');
    ylabel('im\{H(f)\}');
    title(sprintf('f = %0.1fMHz',f(m)/1e6));
    if b == 1 	 
	 text(-sum(abs(h(:,1))),sum(abs(h(:,1)))*1.05,'Note: Negative input values are inverted for \alpha_{dB}, v_{ph}, L_i', 'color','r');   
	 end
    fn = strcat(filename,'_',num2str(fn_counter));
	 %print(fn,'-dpng',"-S640,480");
	 print(fn,'-dgif',"-S640,480", "-r0");
	 fn_counter = fn_counter + 1;
end;