% Finding Earthquake Location from Recorded Data

% Import Data Files
loc = importdata('loctim.txt');
bound = importdata('border.xy');
lat = loc(:,1);     % latitude
lon = loc(:,2);     % longitude
ele = loc(:,3);     % elevation
tim = loc(:,4);     % time of arrival 

% Plot Switzerland
plot(bound(:,1),bound(:,2))
hold on;
% Plot Stations
plot(lon,lat,'*')
% Initial Guess
lat0 = 46.9; lon0 = 9;
plot(lon0,lat0,'ro')
plot(lon(6),lat(6),'go')

% Parameters to change degrees to kilometers
latkm = 111.19; lonkm = 75.82; vp = 5.8;


% Initial Calculation
t0 = sqrt( ((lon(6)-lon0)*lonkm).^2 + ((lat(6)-lat0)*latkm).^2 )/vp


m(1) = t0; m(2)= lon0; m(3) = lat0;

for i = 1:6
    for i=1:length(lat)
        diffx = (lon(i)-m(2))*lonkm;
        diffy = (lat(i)-m(3))*latkm;
        D0(i) = sqrt(diffx^2+diffy^2);
        d(i) = tim(i) - D0(i)/vp - m(1);
        G(i,3) = -(diffy/D0(i))/vp; G(i,2) = -(diffx/D0(i))/vp; G(i,1) = 1;
    end
% Least Squares Solution
dm = inv(G'*G)*G'*d'
dm(2) = dm(2)/lonkm;
dm(3) = dm(3)/latkm;
% Update 'm' Vector
m = m+dm'
% Iterative Solution
plot(m(2),m(3),'o')
end

% Final Solution
plot(m(2),m(3),'kd')