三站平面时差定位的MATLAB仿真和GDOP图

clc

clear

x2=-12.99x3=12.99xt=0x1=0%%%星型d=15km

y2=7.5y3=7.5yt=-15y1=0

z2=0.01z3=0zt=0.01 z1=0.01

%

% x1=-25.98x2=25.98x3=0xt=0%%%星型d=30km

% y1=15y2=15y3=-30yt=0

% z1=0.2z2=0.2z3=0.2 zt=0.25

% xt=-12.99x2=12.99x3=0x1=0%%%菱形型d=15km

% yt=7.5y2=7.5y3=15y1=0

% zt=0z2=0z3=0 z1=0.1

% x1=-25.98x2=25.98x3=0xt=0%%%菱形型d=30km

% y1=15y2=15y3=30yt=0

% z1=0z2=0z3=0 zt=0.1

% x1=-10.6x2=10.6x3=21.2xt=0%%%平行角型形型d=15km

% y1=10.6y2=10.6y3=0yt=0

% z1=0z2=0z3=0 zt=0.1

% x1=-21.2x2=21.2x3=42.4xt=0%%%平行型形型d=30km

% y1=21.2y2=21.2y3=0yt=0

% z1=0z2=0z3=0 zt=0.1

%

% x1=-30x2=0x3=30xt=0%%%倒三形型d=30km

% y1=30y2=30y3=30yt=0

% z1=0z2=0z3=0 zt=0.1

z=15

y=-200:1:200x=-200:1:200

%

for i=1:401

for j=1:401

% for i=-70:70

% for j=-70:70

m=x(i)

n=y(j)

r1=((m-x1).^2+(n-y1).^2+(z-z1).^2).^(1/2)

r2=((m-x2).^2+(n-y2).^2+(z-z2).^2).^(1/2)

r3=((m-x3).^2+(n-y3).^2+(z-z3).^2).^(1/2)

%r4=((m-x4).^2+(n-y4).^2+(z-z4).^2).^(1/2)

rt=((m-xt).^2+(n-yt).^2+(z-zt).^2).^(1/2)

c11=(m-x1)/r1c21=(m-x2)/r2c31=(m-x3)/r3ct1=(m-xt)/rt

c12=(n-y1)/r1c22=(n-y2)/r2c32=(n-y3)/r3ct2=(n-yt)/rt

c13=(z-z1)/r1c23=(z-z2)/r2c33=(z-z3)/r3ct3=(z-zt)/rt

c=[(-ct1+c11) (-ct2+c12) (-ct3+c13)(-ct1+c21) (-ct2+c22) (-ct3+c23)(-ct1+c31) (-ct2+c32) (-ct3+c33)]

b=inv(c'*c)*(c')

b11=b(1)

b12=b(4)

b13=b(7)

b21=b(2)

b22=b(5)

b23=b(8)

b31=b(3)

b32=b(6)

b33=b(9)

% b14=b(10)

% b24=b(11)

% b34=b(12)

sigmap=0.5%测量误差的标准差

sigmas=0.0075%测时误差

eta=0.35%Ri站距离和测量之间的相关系数

sigma11=sigmap.^2+sigmas.^2*2

sigma22=sigmap.^2+sigmas.^2*2

sigma33=sigmap.^2+sigmas.^2*2

sigma44=sigmap.^2+sigmas.^2*2

sigma12=eta*sigmap.^2+sigmas.^2

sigma13=eta*sigmap.^2+sigmas.^2

sigma14=eta*sigmap.^2+sigmas.^2

sigma21=eta*sigmap.^2+sigmas.^2

sigma23=eta*sigmap.^2+sigmas.^2

sigma24=eta*sigmap.^2+sigmas.^2

sigma31=eta*sigmap.^2+sigmas.^2

sigma32=eta*sigmap.^2+sigmas.^2

sigma34=eta*sigmap.^2+sigmas.^2

sigma41=eta*sigmap.^2+sigmas.^2

sigma42=eta*sigmap.^2+sigmas.^2

sigma43=eta*sigmap.^2+sigmas.^2

sigmax2=b11*b11*sigma11+b11*b12*sigma12+b11*b13*sigma13+b12*b11*sigma21+b12*b12*sigma22+b12*b13*sigma23+b13*b11*sigma31+b13*b12*sigma32+b13*b13*sigma33

sigmay2=b21*b21*sigma11+b21*b22*sigma12+b21*b23*sigma13+b22*b21*sigma21+b22*b22*sigma22+b22*b23*sigma23+b23*b21*sigma31+b23*b22*sigma32+b23*b23*sigma33

sigmaz2=b31*b31*sigma11+b31*b32*sigma12+b31*b33*sigma13+b32*b31*sigma21+b32*b32*sigma22+b32*b33*sigma23+b33*b31*sigma31+b33*b32*sigma32+b33*b33*sigma33

gdop(j,i)=(sigmax2+sigmay2+sigmaz2).^(1/2)

end

end

figure(1)

[c,handle]=contour(x,y,gdop,25)%[c,handle]=contour(gdop,25)

clabel(c,handle)

% hold on

% plot(0,0,0.1,'rpentagram','linewidth',2)

% text(0,0,1,'T')

% plot(-30,0,0,'rpentagram','linewidth',2)

% text(-31,0,0,'R1')

% plot(30,0,0,'rpentagram','linewidth',2)

% text(31,0,0,'R2')

% plot(0,-5,0,'rpentagram','linewidth',2)

% text(0,-6,0,'R3')

% plot(0,5,0,'rpentagram','linewidth',2)

% text(0,6,1,'R4')

% hold off

xlabel('x方向(单位:km)')

ylabel('y方向(单位:km)')

% title('GDOP图')

知道怎么用蒙特卡洛仿真求圆的面积吗？下面有一个pudn上的代码，使用蒙特卡洛方法求圆心在原点，半径为1的圆的面积。你可以用类似的方法解决你的问题，只要把条件改成落点同时在三个圆里面就行了

sita=0:0.01:2*pi

x=sin(sita)

y=cos(sita)% 计算半径为1的圆周上的点，以便作出圆周观察

m=0 % 在圆内在落点计数器

x1=2*rand(1000,1)-1% 产生均匀分布于[-1, +1]直接的两个独立随机数x1,y1

y1=2*rand(1000,1)-1

N=1000% 设置试验次数

for n=1:N % 循环进行重复试验并统计

p1=x1(1:n)

q1=y1(1:n)

if (x1(n)*x1(n)+y1(n)*y1(n))<1 % 计算落点到坐标原点的距离,判别落点是否在圆内

m=m+1 % 如果落入圆中,计数器加1

end

plot(p1,q1,'.',x,y,'-k',[-1 -1 1 1 -1],[-1 1 1 -1 -1],'-k')

axis equal% 坐标纵横比例相同

axis([-2 2 -2 2]) % 固定坐标范围

text(-1,-1.2,['试验总次数 n=',num2str(n)])% 显示试验结果

text(-1,-1.4,['落入圆中数 m=',num2str(m)])

text(-1,-1.6,['近似圆面积 S_c=',num2str(m/n*4)])

set(gcf,'DoubleBuffer','on') % 双缓冲避免作图闪烁

drawnow % 显示结果

end

int min_figure(int &min,int a[],int n,int i){ //因为返回是整形

min=a[0]

for(i=0i<ni++)

{

if(a[i]<min)//找最小的 不是if(a[i]<a[0])

min=a[i]

}

return min

}

---