% plotting voltage standing waves on a transmission line
% terminated with a complex load ZL
clear,clf
L=2.5;   % line length
dz=0.02;
freq=.1e9;  % frequency
wavel=3e8/freq;
alpha=0;    % attenuation constant
gam=alpha+j*2*pi/wavel;  % prop. constant
omega=2*pi*freq;
Zo=50;     % characteristic impedance of line
ZL=10;      % load impedance
Refl=(ZL-Zo)/(ZL+Zo);    % reflection coefficient
t1=0;
tmax=3*1e-9;  % time window 
t2=tmax;
% time step should not be > 0.1 wavelength
dt=1/freq*0.05;
t=[0:dt:tmax];
Nt=length(t);
disp(['number of time steps in window= ',num2str(Nt)])
for it=1:Nt
   ejt=exp(j*omega*t(it));
   ct=0;
   for z=0:dz:L
	ct=ct+1;
	Z(ct)=-z;
	vinc=exp(gam*(z));   % incident wave phasor
   vref=Refl*exp(-gam*(z));   % reflected wave phasor
	V(ct)=real((vinc+vref)*ejt)*exp(-alpha*L);
   Vi(ct)=real(vinc*ejt)*exp(-alpha*L);
   Vr(ct)=real(vref*ejt)*exp(-alpha*L);
 end
 plot(Z,V,Z,Vi,'r-',Z,Vr,'g-')
 legend('Total','Incident','Reflected')
 axis([min(Z),max(Z),-2,2])
 xlabel('z'),ylabel('V(t)')
 title(['ZL=',num2str(ZL),',  REFL COEFF=',num2str(Refl)])
 text(min(Z),-1.8,'SOURCE')
 text(max(Z)-.08*L,-1.8,'LOAD')
 pause(0.01)
 M(:,it)=getframe;
 end

% play the movie
 movie(M,1)