function plotLTP(b,a,tmax,fmax)
% plotLT(b,a,tmax,fmax)
% Plot Laplace transforms
% First, plot a pole zero diagram of H(s)=B(s)/A(s).
% Then plot the inverse Laplace transform h(t) from 0 to tmax,
% and the phase and magnitude of H from 0 to wmax.
% Note: any impulsive parts in h(t) will not be plotted.
    figure
	axis equal tight
	
%  pole-zero plot
	subplot(2,2,2)
    pzd(b,a)
	
% do the frequency domain
	f=[0:fmax/10000:fmax];
    w=2*pi*f;
	H=freqs(b,a,w);
	subplot(2,2,3)
    hold on
    axis square
    plot(f,zeros(1,length(f)),'k')
	plot(f,(1/pi)*(angle(H)),'r')
    grid
	title('phase of H(j\omega)')
	xlabel('f in Hertz')
	ylabel('multiples of pi')
    subplot(2,2,1)
    hold on
    axis square
    plot(f,zeros(1,length(f)),'k')
	plot(f,abs(H),'b')
    grid
	title('magnitude of H(j\omega)')
	xlabel('f in Hertz')
	
%  do the time domain
    subplot(2,2,4);
	t=[0:tmax/500:tmax];
	% The MATLAB function 'impulse' has trouble
	% with poles near the jw axis. 
	L=1024;ts=tmax/L;t=ts*[0:L];
	h=impulse(b,a,t);
	mn=min(0,min(h));
    [hm,km]=max(h);
    hm=max(0,hm);
    d=.1*(hm-mn+.01);
    mn=mn-d;mx=hm+2*d;
	
    plot(t,h,'b',[0 tmax],[0,0],'k',[0 0],[mn,mx],'k')
	hold on
	plot(ts*(km+[-10 40]),[hm,hm],'r',ts*[km km],[mn,mx],'k')
	text(ts*(km+10),hm+d,num2str(hm))
	text(ts*(km+3),mn,num2str(ts*km))
	axis([0 tmax mn mx])
	ttl='Time domain signal, h(t)';
	if length(b)==length(a)
		ttl=[ttl,'  (impulse at t=0 is not shown)'];
	end	
	hold off,xlabel('t in seconds'),axis('off'),title(ttl)
    figure(gcf)