求一个双基站的gdop的matlab仿真程序。。。。。。。

知道怎么用蒙特卡洛仿真求圆的面积吗？下面有一个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

1 x'-- x的转置矩阵

x*x'的结果是矩阵或数值，由x的维数（x若为向量，维数就是其长度）决定

2 inv 求逆矩阵。

miu=1/(ad'*inv(R)*ad)

既然程序没有报错，那么(ad'*inv(R)*ad)就是一个数值，而

w = miu*inv(R)*ad = 1/(ad'*inv(R)*ad) * inv(R)*ad 是矩阵运算，

不能按照代数四则运算进行简单的“约分”败饥。

3 y(1,v)=w.'*a

v=v+1 这是为了记录向量y，v表液枯耐示当前值在向量中的位置。算完这一个后，位置向后移动，从v到（v+1），继续循环，算出的y值闹春就记录在（v+1）位置处

对于光学的东西不是太懂，能做到的只是帮助你更容易地理解程序每一句话，希望对你有帮助

最好将下面的程序复制粘贴到matlab的编辑窗口中再阅读

clear % 清除matlab工作内存

lam=500e-9 %定义一个常数λ： 500乘以10的-9次方

a=1e-3D=1 %定义常数a：10的-3次方；常数D=1

ym=3*lam*D/a %定义ym

ny=51 %定义取点个数 ny=51

ys=linspace(-ym,ym,ny) %定义一个数组ys，从-ym到孙胡ym之间取郑凯和ny个点

np=51 %定义取点个数 np=51

yp=linspace(0,a,np) %定义一个数组yp，从0到a之间取np个点

for i=1:ny%循环语句，一共循环51次

sinphi=ys(i)/D %计算得到 sinΦ，

alpha=2*pi*yp*sinphi/lam % 计算得到 α=2Δ/λ=πypys/λD（由于yp是一个向量，这里得到的α也是一个向量，包含51个数）

sumcos=sum(cos(alpha)) % 计算得到 Σcosα

sumsin=sum(sin(alpha)) % 计算得到 Σsinα

B(i,:)=(sumcos^2+sumsin^2)/np^2 % 计算得到光强 I=I0(Σcosα)2+(Σsinα)2，并存储为矩阵B的第i行

end

N=255

Br=(B/max(B))*N%计算得到光强分布矩阵

subplot(1,2,1) %将绘图区分成左右两部分，首先在左侧绘图区绘图

image(ym,ys,Br)%绘制光强分布图

colormap(gray(N))%更改图像模式为灰度图像

subplot(1,2,2) %然后选择在右侧绘图区绘图

plot(B,ys) %绘制光强喊盯分布曲线

---