% lossytransmline.m Last modified 7/18/2008 % 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') hold off 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) hold off 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') hold off 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'); hold off pause text(-1,1.1*lim,' ') set(thnd1,'String',' ') for i=1:801 set(L1,'Ydata',framedata(i,:)) set(M,'Ydata',[framedata(i,1),framedata(i,101)]) pause(.01) end text(-.9,1.1*lim,'Press Enter to Continue') hold off 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) hold off pause