% Lossy Transmission Line Animation
% Input data includes the length L in meters, the inductance per unit length Lp, 
% resistance per unit length Rp, capacitance per unit length Cp, conductance
% per unit length Gp, the load impedance ZL, the source impedance Zg, the line
% length L in meters, the phasor source voltage Vg in rectangular form, and the
% frequency of operation f in hertz.

% This m-file was written at the University of Wyoming in the Electrical
% and Computer Engineering Department and is to be distributed without
% cost.
clear all
set(0,'DefaultAxesFontSize',12);
set(0,'DefaultTextFontSize',12);
L=1000;%Line length in meters, input data
Vg=10;%Source voltage in phasor form (rectangular form), input data
f=.3e6;%frequency of driving source in Hz., input data
omega=2*pi*f;
% The following data is for RG58C/U 
% Instrumentation cable.Data taken from the text by Gayle
% Miner, "Lines and Electromagnetic Fields for Engineers,"
% Oxford Press, 1996.
Lp=2.52e-7;% inductance Henrys per meter, input data
Rp=.028;% resistance Ohms per meter, input data
Cp=101e-12;% capacitance Farads per meter, input data
Gp=5.91e-14;%1.8e-10;% conductance Siemens per meter, input data
disp(['   Line Length = ' num2str(L) ' meters'])   
disp(['   Lprime = ' num2str(Lp) ' henrys/m'])
disp(['   Cprime = ' num2str(Cp) ' farads/m'])
disp(['   Rprime = ' num2str(Rp) ' ohms/m'])
disp(['   Gprime = ' num2str(Gp) ' siemens/m'])
disp(['   f = ' num2str(f) ' hertz'])
ZZ=Rp+j*omega*Lp;
YY=Gp+j*omega*Cp;
up=1/sqrt(Lp*Cp);%Lossless propagation velocity, Gp=Rp=0
disp(['   Lossless Propagation Velocity = ' num2str(up) ' m/s'])
Z0=sqrt(ZZ/YY);% Characteristic impedance
disp(['   Z0 = ' num2str(Z0) ' ohms'])
gamma=sqrt(ZZ*YY);%propagation constant
lambda=2*pi/(imag(gamma));% wavelength in meters
disp(['   Wavelength = ' num2str(lambda) ' m'])
disp('    ')
ZL=input('   The Load Impedance is ZL = ');
%ZL=100;%Load impedance (rectangular form)
Zg=50;%generator impedance (rectanggular form)
%Calculate these from input data
zoverL=[-1:.01:0];
Zin=Z0*((ZL+Z0*tanh(gamma*L)))/(Z0+ZL*tanh(gamma*L));% input impedance
Gamma=(ZL-Z0)/(ZL+Z0);%reflection coefficient
Strg1=['   The Reflection Coefficient = ',num2str(Gamma)];
disp(Strg1)
Strg3=['   The Input Impedance to the Line = ', num2str(Zin) ' ohms'];
disp(Strg3)
Fact=(Vg*Zin/(Zg+Zin))/(exp(gamma*L)+Gamma*exp(-gamma*L));
Vofz=Fact*(exp(-gamma*L*zoverL)+Gamma*exp(gamma*L*zoverL));
VSWR=(max(abs(Vofz)))/(((min(abs(Vofz)))));
Strg2=['   The Voltage Standing Wave Ratio = ',num2str(VSWR)];
disp(Strg2)
lim=max(abs(Vofz));

figure(1);clf;%Plot the so-called standing wave pattern
lim=max(abs(Vofz));
axis([-1 0 0 1.2*lim]) 
hold on
plot(zoverL,abs(Vofz));
box on
xlabel('Dimensionless Distance, z/L')
ylabel('Standing Wave Amplitude')
text(-.9,.1*lim,Strg1)
text(-.9,.2*lim,Strg2)
text(-.9,1.1*lim,'Press Enter to Continue')
pause
mag=abs(Vofz);
angl=angle(Vofz) ;
omegat2=linspace(0,4*pi,101);
Vin=Vg*Zin/(Zin+Zg);
VL=Vin*((1+Gamma)/(exp(gamma*L)+Gamma*exp(-gamma*L)));
Vint=abs(Vin)*cos(omegat2+angle(Vin)*ones(1,101));
VLt=abs(VL)*cos(omegat2+angle(VL)*ones(1,101));
Vgt=abs(Vg)*cos(omegat2+angle(Vg)*ones(1,101));

