function [psiout,x,y,z,M,psi1,Ez,Af]=propWAPErev11(ny,nz,numstep,fr,c0,envfile,wid,steplength,aspect,depP,zs,ys,filen,nout,rdcw_flag,range_update_steps)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Parabolic equation simulation of 3-D (x,y,z) acoustic propagation.
% An initial condition psi(x_0,y,z) is marched forward using the split-step
% method to a final range psi (range, y, z). 
%
% T. Duda (WHOI)  originally vertical prop version, march 2006
%
% rev 9 has proper half step
% rev 10 loops in k domain
% rev 11 deletes obsolete variables, is set up for range-independent medium.
%    The medium can have only (y,z) variation, but this code can be modified to 
%    change the medium as it steps along (Read in from a saved enf file, for
%    example.)
%
%  psiout: output psi, nout total number, evenly spaced.
%  x: final range, propagation axis direction
%  y: transverse grid position
%  z: z grid
%  M: used for movles, not typically used, obsolete.
%  psi1: start file at range x=0
%  Ez: Flat (x,y) depth-averaged intensity, integral dz psi^2
%  Af: Complex acoustic field psi(x,z) at the source y position
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic
psiout=[];
nnout=round(numstep*(1:nout)/nout)
% persistent rhosq screen Y Z atten fr
lambda=c0/fr ;

nss=60; % plot once every nss steps

% -----------set up grid ....................................................
dx=steplength;
dy=wid/(ny-0.5);
dz=dy/aspect;

ntransz=nz; ntransy=ny;
nstep=numstep;
%z     =  dz/2+(0:(nz-1))*dz;
%zt     =  [-z(nz:-1:1) z];
y     =  (-0.5*ny*dy):dy:(0.5*ny*dy)-dy;
z     =  (-0.5*nz*dz):dz:(0.5*nz*dz)-dz;

[Y,Z]=meshgrid(y,z);
zso1=min(find(z > zs));
zso2=max(find(z < -zs));
yso=min(find(y>ys))
rhosq=(Y.^2 +Z.^2);
d1=( (Y-Y(1,yso)).^2 + (Z-Z(zso1,1)).^2);
d2=( (Y-Y(1,yso)).^2 + (Z-Z(zso2,1)).^2); % verified to be symmetrc about zero

Ly=max(y)-min(y);
Lz=max(z)-min(z);

% --- Green wide angle starting field from FOR3d
k0= 2*pi*fr/c0;  % 2 pi   over lambda = ko
psi=0.5 * k0 * (1.4467-.04201*k0*k0*d1 ) .* exp(-(k0*k0*d1)/3.0512 ) ;
psi=psi-( 0.5 * k0 * (1.4467-.04201*k0*k0*d2 ) .* exp(-(k0*k0*d2)/3.0512 ) ) ;
%   d1 and d2 are distances already squared....

clear d1 d2 rhosq
clear psiz psiy ido dum

% ----------------------          Free Propagator set up ------------------
ky=(2*pi/(Ly))*[0:ny/2];
ky(ntransy:-1:ny/2+1)=-ky(2:ny/2+1);
ky=ones(ntransz,1)*ky;

%kz= aspect*ky';       % make symmetric
kz=(2*pi/(Lz))*[0:nz/2];
kz(ntransz:-1:nz/2+1)=-kz(2:nz/2+1);
kz=ones(ntransy,1)*kz;
kz=kz';

% ----- run environment design code -- produces attbo, cin  -- requires Y and Z

dista=0;
eval(envfile);

hh=max(abs(z));
la400=3.75; % wavelength at 400 Hz (testing done at 400 Hz, domain scales with lambda)

% ------------------        Attenuation function (Sponge layers all around)
alpha=lambda/(50*lambda/la400) ; beta=4*4/dx;  % alpha 4/h  for 4 kHz    % was lamb/2 for 10k hz
attb= exp( -beta*dx*exp( -((z(nz:-1:1)-hh)/(alpha*hh)).^2));
attenz=attb.*fliplr(attb);
attenz=attenz'*ones(1,ny);

hh=max(abs(y));

alpha=lambda/(50*lambda/la400) ; beta=4*4/dx;  % alpha 4/h  for 4 kHz    % was lamb/2 for 10k hz
att= exp( -beta*dx*exp( -((y(ny:-1:1)-hh)/(alpha*hh)).^2));
atteny=att.*fliplr(att);
atteny=atteny'*ones(1,nz);
atteny=atteny';
atten=attenz.*atteny;
clear atteny attenz
atten=atten.*attbo;

% ---------------------         Phase screen
U=(cin-c0)./cin;         
if max(max(isinf(U)))
    max(max(isinf(U)))
    disp('infs ... won''t run')
    return
