% This script solves the steady state response of a lossless
% transmission line driven by a +/- 10V source square wave.
% The source impedance is 30 ohms, the line characteristic impedance 
% is 50 ohms the line length is 5500 m.  The propagation is 3e8 m/s and
% the fundamental frequency of the square wave is 1e5 Hz.  Forty one
% harmonics are carried in the solution which is accomplished by way
% of the phasor method term-by-term.
% transmission Line with +/-10V square wave source
set(0,'DefaultAxesFontSize',12);
set(0,'DefaultTextFontSize',12);
clear all
clc
Zg=30;
Z0=50;
L=5500;
up=3e8;
nharm=41;
Vg=zeros(1,nharm);Gamma=zeros(1,nharm);betaL=zeros(1,nharm);Zin=zeros(1,nharm);
VL=zeros(1,nharm);Vin=zeros(1,nharm);
f0=1e5;
omega0=2*pi*f0;
disp(['    The Characteristic Impedance Z0 =' num2str(Z0) ' Ohms']);
disp(['    The Source Impedance Zg =' num2str(Zg) ' Ohms']);
disp(['    The Line Length L =' num2str(L) ' m']);
disp(['    The Fundamental Frequency f0 = ' num2str(f0) ' Hz']);
disp(['    The Propagation Velocity up = ' num2str(up) ' m/s']);
disp('      ');
ZL=input('    The Load Resistance in ohms ZL = ');
omega=linspace(1,nharm,nharm)*omega0;up=3e8;
betaL=(L/up)*omega;
omega0t=linspace(0,6*pi,301);
for n=1:2:nharm;
    if ((n+1)/2)==fix((n+1)/2);
        Vg(1,n)=-j*40/(n*pi);% changed sign on this term
    end
    Gamma(1,n)=(ZL-Z0)/(ZL+Z0);
    Zin(1,n)=Z0*((ZL+j*Z0*tan(betaL(1,n)))/(Z0+j*ZL*tan(betaL(1,n))));
    Vin(1,n)=Vg(1,n)*(Zin(1,n)/(Zin(1,n)+Zg));
    VL(1,n)=Vin(1,n)*((1+Gamma(1,n))/(exp(j*betaL(1,n))+Gamma(1,n)*exp(-j*betaL(1,n))));
end
sum1=zeros(1,301);sum2=zeros(1,301);
for n=1:2:nharm
    sum1=sum1+abs(Vin(1,n))*cos(n*omega0t+angle(Vin(1,n))*ones(1,301));
    sum2=sum2+abs(VL(1,n))*cos(n*omega0t+angle(VL(1,n))*ones(1,301));
end
vgoft=10*[0 ones(1,49) 0 -ones(1,49) 0 ones(1,49) 0 -ones(1,49) 0 ones(1,49) 0 -ones(1,49) 0];
vinoft=sum1;
vLoft=sum2;
lim1=1.1*max([vLoft vgoft]);

figure(1);clf;
axis([0 6*pi -lim1 lim1]);
hold on
plot(omega0t,vgoft,'r')
hold on
plot(omega0t,vinoft,'b')
hold on
plot(omega0t,vLoft,'k')
hold on
box on
xlabel('Dimensionless Time, \omega_0t')
ylabel('Voltages, v_g(t), v(_i_n(t), v_l(t)')
text(1,-.9*lim1,'Red = v_g(t)     Blue = v_i_n(t)     Black = v_L(t)');
text(2,.9*lim1,'Press Enter to Continue');
pause
% Now go after the animation
Vofz=zeros(nharm,201);
zoverL=linspace(-1,0,201);
for n=1:2:nharm
    Vofz(n,:)=Vin(1,n)*(exp(-j*betaL(1,n)*zoverL)+Gamma(1,n)*exp(j*betaL(1,n)*zoverL));
    Vofz(n,:)=Vofz(n,:)/(exp(j*betaL(1,n))+Gamma(1,n)*exp(-j*betaL(1,n)));
end
vlineoft=zeros(301,201);
for t=1:301
    sum3=zeros(1,201);
    for n=1:2:nharm
        sum3=sum3+abs(Vofz(n,:)).*cos(n*omega0t(1,t)*ones(1,201)+angle(Vofz(n,:)));
    end
    vlineoft(t,:)=sum3;
end

figure(2);clf;
box on
axis([-1.1 .1 -lim1 lim1]);
hold on
plot([-1 -1],[-15 15]);
hold on
plot([0 0],[-15 15]);
hold on
xlabel('Dimensionless Distance along the Line, z/L')
hold on
ylabel('Line Voltage v(z,t), Volts');
Vplot=plot(zoverL, vlineoft(1,:),'EraseMode','xor');
T=text(-.9,.9*lim1,'Press Enter to Animate');
pause
for t=1:301
    set(Vplot,'Ydata',vlineoft(t,:));
    set(T,'String','                       ');
    pause(.04)
end
text(-.9,.9*lim1,'Press Enter to Continue');

pause

figure(3);clf;
[Z,T]=meshgrid(zoverL,omega0t);
mesh(Z,T,vlineoft)
view(35,30)
text(-.8,0,-23,'Location, z/L','Rotation',-12);
text(0,7,-23,'Time, \omega_0t','Rotation',30);
zlabel('Voltage, v(z,t)')