%dipping reflector

function [reflections] = diprefl(S,receivers,time,layer,velocity1)

segments = length(layer)-1;

for i = 1:segments
    subplot(2,1,1), plot (S(1),S(2) ,'rx');
    hold on;
    Lt = layer(i,:);
    Rt = layer(i+1,:);
    d = Rt - Lt; 
    a = Lt - S;
    gamma = -((dot(a,d))/(dot(d,d))); 
    SP = (a + d*gamma);
    
    intersect = SP + S;
    Sprime = 2*SP + S;
    
    %A1, B1, C1 represent the reflective layer
    A1 = Rt(2) - Lt(2);
    B1 = Lt(1) - Rt(1);
    C1 = A1 * Lt(1) + B1 * Lt(2);
    
  for j = 1:length(receivers)
    R = [receivers(j) 0];
    subplot(2,1,1), plot (R(1),R(2) ,'go');
    hold on;
    
    %parameterize reflection as Ax +By = C
    A2 = R(2) - Sprime(2);
    B2 = Sprime(1) - R(1);
    C2 = A2 * Sprime(1) + B2 * Sprime(2);
    
    det = A1*B2 - A2*B1;
    % This is broken down to handle the lines are parallel case which I
    % think could still occur when the source was on a line through the
    % reflector
    if(det == 0)
        %Lines are parallel so set intercept2 to Receiver - wavelet will
        %not plot a zero time reflection
        intercept2(1) = R(1);
        intercept2(2) = R(2);
    else
        intercept2(1) = (B2*C1 - B1*C2)/det;
        intercept2(2) = (A1*C2 - A2*C1)/det;
    end;    
   
    %if reflection actually occurs 
    if (intercept2(1) > Lt(1) && intercept2(1) < Rt(1) )
        subplot(2,1,1), plot([S(1) intercept2(1)  ],[S(2) intercept2(2)  ] , 'r-');
        subplot(2,1,1), plot([ intercept2(1) R(1) ],[ intercept2(2) R(2) ] , 'g-');
        hold on;
        lengthSprime = ((intercept2(1) - R(1))^2 + (intercept2(2) - R(2))^2)^0.5;
        ref_time(j) = lengthSprime / velocity1;
        
    else
        ref_time(j) = 0;
    end;
end;
reflections(i,:) = ref_time;
end;