function [pd]=plot4(NX,NZ,NT,dx,dt,c,RES,BVEL,res,REALDATA,SYNDATA,nit,pd,od)
  subplot(311);XL='Offset (ft)';ZL='Depth (ft)';
plot([2:nit],res(2:nit),'*',[2:nit],res(2:nit));title('L2 Seismogram Residual vs Iteration Number');
 xlabel('Iteration Number');ylabel('Residual'); subplot(312);
 plot1(NX,NZ,dx,dx,BVEL,['Reconstructed Velocity Model at Iteration = ',num2str(nit)],XL,ZL);subplot(313);
ts=1;te=NT;time=sqrt([ts:te]*dt).^(5.5);
pl1=min([REALDATA(NX/2,ts:te)/max(REALDATA(NX/2,ts:te)).*time SYNDATA(NX/2,ts:te)/max(SYNDATA(NX/2,ts:te)).*time]);
pl2=max([REALDATA(NX/2,ts:te)/max(REALDATA(NX/2,ts:te)).*time SYNDATA(NX/2,ts:te)/max(SYNDATA(NX/2,ts:te)).*time]);
plot([ts:te]*dt,REALDATA(NX/2,ts:te)/max(REALDATA(NX/2,ts:te)).*time ,'--');hold on;ylabel('Amplitude');
plot([ts:te]*dt,SYNDATA(NX/2,ts:te)/max(SYNDATA(NX/2,ts:te)).*time,'-');ylabel('Amplitude'); xlabel('Time (s)')
hold off
pd(:,nit)=SYNDATA(NX/2,:)'.*time';
od(:,nit)=REALDATA(NX/2,:)'.*time';
title('Observed (--) vs Predicted (-) Gained Seismograms');  
axis([ts*dt te*dt pl1 pl2]);
hold off;
