function M=xlineani(R_source,R_load,vin,vlineinit,vfinit,steps)
%Time-domain transmission line animation.   xlineani(R_source,R_load,vin);
%Source and load resistances are normalized to line characteristic impedance.
%Vin is applied at the left through R_source; on the right it is terminated by R_load.
%Vin is a vector of time points up to 200 points long. 
%(If Vin is shorter than 200 points, it will be padded to zero).
%
%OPTIONAL INITIAL CONDITION PARAMETERS
%xlineani(R_source,R_load,vin,VLINEINIT,VFINIT)
%vlineinit is an optional argument for initial voltage distribution on the line, and 
%vfinit is an optional argument for initial forward travelling wave voltage.  
%Making vlineinit=vfinit will produce initial conditions with no reverse travelling wave; 
%specifying only vlineinit will produce equal waves in both directions.
%Both are up to 40 points, and will be padded with zeros to fill as needed.
%
%OPTIONAL NUMBER OF STEPS TO RUN
%xlineani(R_source,R_load,vin,0,0,Steps)
%Steps is the number of animation steps, which defaults to 200.
%xlineani(R_source,R_load,vin,VLINEINIT,VFINIT,Steps)
%allows intial conditions on the line too, as discussed above.
%
%Note that you can:
%    -drag the legend out of the way.
%    -hit cntl-C from the matlab command window to stop the animation.



%To do: optional speed, optional lengths, legends on plot
%Charles R. Sullivan, 2/18/00
%charles.r.sullivan@dartmouth.edu
%% must fix :  doesn't work with both plots
%OPTIONAL OUTPUT M=xlineani(....)
%The animation can then be run, more smoothly, with 
%>>movie(M)

if nargin >= 6
   anilen=steps;
else
   anilen=200;
end
llen = 40;
x=1:llen;


%Fix length of vin to match anilen
vin = trimpad(vin,anilen);

% Set up vf and vr according to which inputs specified.
if nargin < 4  % zero line init condit
   vf = zeros(1,llen);
   vr = zeros(1,llen);
else
   if nargin >= 5 % have vf and vlineinit
      vf=trimpad(vfinit,llen);
      vlineinit=trimpad(vlineinit,llen);
      vr = vlineinit-vf;
   else   % have only vlineinit
      vf = trimpad(vlineinit,llen)/2;
      vr = vf;
   end
end

   
%Calculate reflection parameters
rhol = (R_load-1)/(R_load+1);
rhos = (R_source-1)/(R_source+1);
div = 1/(1+R_source); %voltage divider coef at input.

%set up initial plots
subplot(2,1,1)
h=plot(x,vf,'>:',x,vr,'<:',x,vf+vr,'+-');
axis([1 llen -2 2])
ylabel('V')
subplot(2,1,2)
h2=plot(x,vf,'>:',x,-vr,'<:',x,vf-vr,'+-');
ylabel('I')
legend('forward','backward','total')
axis([1 llen -2 2])
   
   
%Loop to generate animation (this is it!)
for ti=1:anilen
   set(h(1),'Ydata',vf);
   set(h(2),'Ydata',vr);
   set(h(3),'Ydata',vf+vr);
      set(h2(1),'Ydata',vf);
   set(h2(2),'Ydata',-vr);
   set(h2(3),'Ydata',vf-vr);   
   %Show and/or collect plots
   if(nargout > 0)
      M(ti)=getframe;
    else
       pause(.001)
    end
    
   %Advance a time step
   vfnew = [rhos*vr(2) vf(1:(llen-1))];
   vr = [vr(2:llen) rhol*vf(llen-1)];
   vf = vfnew;
   vf(1) = vfnew(1) + div*vin(ti);
end

%%%%%%%%%%%%%%%%%%%END OF MAIN FCN.

function trimmed = trimpad(v,len)
%trim to vector v length len, or pad with zeros

if length(v)>=len
             trimmed=v(1:len);
else
             trimmed=[v zeros(1,len-length(v))];
end