最近碰到需要使用蒙地卡羅法找出圓錐的體積
示意圖:https://www.dropbox.com/s/gemf3q3yqrmg1lt/cone.jpg
想法: 利用圓柱來抓出1/3的圓錐，但投出機率是錯誤的...
以下是我的程式
nh=0
do i=1, n
call random_number (rand)
r=rand
call random_number (rand)
a=2*3.1415926*rand !角度一 0~2*pi
call random_number (rand)
b=2*3.1415926*rand !角度二 0~2*pi
call random_number (rand)
z=rand !高度 0~1
x=r*cos(a)
y=r*sin(b)
if (z>=f1(x,y)) then
nh=nh+1
endif
enddo
pr=(real(nh)/real(n)) !←主要是這裡的機率一直抓不到
area1=(real(nh)/real(n))*3.1415926*1**2
write(*,*) 'con area',area1
!external function
real function f1(x,y)
f1=sqrt(x**2+y**2)
end function
不曉得邊界的界定是錯在哪邊? 懇請大大幫忙~謝謝!