end
screenm = exp(-i*k0*dx*U);
 % clear U cin

%     standard PE     fr0=exp(-i*0.5*dx*(kz.^2+ky.^2)/k0); 
%WAPE Thomson Chapman 1983 JASA  propagator, 2 half-steps per "range step"
kappasq = (kz.^2+ky.^2);
fr0=exp(-i*0.5*dx*kappasq./ (  (k0^2 -kappasq).^(0.5) + k0 ) ); 
disp('____________ WAPE rev 10 _________________')
clear kappasq

%dista=0;  % moved to before eval envfile
psi1=psi;
% [j1,j2]=size(screen);
% midl=find( abs(yin)< 1e-8 );
% nrng=(midl-(ny/2) ):(midl+ (ny/2)-1);
Ez=ones(nstep,ny);
nzhalf=nz/2;
Af=ones(nzhalf,nstep);

figure(1);clf

tempout=psi;

chat_interval=2*round(.5*nstep/50); % list progress at 2% intervals

if(rdcw_flag)
disp([ 'Environment update interval '   num2str(range_update_steps) ' steps'])
disp([ 'Environment update interval '   num2str(dx*range_update_steps) ' meters'])
else
    disp('* Range (x) independent environment *')
end

for jj=1:nstep   % main prop loop -----------------------------------------------
    
   % screenm=screen(:,nrng);    
    
    %   medium fixed in march-direction x in this code, need something in
    %   this loop to incorporate x-dependent enviromnent. 
    % Can call the  envfile function as one option , i.e.
    %        eval(envfile) 
    % where different things happen for differing jj
    
    % ---------------------         Phase screen renewal.. new every 20
    % steps
    if(mod(jj-1,range_update_steps)==0 & jj>1 )  % first step calculates at zero range.
        eval(envfile)
        U=(cin-c0)./cin;         
        if max(max(isinf(U)))
            max(max(isinf(U)))
            disp('infs ... won''t run')
            return
        end
        screenm = exp(-i*k0*dx*U);
       % clear U cin
    end
    
    
%     if(mod(jj,25) ==0)
%         disp('randomize seafloor')
%         eval(envfile);  %randomize bottom 
%     end  
    
    if(jj==1 | mod(jj,nss)==0)  % plot every nss steps
       wat= find(  z > -(max(dep(:))+15) & z < .5);
       subplot(211)
       pcolor(y(1:2:end),z(wat), abs(tempout(wat,1:2:end)));shading flat;
     %caxis([ .4 1.2]);
       caxis([ max(max(abs(tempout)))*.1 max(max(abs(tempout))) ])
       colorbar('southoutside')
       title('|field(y,z)|')
       subplot(212)
       %pcolor(y(1:2:end),z(wat),angle(screenm(wat,1:2:end)));shading flat
       pcolor(y(1:2:end),z(wat),(U(wat,1:2:end)));shading flat
       xlabel('y (m)')
       ylabel('z (m)')
         % axis equal
        title([ 'Medium exp(-ik0dxU) @ x= ' num2str(dista) ' m'] );
        title([ 'Medium U @ x= ' num2str(dista) ' m'] );
        colorbar('southoutside')
        caxis( [ -.005 0] )
        pause(0.4)
        
    end
   
    % rev 10 prop -------------------------------
    if(jj ==1)
        psi= fft2( psi );
    end
        %psi=ifft2(fr0.*fft2(atten.*screenm.*ifft2(fr0.*fft2(psi))));  
                 % rev 9 yt syntax above, true half=step, without above if 
    psi=fr0.*fft2(atten.*screenm.*ifft2(fr0.*psi));
    tempout=ifft2(psi);
    Ez(jj,:)=sum( abs(tempout).^2 );
    Af(:,jj)=tempout(1:nzhalf,yso);    % field at source y
    if(jj==nstep)
        psi=tempout;
    end
    % ----------------------------------------
    
    dista  = dista + dx;
    seco=toc;
    
    if( mod(jj,chat_interval)==0)
     disp(['to ' num2str(dista) ' m ; percent ' num2str(100*jj/nstep) ';  time: ' num2str(seco) ' sec']) 
    end
        
    ccc=jj-nnout;
    if( min(abs(ccc))==0)
        nuu=find(ccc==0)
        if(nuu==1)
           psiout=shiftdim( psi(1:(1+nz/2),:) ,-1);
        else
           psiout(nuu,:,:)=psi(1:(1+nz/2),:) ;
        end
    end  
end  % main prop loop

x=dista;
% psifinal=psi;
M='off';  % for making movies if you want (and the PC is able) 

%whos
