%%
%Author Edwin CASTRILLO
%assigning values
A = importdata('arrival_times.txt', '\t', 1);
load val.txt
tp=A.data(:,4);
ts=A.data(:,5);
v2=[A.data(:,7);A.data(:,8)];
name = A.textdata(:,1);
xco = [A.data(:,1);A.data(:,1)];
yco = [A.data(:,2);A.data(:,2)];
zco = [A.data(:,3)/1000;A.data(:,3)/1000];
%Definition of data vector matrix, Kernel matrix, and others
G=zeros(length(xco),4);
dobs=zeros(length(xco),1);
dtdx=zeros(length(xco),1);
dtdy=zeros(length(xco),1);
dtdz=zeros(length(xco),1);
dtdt=zeros(length(xco),1);
arrvt=zeros(2*length(xco),1);
ttt=zeros(length(xco),1);
dpred=zeros(length(xco),1);
d=zeros(length(xco),1);
res=zeros(length(xco),1);
mic=zeros(100000,4);
dmn=zeros(4,1);
mestn=zeros(4,1);
OT=187.2;
%xg=150;
%yg=80;
xg=val(1,1);
yg=val(1,2);
zg=10;
tg=OT;
mi=[xg;yg;zg;tg];
cont=0;
arrvt=[tp;ts];
% Plotting stations locations
figure(1);
plot(xco(:,1),yco(:,1),'b*')
% Plotting the initial guess as a green circle
figure(2);
plot(xco(:,1),yco(:,1),'b*')
hold on
plot(xg,yg,'go')
% definition of the intial minimum error
E=10;
%begin the iterations
while min(abs(E))>(1*10^0)
cont=cont+1;
if cont<=1
dobs(:)=[tp(:);ts(:)];
for i=1:18
alpha=v2(i);
% calculating the partial derivatives
dtdx(i)=(-(xco(i)-mi(1)));
dtdx(i)=dtdx(i)/(alpha^2*(arrvt(i)-mi(4)));
dtdy(i)=(-(yco(i)-mi(2)));
dtdy(i)=dtdy(i)/(alpha^2*(arrvt(i)-mi(4)));
dtdz(i)=(-(zco(i)-mi(3)));
dtdz(i)=dtdz(i)/(alpha^2*(arrvt(i)-mi(4)));
dtdt(i)=1;
% travel time calculation
ttt(i)=sqrt(((xco(i)-xg)*(xco(i)-xg))+((yco(i)-yg)*(yco(i)-yg))+((zco(i)-zg)*(zco(i)-zg)))/alpha;
dpred(i)=OT+ttt(i);
d(i)=dobs(i)-dpred(i);
end
G=[dtdx dtdy dtdz dtdt];
% Getting Eigenvalues and Eigenvectors of G
[U,S,V] = svd(G);
% Least square solution
dm=inv(G'*G)*G'*d;
% updating the 'mi' vector
mi=mi+dm;
mic(cont,:)=(mi(:))';
dc=G*mi;
E=(sum(d(:)).^2);
% Plot first solution as a cyan circle
figure(3);
plot(xco(:,1),yco(:,1),'b*')
hold on
plot(xg,yg,'go')
hold on
plot(mi(1),mi(2),'co')
% Creating a new figure for the next stage and plotting
% and Plot first solution as a cyan circle
figure(4);
plot3(mi(1),mi(2),-mi(3),'co')
hold on
% evalauting condition before starting next iterations
elseif cont>1
dobs(:)=[tp(:);ts(:)];
for i=1:18
alpha=v2(i);
% calculating the partial derivatives
dtdx(i)=(-(xco(i)-mi(1,1)));
dtdx(i)=dtdx(i)/(alpha^2*(arrvt(i)-mi(4)));
dtdy(i)=(-(yco(i)-mi(2,1)));
dtdy(i)=dtdy(i)/(alpha^2*(arrvt(i)-mi(4)));
dtdz(i)=(-(zco(i)-mi(3,1)));
dtdz(i)=dtdz(i)/(alpha^2*(arrvt(i)-mi(4)));
dtdt(i)=1;
% travel time calculation
ttt(i)=sqrt(((xco(i)-mi(1,1))*(xco(i)-mi(1,1)))+((yco(i)-mi(2,1))*(yco(i)-mi(2,1)))+((zco(i)-mi(3,1))*(zco(i)-mi(3,1))))/alpha;
dpred(i)=mi(4)+ttt(i);
d(i)=dobs(i)-dpred(i);
end
G=[dtdx dtdy dtdz dtdt];
% Getting Eigenvalues and Eigenvectors of G
[U,S,V] = svd(G);
% Least square solution
dm=inv(G'*G)*G'*d;
% updating the 'mi' vector
mi=mi+dm;
mic(cont,:)=(mi(:))';
dc=G*dm;
%E=sum((dobs-dpred).^2);
E=(sum(d(:)).^2);
% Plot Each Iterative Solution as a blue circle
%plot3(mi(cont2,1),mi(cont2,2),-mi(cont2,3),'bo')
%hold on
end
%plot3(mi(1),mi(2),-mi(3),'ro');
end
%deleting zeros rows in the mic matrix
av=mic;
av(~any(mic,2),:)=[];
mic=av;
clear av
for i=2:cont-1
% Plot each iteratively calculated solution as a blue circle
plot3(mic(i,1),mic(i,2),-mic(i,3),'o','Color',[0 0 1]);
end
sprintf(' # of iterations %d',cont)
% Plot the best model solution as a red circle
hold on
plot3(mic(cont,1),mic(cont,2),-mic(cont,3),'o','Color',[1 0 0]);
xlabel('X'),ylabel('Y'), grid;
hold on
plot3(xco(:,:),yco(:,:),-zco(:,:),'*','Color',[0 0 1])
(
geocities.com/mx)