clear
clf

printconstant = 0;
n = 100;

A = zeros(n,n);

R = n/2;
R_earth = n;

a = 1;
b = 1;

% Initialize Earth model with core at A(1,1), or top left corner.

for i = 1:n
    for j = 1:n
        if sqrt(i.^2 + j.^2) < R
            A(i,j) = 1;
            %         elseif round(sqrt(i.^2 + j.^2)) = R;
            %             A(i,j) = 2;
        elseif sqrt(i.^2 + j.^2) > R_earth;
            A(i,j) = 1;
            %         elseif round(sqrt(i.^2 + j.^2)) = R;
            %             A(i,j) = 2;
        end
    end
end

% Initialize wave at A(n,1), or the bottom left corner.
% wave_vec = [i;j];
wave_vec = [+1;-3];


current_x = 1;
current_y = n-1;

dt = 0.2 % time step

for iter= 1:1000
    A(round(current_x), round(current_y)) = 2;
    
    %if mod(iter,10)==0
    if mod(iter,1)==0 % always true (for now)
        imagesc(A');
        set(gca,'YDir','normal');
        %   colormap(1-gray);
        axis equal;
        axis tight;
        pause(0.01)
        printconstant = printconstant + 1;
        saveas(gcf, ['wave_',sprintf('%03d',printconstant),'.png'])
        
        %        figure
        %        plot(zeros(n,n),ones(n,n))
        %        axis([0 n 0 n])
        %        pause(0.)
        %     printconstant = printconstant + 1;
        %     saveas(gcf, ['wave_',sprintf('%03d',printconstant),'.png'])
    end
    
    %    current_y = n+1-j;
    %    current_x = i+n-1;
    %    current_x = current_x + wave_vec(1)*dt;
    %    current_y = current_y + wave_vec(2)*dt;
    
    magnitude = (current_x.^2 + current_y.^2);
%     fprintf('%f\n',sqrt(magnitude)/R_earth);
    
    normal_vec_surface = (norm(wave_vec))*[current_x; current_y]/(magnitude);
    theta = acosd(norm(normal_vec_surface)/norm(wave_vec));
    
    if sqrt(magnitude) > R && sqrt(magnitude) < R_earth
        current_y = current_y + wave_vec(2)*dt;
        current_x = current_x + wave_vec(1)*dt;
        fprintf('x= %f y= %f vx= %f vy= %f\n',current_x,current_y,wave_vec(1),wave_vec(2))
        
    elseif sqrt(magnitude) <= R;
        
        new_wave_vec = 100*[current_x; current_y]/(magnitude);
        wave_vec = new_wave_vec;
        current_y = current_y + wave_vec(2)*dt;
        current_x = current_x + wave_vec(1)*dt;
        fprintf('Reflection A: x= %f y= %f vx= %f vy= %f\n',current_x,current_y,wave_vec(1),wave_vec(2))
        
    elseif sqrt(magnitude) >= R_earth;
        new_wave_vec = -75.*[-current_x; current_y]/(magnitude);
        wave_vec = new_wave_vec;
        current_y = current_y + wave_vec(2)*dt;
        current_x = current_x + wave_vec(1)*dt;
        fprintf('Reflection B: x= %f y= %f vx= %f vy= %f\n',current_x,current_y,wave_vec(1),wave_vec(2))
    end
    
end



% if mod(iter+9,10)==0
%     f = f+1;
%     figure
%     plot(zeros(n,n),ones(n,n))
%     axis([0 n 0 n])
%     pause(0.01)
%     %     printconstant = printconstant + 1;
%     %     saveas(gcf, ['wave_',sprintf('%03d',printconstant),'.png'])
% end
%
% figure
% plot(zeros(n,n),ones(n,n),'b-')
% axis([0 n 0 n])
