% seis1d program (L. Braile 6/23/95, modified 1/5/02 and 3/2/03)
%
% seis1d calculates vertical incidence seismograms including
% primary reflections, first order multiples and primaries with
% a multiple path in the first layer only. 
% Reflection coefficients
% are calculated from velocity only.  No transmission losses,
% geometrical spreading or anelasticity are included.
% Multiples can be included or not included.
% Either a ricker or signal wavelet can be used.
%
% input model parameters (nl = number of layers
%	h = layer thickness, v = velocity)
% thicknesses should be in m or km, velocity in m/s or km/s to match.
%
seis1ddata  % input data change for other input or rename input file
% sesi1ddata.m
if nwavelet ~= 1
   nw=nw1;
end
%
ampshift=0.2;	% amplitude separation on plot of adjacent seismograms
%
% set initial parameters
%
n1=nl-1;
rc=zeros(npts,1);
for i=1:n1		
% calculate reflection coefficients and 
% layer travel times
	R(i)=(v(i+1) - v(i))/(v(i+1) + v(i));
	tl(i)=h(i)/v(i);
end
tt=cumsum(tl)*2;	
% calculate two way travel time to top of layers
% 
% plot model
%
d=cumsum(h);	
% depth to top of each layer
d=[0 d];		
for i=1:nl
	z(2*i-1)=d(i);
	z(2*i)=d(i+1);
	vel(2*i-1)=v(i);
	vel(2*i)=v(i);
end
plot(vel,z,'k-')
axis('ij')
axis([vmin vmax 0 maxz])
xlabel('Velocity (km/s)')
ylabel('Depth (km)')
title(mtitle)
pause
figure
%
% calculate reflection coefficient series
%
tt2=tt*2;		
% times for first order multiples
ttf=tt + tt(1);	
% times for primary plus multiple in first layer
for i=1:n1
	np=tt(i)/dt;
	nint=round(np);
	if nint==0
		nint=1;
	end
	if nint<=npts
		rc(nint)=rc(nint) + R(i);	%Coeff for primary reflection
   end
   %
   % Switch to turn multiples on or off
   %
   if nmult == 1
	np=tt2(i)/dt;
	nint=round(np);
	if nint==0
		nint=1;
	end
	if nint<=npts
		rc(nint)=rc(nint) - R(i)^2;	%Coeff for first order multiple
	end
%
	np=ttf(i)/dt;
	nint=round(np);
	if nint==0
		nint=1;
	end
	if nint<=npts
		rc(nint)=rc(nint) - R(i)*R(1);	%Coeff for primary and mult in first layer
   end
   end
end
%Calculate wavelet
w=signal(ts,dt,next,0,nw,0,thalf); % Signal wavelet
if nwavelet == 1
   w=ricker(nw,freq,dt,nsw);  % This is the old Ricker subroutine
end
%
plot(rc)
a=conv(rc,w);
nshift=round(tshift/dt);
%
plot(a(nshift:npts+nw-1))
A(:,1)=rc(1:npts);
for i=1:5
	A(:,i+1)=a(nshift:npts+nshift-1)+ampshift*(i+3);
end
ts=0:dt:(npts-1)*dt;
plot(ts,A,'k-')
axis([0 maxt -2*ampshift 10*ampshift])
xlabel('Two-way Traveltime (s)')
ylabel('Reflection Coefficient Series and Seismograms')
title(ftitle)