function P=shot(W,D,T,fs,iifig)
%
% Call P=shot(W,D,T,fs,iifig);
%
% W     = weight in kg ([] for MK82 weight of 1.8 lb);
% D     = depth in m ([] for D=60 ft);
% T     = Max time in s;
% fs    = sample rate in Hz;
%
%
% COPYRIGHT (c) 2008
% 
% Applied Research Laboratories, The University of Texas at Austin
% (ARL:UT)
% All rights reserved.
% 
% This material may be reproduced by or for the U.S. Government
% only. This notice must appear in all copies of this file and its
% derivatives. Further distribution outside the U.S. Government is
% prohibited without the explicit written authorization from the
% Director of the Environmental Sciences Laboratory, ARL:UT.
% 
% THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND,
% EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO THE
% IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
% PURPOSE, OR NONINFRINGEMENT.
% 
% IN NO EVENT WILL ARL:UT BE LIABLE FOR ANY DAMAGES WHATSOEVER ARISING
% OUT OF THE USE OR INABILITY TO USE THIS SOFTWARE, EVEN IF ARL:UT HAS
% BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.
if isempty(W), W='mk82'; end
if isempty(D), D=60*.3042; end
if ischar(W),
    if strcmpi(W,'mk82'),
        W=1.8/2.2;      % 1.8 pounds
    else
        disp('Did not recognize W'); return
    end
end
%
W13=W^(1/3); r0=100; 
PP(1)=3.74e12 * (W13/r0).^1.13;
PP(2)=9.02e11 * W13/r0;
PP(3)=0.22*PP(2);
PP(4)=0.10*PP(2);
PP(5)=0.03*PP(2);
%
t_decay(1)=1.75e-5 * W13 * (W13/r0).^(-0.26);
t_decay(2)=1.48e-4 * W13 / (D + 10.1).^(1/6);
t_decay(3)=1.91*t_decay(2);
t_decay(4)=2.79*t_decay(2);
t_decay(5)=2.79*t_decay(2);

t_bubble(1)=0;
t_bubble(2)=0.21 * W13 / (D + 10.1).^(5/6);
t_bubble(3)=1.71*t_bubble(2);
t_bubble(4)=2.28*t_bubble(2);
t_bubble(5)=2.81*t_bubble(2);
%
dt=1/fs;
t=[0:dt:T]';
P.ndim=1; P.dims=length(t);
P.scales={t};
P.labels={'Time - s'};
P.dlab='\muPa';
P.data=PP(1)*exp(-t/t_decay(1));
P.parms=[W D r0 fs]
for jb=1:4,
    jb1=jb+1;
    P.data = P.data + PP(jb1) * exp(-abs((t - t_bubble(jb1))/t_decay(jb1)));
    iit=find(t>=t_bubble(jb) & t<=t_bubble(jb1));
    PII=-(pi/2)*(PP(jb)*t_decay(jb) + PP(jb1)*t_decay(jb1))/(t_bubble(jb1)-t_bubble(jb));
    P.data(iit)=P.data(iit) + PII*sin(pi*((t(iit)-t_bubble(jb))/(t_bubble(jb1)-t_bubble(jb))));
end
if iifig, view5(P); end
