% signal.m  calculate a source pulse
function [s] = signal(t,dt,n,v1,npts,na,thalf)
% t = length of wavelet in s
% dt = sample interval
% The approximate predominant frequency is 2*t/n
% Be sure to set dt small enough to avoid aliasing; in general
% dt <= 1/(10*dt)
% n = number of extrema of wavelet (if n<0, wavelet mult by -1)
% v1 = velocity of upper layer or 1
% npts is number of points in wavelet
% na = number of zeros before wavelet (usually = 0)
% To make a wavelet more aysmmetric and more minimum phase, use thalf to
% exponentially taper the later part of the wavelet.
% thalf = if thalf > 0, the wavelet is weighted by an exponential such that at 
%         thalf seconds the amplitude is multipied by 0.5
%         thalf is usually chosen to be ~t or t/2
% s = signal
%
if v1<=0; v1=1; end
factor=1;
if n<0; factor=-1; end
n=abs(n);
s=[];
if na>0; s=[s zeros(1,na)]; end
ns=t/dt + 1;
fn=n;
fm=(fn+2)/fn;
d=fn*pi/t;
tt=-dt;
a=0;
for i=1:ns
   tt=tt+dt;
   t1=d*tt;
   t2=fm*t1;
   c=sin(t1)-sin(t2)/fm;
   s(na+i)=c;
   if c>a; a=c; end
end
a=v1/a;
for i=1:ns
   j=na+i;
   s(j)=a*s(j);
end
se=0;
j=na+ns+1;
for i=j:npts
   s(i)=se;
end
if thalf>0
   n2=na+1;
   for i=n2:npts
      time=(i-n2)*dt;
      x=time*0.696/thalf;
      s(i)=s(i)*exp(-x);
   end
end
smax=max(abs(s));
s=factor*s/smax;