x=1:1:120;
y=1:1:80;
for z1=1:1:80
for z2=1:1:120
if(z1<=(-z2*2+140))
z(z1,z2)=8*exp(-0.025*sqrt((z1-40)^2+(z2-40)^2));
end
if((z1>(-z2*2+140))&(z1<(-z2*2+220)))
z(z1,z2)=8*exp(-0.025*sqrt((z1-40)^2+(z2-40)^2))*(0.5-0.0125*(z1+2*z2-180))+(3.2425+0.00215*z2-0.0045*z1+0.0008*z1*z2)*(0.5+0.0125*(z1+2*z2-180));
end
if(z1>=(-z2*2+220))
z(z1,z2)=3.2425+0.00215*z2-0.0045*z1+0.0008*z1*z2;
end
end
end
subplot(1,2,1);
[f,g]=contour(x,y,z);
clabel(f,g);
%其中13号点是已知基准点
Point=[1 15 35 z(15,35)+0.05*rand(1)
2 15 65 z(15,65)+0.05*rand(1)
3 15 95 z(15,95)+0.05*rand(1)
4 35 20 z(35,20)+0.05*rand(1)
5 35 50 z(35,50)+0.05*rand(1)
6 35 80 z(35,80)+0.05*rand(1)
7 35 110 z(35,110)+0.05*rand(1)
8 55 20 z(55,20)+0.05*rand(1)
9 55 50 z(55,50)+0.05*rand(1)
10 55 80 z(55,80)+0.05*rand(1)
11 55 100 z(55,110)+0.05*rand(1)
12 75 35 z(75,35)+0.05*rand(1)
13 75 65 z(75,65)+0.05*rand(1)];
LINE=[1 2
2 3
1 4
1 5
2 6
3 7
4 5
4 5
5 6
6 7
4 8
5 9
6 10
7 11
8 9
9 10
10 11
8 12
9 12
9 13
10 13];
for i=1:1:21
A_0(i,LINE(i,1))=-1;
A_0(i,LINE(i,2))=1;
end
for i=1:1:21
for j=1:1:12
A(i,j)=A_0(i,j);
end
end
B=A;
for i=1:1:12
for j=1
Q(i,j)=((Point(i,2)-Point(j,2))^2+(Point(i,3)-Point(j,3))^2)^0.5;
end
end
for i=1:1:21
h1(i,1)=P(LINE(i,2),4)-Point(LINE(i,1),4);
end
for i=1:1:21
h2(i,1)=0;
if(LINE(i,1)==13)
h2(i,1)=2*Point(13,4);
end
if(LINE(i,2)==13)
h2(i,2)=-2*Point(13,4);
end
end
for i=1:1:21
l(i,1)=h1(i,1)+h2(i,1);
end
P=eye(21);
N11=A**P*A;
N12=A**P*(B*Q);
N22=(B*Q)**P*(B*Q);
N21=N12*;
Wx=A**P*l;
Walpha=(B*Q)**P*l;
M=N22-N21*N11*N12;
Walpha_2=Walpha-N21*(N11^(-1))*Wx;
alpha=M^(-1)*Walpha_2;
X=N11^(-1)*Wx-N11^(-1)*N12*alpha;
for z1=1:1:80
for z2=1:1:120
z_after(z1,z2)=0;
for i=1
z_after(z1,z2)=z_after(z1,z2)+alpha(i,1)*((((z1-Point(i,2))^2+(z2-Point(i,3))^2))^0.5);
end
end
end
subplot(1,2,2);
[f,g]=contour(x,y,z_after);
clabel(f,g);