figure(2);clf;%Plot vg(t), vin(t) and vL(t)
lim2=max(abs([Vg VL Vin]));
axis([0 4*pi -1.1*lim2 1.1*lim2])
hold on
plot([0 omegat2(1,101)], [0 0])
hold on
box on
plot(omegat2,Vint,'--',omegat2,VLt,'-.',omegat2,Vgt,'.')
xlabel('Dimensionless Time, \omega t')
ylabel('v_g(t), v_i_n(t), v_L(t), volts')
text(5,-.6*lim2,'v_i_n  - - - -','FontSize',12)
text(5,-.8*lim2,'v_g . . . .','FontSize',12)
text(5,-lim2,'v_L . - . - . -','FontSize',12)
text(.4*pi,.9*lim2,'Press Enter to Continue')
pause

% calculate the animation data
omegat=linspace(0,8*pi,801);
framedata=zeros(801,101);
for i=1:801
    framedata(i,:)=mag.*cos(omegat(1,i)*ones(1,101)+angl);
end


figure(3);clf;%Plot v(z) for various times
axis([-1.2 .2 -1.2*lim 1.2*lim])
hold on
box on
xlabel('Dimensionless Distance on the Line, z/L')
ylabel('Line Voltage, v(z,t)')
thnd1=text(-.9,1.1*lim,'Press Enter to Animate');
text(-1,1.1*lim,'                      ')
set(thnd1,'String','                               ')
for i=1:10:201
    %omegat=2*pi*i/25;
    %Voft=mag.*cos(omegat*ones(1,101)+angl);
    plot(zoverL,framedata(i,:),'k');
    hold on
    %set(M,'Ydata',[Voft(1,1),Voft(1,101)])
    pause(.4)
end
text(-.9,1.1*lim,'Press Enter to Animate')
pause



figure(4);clf;% Animate v(z,t)
L1=plot(zoverL,framedata(1,:),'k','EraseMode','xor');
hold on
%hold on
plot([-1 -1],[-1.2*lim 1.2*lim])
hold on
plot([0 0],[-1.2*lim 1.2*lim])
hold on
M=plot([-1 0],[framedata(1,1) framedata(1,101)],'o','EraseMode','xor');
axis([-1.2 .2 -1.2*lim 1.2*lim])
xlabel('Dimensionless Distance on the Line, z/L')
ylabel('Line Voltage, v(z,t)')
thnd1=text(-.9,1.1*lim,'Press Enter to Animate');
pause
text(-1,1.1*lim,'                      ')
set(thnd1,'String','                               ')

for i=1:801
    %omegat=8*pi*i/800;
    %Voft=mag.*cos(omegat*ones(1,101)+angl);
    set(L1,'Ydata',framedata(i,:))
    %surfdata=[surfdata;Voft];
    set(M,'Ydata',[framedata(i,1),framedata(i,101)])
    pause(.01)
end
text(-.9,1.1*lim,'Press Enter to Continue')
pause

figure(5);clf;
tdata=linspace(0,8*pi,101);
threeddata=[];
for i=1:8:801
    threeddata=[threeddata;framedata(i,:)];
end
[X,Y]=meshgrid(tdata,zoverL);
%mesh(X,Y,threeddata','Edgecolor',[0 0 0]);
%newmapdata=colormap;
%newmapdata=ones(size(newmapdata));
%colormap=newmapdata;
mesh(X,Y,threeddata')
%colormap([0 0 0])
axis([0 8*pi -1 0  -1.2*lim 1.2*lim])
view(45,30)
text(8*pi,-.8, -1.6*lim,'Dimensionless Distance, z/L','Rotation', 20)%x-axis label
text(3,-1, -1.6*lim,'Dimensionless Time,  \omegat','Rotation',-18)
zlabel('Line Voltage, v(z,t)')
text(1,0,1.1*lim,'Press Enter to Continue', 'Rotation',-18)
pause