function [] = impulseOnLineFun(r1,r2,alpha,vph,L,T,rectangular,ena_dyn_sizing,filename)
% Created by Thomas Handte, INÜ, Aug. 2013


%filename = 'test';
%
%rectangular = false; % switch between rectangular and triangular impuse
%ena_dyn_sizing = true; % enable dynamic sizing
%r1 = 0.9 % reflection coefficient at input	
%r2 = -1; % reflection coefficient at output
%alpha = 0.0001; % attenuation of line [Np / m]
%vph = 2E8; % phase velocity in m/s
%L = 1000; % line length in m
%T = 1E-6; % rectangular impulse from -T/2 to T/2

% change phase velocity from km/s in m/s
vph = vph.*1e3;

% change T from µs in s
T = T.*1e-6;

% length of the impulse
L_imp = vph.*T;

% compute the number of partial waves to consider
N = ceil(3+L_imp./(4.*L));

% set simulation time according to parameters
if (r1==0) % consider only a forward and a backward travelling wave
	t = linspace(0,2.*(L+0.25.*L_imp)./vph,30);
elseif (abs(r1)<=(0.25)) % consider two forward and two backward travelling waves
	t = linspace(0,4.*(L+0.125.*L_imp)./vph,30);
else % consider three forward and three backward travelling waves
	t = linspace(0,6.*(L+1./12.*L_imp)./vph,30);
end
	
x = 0:0.1:L;
x = x';

fn_counter = 0;
tic
for t0 = t 
	u_x_t0 = zeros(length(x),1);
	fwd_waves = u_x_t0;
	bwd_waves = u_x_t0;

	for n = 0:N-1
		if rectangular
			fwd_waves = fwd_waves + (r1.*r2).^n.*exp(-	alpha.*(2.*L.*n+x)).*rectangular_impulse(t0-(2.*L.*n+x)./vph,T);
			bwd_waves = bwd_waves + (r1.*r2).^n.*exp(-alpha.*(2.*L.*(n+1)-	x)).*rectangular_impulse(t0-(2.*L.*(n+1)-x)./vph,T);
		else
			fwd_waves = fwd_waves + (r1.*r2).^n.*exp(-	alpha.*(2.*L.*n+x)).*triangular_impulse(t0-(2.*L.*n+x)./vph,T);
			bwd_waves = bwd_waves + (r1.*r2).^n.*exp(-alpha.*(2.*L.*(n+1)-	x)).*triangular_impulse(t0-(2.*L.*(n+1)-x)./vph,T);
		end
	end
	
	u_x_t0 = (1-r1)./2.*fwd_waves + (1-r1)./2.*r2.*bwd_waves;
	
	u_x_t0_fwd = (1-r1)./2.*fwd_waves;
	u_x_t0_bwd = (1-r1)./2.*r2.*bwd_waves;

	f  = figure(1);
	set(f, "visible", "off");
	plot(x,u_x_t0_fwd,'r--','Linewidth',2);
	hold on
	plot(x,u_x_t0_bwd,'k--','Linewidth',2);
	plot(x,u_x_t0,'b-','Linewidth',2);
	hold off
	
	xlabel('x [m]')
	ylabel('Amplitude [norm. to A]')
	title(sprintf('t_0 = %f \\mus',t0./1E-6))
	xlim([0 L]);

	if ena_dyn_sizing
		ylim([-(1-r1) (1-r1)]);
	else
		ylim([-1 1])
	end
	grid on
	
	legend('u_{fwd}(x,t_0)','u_{bwd}(x,t_0)','u(x,t_0)');
	
	fn = strcat(filename,'_',num2str(fn_counter));
	%print(fn,'-dpng',"-S640,480");
	print(fn,'-dgif',"-S640,480", "-r0");
	
	fn_counter = fn_counter + 1;
end
toc
end

function [x] = triangular_impulse(t,T)
% length of impulse, means impulse lasts from -T/2 ... T/2
x = zeros(length(t),1);
x(t<T/2) = 1 + 0.5.*(t(t<T./2)-(T./2))./(T./2);
x(t<-T/2) = 0;
end

function [x] = rectagular_impulse(t,T)
% length of impulse, means impulse lasts from -T/2 ... T/2
x = zeros(length(t),1);
x(t<T/2) = 1;
x(t<-T/2) = 0;
end
