function [ output ] = CARPE3( params )
%
% [ output ] = CARPE3( params ) 


% 3D PE set up and function call

lambda0=params.c0/params.f0;
range_update_interval=2*round(.5*params.range_update_in_meters/lambda0);

wid=lambda0 * (600/2^11)*params.ny ;  
   % horizontal width , domain is TWICE this across
 
steplength=(params.steplamb)*lambda0;
numstep=round(params.rng/steplength) 
range=steplength*numstep ;
disp(['x final range (km) :  ' num2str(range/1000)])

freqs=[ params.f0 ];  % can loop over frequancies for broadband, need to 
              %  change the output tabulation and file writing.

for freq=freqs
  tic 
  freq
  
% propagate x direction     
  x=0; 
  [psifinal,x,y,z,M,psi1,Ez,Af] = propWAPErev11(params.ny,params.nz,...
      numstep,freq,params.c0,params.Envfile, wid, steplength , ...
           params.aspect,params.depP,params.zs,params.ys,params.filen, ...
           params.nout, params.rdcw_flag, range_update_interval);
   
  anorm1=max(max(abs(psi1)));
  anormfin=max(max(abs(psifinal)));

  figure(2);clf
  subplot(211)
  
  finaldb=20*log10(abs(squeeze(psifinal(end,:,:))/anorm1));
  finlin=squeeze(abs(psifinal(end,:,:)));
  wat=find( z > -(params.depP+5) & z < .1);
  %pcolor(y(1:2:end),z(wat),finlin(wat,1:2:end));shading flat
  pcolor(y(1:2:end),z(wat),finaldb(wat,1:2:end));shading flat
  caxis(max(finaldb(:))+[-20 0])
  title('|\psi_{final}| (dB)')
  ylabel('depth (m)')
  xlabel('distance (m)')
  %caxis([ max(finaldb(:))-40 max(finaldb(:))+2])
  % axis equal
  colorbar
  subplot(212)
  %figure(2); clf
  pcolor( y(1:2:end),z(wat),angle(squeeze(psifinal(end,wat,1:2:end))) ) ;shading flat
  title('phase angle \psi_{final}')
  ylabel('depth (m)')
  xlabel('distance (m)')
  %axis equal
  colorbar
  pause(4)


  % eval([ 'save ' outdir 'psi' num2str(freq) ' y z x psifinal psi1 freq c0 range steplength -V6' ])
  toc
end  % freq loop

%anorm2=max(max(abs(psiout2)));
%figure(2)
%imagesc(x,z,20*log10(abs(psiout2)/anorm2),[-80 0]);

% save outputfile z x y psi1 psifinal freq dep zs ... 
%    numstep steplength nz range f0 c0 aspect Af Ez
 
 clear finaldb finlin
 %eval(['save ' outfile]);
 
 xg=(1:numstep)*steplength;
 xgg=xg'*ones(1,params.ny);
 ygg=ones(numstep,1)*y;
 rgg=sqrt(xgg.^2 + ygg.^2) ;
 zgg=ones(numstep,1)*(z(1:params.nz/2)-params.zs);
 rzgg=sqrt(xgg(:,1:params.nz/2).^2 + zgg.^2);
 
 if(exist('Ez'))
     figure(3);clf
     refI=max(max(Ez'.*rgg') );
     pcolor(xg(1:3:end),y(1:3:end), ...
         10*log10((Ez(1:3:end,1:3:end)/refI)'.*rgg(1:3:end,1:3:end)'));shading flat
     caxis([ -40 -10])
     colorbar
     xlabel('x (m)')
     ylabel('y (m)')
     title('depth average Ez','interpreter','none')
 end
 
 if(exist('Af') )
     figure(4);clf
     %pcolor(xg',z(1:nz/2),20*log10(abs(Af).*sqrt(rzgg')));shading flat
     pcolor(xg(1:3:end)',z(1:2:params.nz/2),20*log10(abs(Af(1:2:end,1:3:end))));shading flat
     caxis([ -20*log10(range) -30])
     colorbar
     xlabel('x (m)')
     ylabel('z (m)')
     title('|Af| dB at ys, the source location in the transverse dim.',...
         'interpreter','none')
 end
 
 output.psifinal    =psifinal;
 output.xgg     = xgg ;  
 output.ygg     = ygg;
 output.rgg     = rgg; 
 output.zgg     = zgg;
 output.Ez      = Ez;
 output.Af      = Af;
 
end
 

