? 
% Import Data Files
loc = importdata('loctim.txt');
bound = importdata('border.xy');
lat = loc(:,1);     %緯度
lon = loc(:,2);     %經(jīng)度
ele = loc(:,3);     %拔高度
tim = loc(:,4);     %地震運(yùn)動(dòng)到達(dá)時(shí)間
% Plot Switzerland
plot(bound(:,1),bound(:,2))
hold on;
% Plot Stations
plot(lon,lat,'*')
time = distance/velocity
t-t0 = sqrt((x-x0)^2 + (y-y0)^2)/v
% 初始猜想
lat0 = 46.9; lon0 = 9;
% 使用紅色圓圈標(biāo)出初始猜想震中
plot(lon0,lat0,'ro')
% 使用綠色圓圈標(biāo)記地震臺(tái)G
plot(lon(6),lat(6),'go')
% 將度數(shù)轉(zhuǎn)換成公里的參數(shù)
latkm = 111.19; lonkm = 75.82; vp = 5.8;
% 初始計(jì)算
t0 = sqrt( ((lon(6)-lon0)*lonkm).^2 + ((lat(6)-lat0)*latkm).^2 )/vp
%初始猜想存儲(chǔ)到數(shù)組m
m(1) = t0; m(2)= lon0; m(3) = lat0;
%循環(huán)遍歷所有地震臺(tái)
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);        % 時(shí)間
        % 存儲(chǔ)結(jié)果到矩陣G
        G(i,3) = -(diffy/D0(i))/vp; G(i,2) = -(diffx/D0(i))/vp; G(i,1) = 1;
end
G*m = d
dm = inv(G'*G)*G'*d'
% 選擇迭代次數(shù)
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
% 最小二乘解
dm = inv(G'*G)*G'*d'
% 將dm變?yōu)槎葦?shù)，因?yàn)樽鴺?biāo)系中用度數(shù)
dm(2) = dm(2)/lonkm;
dm(3) = dm(3)/latkm;
% Update 'm' Vector
m = m+dm'
% 使用藍(lán)色圓圈標(biāo)識(shí)每次迭代結(jié)果
plot(m(2),m(3),'o')
end
%使用黑色菱形框線標(biāo)識(shí)最終結(jié)果
plot(m(2),m(3),'kd')
dm =
1.0e-006 *
-0.0262
0.4140
0.9619
m =
40.2578 7.5580 47.2169