function [dxmax, nxmin, thetamax, Dillum, memin, Dsp, Drp] = wtpropgeom(lambda, z, Ds, Dr, r0sphs, r0sphr)
% Usage: [dxmax, nxmin, thetamax, Dillum, memin, Dsp, Drp] = wtpropgeom(lambda, z, Ds, Dr, r0sphs, r0sphr)
%
% Description:
%   This function calculates limits on propgation sampling for a given
%   propagation scenario. In particular, given the propagation geometry and
%   information about the intervening propgation medium, this routine
%   returns the maximum allowable grid spacing and the number of grid
%   points at that grid spacing necessary to cover the regions of interest.
%   Comments at the end of this m-file comprise the informal write-up which
%   motivates the equations used.
%
% Inputs:
%   lambda: Wavelength of light (m)
%   z: Propagation distance (m)
%   Ds: Diameter of region  of interest (usually the aperture diameter) in the source plane (m)
%   Dr: Diameter of region of interest in the receiver plane (m)
%   r0sphs: Spherical r0 at the source (m)
%   r0sphr: Spherical r0 at the receiver (m)
%
% Outputs:
%   dxmax: Maximum grid spacing (m)
%   nxmin: Minimum number of grid points assuming dx spacing above
%   thetamax: Maximum half-angle at which light is scattered in the path from source to receiver (rad.)
%   Dillum: Diameter of the illuminated region at the receiver (m)
%   memin: Minimum allowable extent of the mesh (m)
%   Dsp: Diameter of the effective region of interest in the source plane (m)
%   Drp: Diameter of the effective region of interest in the receiver plane (m)
%
if (~exist('lambda','var')) lambda = []; end;
if (isempty(lambda))
   lambda = 1.0e-06;
end;
if (~exist('z','var')) z = []; end;
if (isempty(z))
   z = 1.0e+06;
end;
if (~exist('Ds','var')) Ds = []; end;
if (isempty(Ds))
   Ds = 1.0;
end;
if (~exist('Dr','var')) Dr = []; end;
if (isempty(Dr))
   Dr = 1.0;
end;
if (~exist('r0sphs','var')) r0sphs = []; end;
if (isempty(r0sphs))
   r0sphs = Inf;
end;
if (~exist('r0sphr','var')) r0sphr = []; end;
if (isempty(r0sphr))
   r0sphr = Inf;
end;
%
if (~isinf(r0sphr))
   Dsp = Ds + (lambda*z/r0sphr);
else
   Dsp = Ds;
end;
if (~isinf(r0sphs))
   Drp = Dr + (lambda*z/r0sphs);
else
   Drp = Dr;
end;
%
dxmax = (lambda*z)/(Dsp+Drp);
thetamax = lambda/(2*dxmax);
Dillum = Ds+(2*thetamax*z);
memin = (Dr+Dillum)/2;
nxmin = ceil(memin/dxmax);
%
% In the case of plane wave propagation, the mesh spacing requirement is as
% follows:
%  
% dx <= lambda*z/(Dsrc' + Drcvr')
%  
% where the primes refer to the "apparent" diameters of the source and
% receiver as seen from one another, taking into account blurring.  The
% amount of blurring in each case is inversely proportional to the
% spherical wave coherence diameter, corresponding to the given viewing
% direction:
%  
% Dsrc' ~= Dsrc + lambda*z/r0_sph_rcvr Drcvr' ~= Drcvr +
% lambda*z/r0_sph_src
%  
% (although of course in reality the blur does not have a hard edge)
%  
% Note that when blurring is strong and the source and receiver diameters
% are small (by comparison to the blurring), the constraint on dx
% simplifies as follows:
%  
% dx <= lambda*z/(Dsrc' + Drcvr') ~= r0_sph_rcvr + r0_sph_src
%  
% Now with respect to the size of the mesh, it just has to be big enough to
% ensure that we won't see any aliasing effects inside the receiver
% aperture (it's OK if there are some outside).  In general, unless we've
% pre-filtered the field, energy will spread outward from the edges of the
% source aperture at the maximum angle representable on the mesh:
%  
% theta_max = +/- lambda/(2*dx)
%  
% therefore at the rcvr the illuminated region will be of size:
%  
% D_illum = Dsrc + 2*theta_max*z = Dsrc + lambda*z/dx
%  
% If the mesh size is at least that big, there should be no wrap around.
% You can go smaller, in which case you'll get wraparound , but that's OK,
% so long as the the following holds:
%  
% mesh_extent = nx*dx >= (D_illum + Drcvr)/2
%  
% Alternately, you can apply a spatial filter and thus reduce theta_max,
% and D_illum, and thus the required number of mesh points for a given dx
% BUT THIS IS ONLY ALLOWED IF DX IS SMALLER THAN NECESSARY, SO THAT THE
% MAXIMUM ANGLE, LAMBDA/DX, IS BIGGER THAN NECESSARY, SO WE CAN AFFORD TO
% APPLY A FILTER.
%  
% If you break a propagation into two or more pieces, all of the above
% applies to each piece individually, where you would calculate the
% spherical wave coherence diameters in each direction just for that piece.
% It's a little tricky to define the "sizes" of the intermediate
% pseudosources and pseudoreceivers, but it can be done, bearing in mind
% what we know about the sizes of the actual sources and receivers, and the
% angular spectrum of the light, which of course broadens each time we hit
% a phase screen.  (note however that it doesn't change during a vacuum
% propagation)
%  
% This fact - that you can break up the path and apply the same analyses -
% was something I was aware of, but hadn't thought about much, until a
% recent problem was posed.  In this particular case, the highly nonuniform
% distribution of turbulence along the path means that breaking up the
% problem can make a very big difference, as Bill Bradford has pointed out.
%
% Steve Coy & Bob Praus
% MZA Associates Corporation
% voice: (505)245-9970, ext. 11
% fax: (505)245-9971



