load wavelet.mat
%thickness (m)
z=[100 200 200];
%velocity of each layer (m/ms)
v=[0.45 2 3 ];
v4=6;
%critical ray parameter
p_crit_1=1/v(2);
p_crit_2=1/v(3);
p_crit_3=1/v4;

[t_int(1) xcrit(1) t0(1)]=one_layer(p_crit_1,z(1),v(1));

[t_int(2) xcrit(2)  t0(2)]=one_layer(p_crit_2,z(1:2),v(1:2));
[t_int(3) xcrit(3) t0(3)]=one_layer(p_crit_3,z,v);

%root mean square
v_rms_1= sqrt( sum( v(1).^2 .* t0(1)) ./ sum(t0(1)) )
v_rms_2= sqrt( sum( v(1:2).^2 .* t0(1:2)) ./ sum(t0(1:2)) )
v_rms_3= sqrt( sum( v(1:3).^2 .* t0(1:3)) ./ sum(t0(1:3)) ) 

%start task 2

x=-1000:10:1000;

t_direct=abs(x)/v(1);  % direct wave travel time
t_reflect_1 = trefl(t0(1),v_rms_1,x);
t_refract_1=trefr(t_int(1),v(2),x);

t_reflect_2 = trefl(t0(2),v_rms_2,x);
t_refract_2=trefr(t_int(2),v(3),x);

t_reflect_3 = trefl(t0(3),v_rms_3,x);
t_refract_3=trefr(t_int(3),v4,x);

figure(1)
hold on
plot(x,t_direct,'b-')
plot(x,t_reflect_1,'r-')
plot(x,t_refract_1,'g-')
plot(x,t_reflect_2,'r-')
plot(x,t_refract_2,'g-')
plot(x,t_reflect_3,'r-')
plot(x,t_refract_3,'g-')

 amp=18;
 amp_refr=10; % amplitude to plot refraction
  t=0:2:1200;
  
  figure(2)
  hold on
  
for i=1:length(x)       % number of geophone
  
  %put direct wave into new trace
  trace1 = wavelet(amp,t_direct(i),t);
  
  %add reflections from each interface
  trace1=trace1+wavelet(amp,t_reflect_1(i),t);
  trace1=trace1+wavelet(amp,t_reflect_2(i),t);
  trace1=trace1+wavelet(amp,t_reflect_3(i),t);
  
  %add refractions from each interface
  if abs(x(i))>xcrit(1)
       trace1=trace1+wavelet(amp_refr,t_refract_1(i),t);
  end
   if abs(x(i))>xcrit(2)
       trace1=trace1+wavelet(amp_refr,t_refract_2(i),t);
   end
   if abs(x(i))>xcrit(3)
       trace1=trace1+wavelet(amp_refr,t_refract_3(i),t);
   end
  
    plot(trace1+x(i),t)
    hold on
    
  section(:,i) = trace1;
  
end   
    title('Shot Gather')
    xlabel('Geophone Spread (m)')
    ylabel('Time(ms)')
    
  
save lab4.mat section