关于圆内均匀分布模拟的实现.txt

这是我在网上收集的用于计算圆内均匀随机数的一些方法

下面这个是我写的(已修正其中考虑不周之处），请大家测试：

Module Random2d
Implicit None 

    Type Point
      Real(8) x,y
    End Type

Contains
  Subroutine Random_Circle(Num,R,Pnt)
    Integer,Intent(in)::Num
    Real(4),Intent(in)::R
    Type(Point),Dimension(Num),Intent(Out)::Pnt
    Real(4) x,y
    Integer i
      Do i=1,Num
15        CALL RANDOM(x)
        CALL RANDOM(y)
        If (x*(1-x)+y*(1-y)<0.25) Goto 15
        Pnt(i).x=R*x
        Pnt(i).y=R*y
      End Do
      Return
   EndSubroutine
End Module

Program Main
  Use Random2d
  Implicit None
  Integer Num
  Parameter(Num=10000)
  Real(4) R
  Integer i
  Type(Point),Dimension(Num)::Pnt
    Print*, "Please Input the radius:"
    Read(*,*) R
    Call Random_Circle(Num,R,Pnt)
    Open(1,File='Rd2d_Circle.dat',Status='Unknown')
      Do i=1,Num
        Write(1,'(2(f12.6,2x))') Pnt(i).x,Pnt(i).y
      EndDo
    Close(1)
End Program Main






不懂就不要死撑，我没有说你的方法不能实现，只是一直在强调，你的方法远非最佳！
看看这个吧。

! this module is to the date type Operating System independent
module TypeKind
  implicit none
  integer(kind = kind(1)), parameter:: IP = kind(1)
  integer(kind = IP), parameter:: LP = kind(.true.)
  integer(kind = IP), parameter:: SP = kind(1.0) 
  integer(kind = IP), parameter:: DP = kind(1.d0)
end module TypeKind



! this subroutine generate Circle-Uniform-Distribution
! sample efficient 100%
subroutine CircleUniform(NO, radii, r)
  use TypeKind, WP => DP
  implicit none
  ! the number of point you want to simulation
  integer(kind = IP), intent(in):: NO
  ! the radii of the cycle
  real(kind = WP), intent(in):: radii   
  real(kind = WP), intent(out):: r(NO)
  
  integer(kind = IP):: i   ! for loop
  real(kind = WP):: r1(NO)
  
  call Random_Number( r1 )
  call Random_Number( r )

  forall( i = 1 : NO : 1 )
    r(i) = max( r(i), r1(i) )
  end forall
  r = radii * r

  return
end subroutine


program main
  use TypeKind, WP => DP
  implicit none
  ! the number of point you want to simulation
  integer(kind = IP), parameter:: NO = 10000
  ! the radii of the cycle
  real(kind = WP):: radii
  real(kind = WP):: r(NO)
  real(kind = WP):: theta(NO)
  integer(kind = IP):: i   ! for loop

  ! parameter
  real(kind = WP), parameter:: Pi = 3.14159265358979_WP
  
  call Random_Seed()
  radii = 1._WP
  call CircleUniform(NO, radii, r)

  ! 至此，已经得到圆均匀分布，但为了画出图形，还是把极角也
  ! 用 Monte-Carlo 撒点吧。
  call Random_Number( theta )
  theta = 2 * theta * Pi


  open(5, file = "Out06.dat")
  do i = 1, NO, 1
    !write(5, *) r(i), theta(i)
    ! 若是 不是去画图，这些完全没有必要。
    write(5, *) r(i) * cos(theta(i)), r(i) * sin(theta(i))
  end do
  close(5)

  stop
end program main






! this module is to the date type Operating System independent
module TypeKind
  implicit none
  integer(kind = kind(1)), parameter:: IP = kind(1)
  integer(kind = IP), parameter:: LP = kind(.true.)
  integer(kind = IP), parameter:: SP = kind(1.0) 
  integer(kind = IP), parameter:: DP = kind(1.d0)
end module TypeKind



! this subroutine generate Circle-Uniform-Distribution
! sample efficient 100%
subroutine CircleUniform(NO, radii, r)
  use TypeKind, WP => DP
  implicit none
  ! the number of point you want to simulation
  integer(kind = IP), intent(in):: NO
  ! the radii of the cycle
  real(kind = WP), intent(in):: radii   
  real(kind = WP), intent(out):: r(NO)
  
  integer(kind = IP):: i   ! for loop
  real(kind = WP):: r1(NO)
  
  call Random_Number( r1 )
  call Random_Number( r )

  forall( i = 1 : NO : 1 )
    r(i) = max( r(i), r1(i) )
  end forall
  r = radii * r

  return
end subroutine


program main
  use TypeKind, WP => DP
  implicit none
  ! the number of point you want to simulation
  integer(kind = IP), parameter:: NO = 20000
  ! the radii of the cycle
  real(kind = WP):: radii
  real(kind = WP):: r(NO)
  real(kind = WP):: theta(NO)
  integer(kind = IP):: i   ! for loop

  ! parameter
  real(kind = WP), parameter:: Pi = 3.14159265358979_WP
  Integer(2) itime(2,4),intet(4),ihr,imin,isec,i100th

  call gettim(ihr,imin,isec,i100th)    
        itime(1,1)=ihr
        itime(1,2)=imin
        itime(1,3)=isec
        itime(1,4)=i100th
  
  call Random_Seed()
  radii = 1._WP
  Do i=1,No
    call CircleUniform(NO, radii, r)
  ! 至此，已经得到圆均匀分布，但为了画出图形，还是把极角也
  ! 用 Monte-Carlo 撒点吧。
    call Random_Number( theta )
    theta = 2 * theta * Pi
  EndDo


!  open(5, file = "Out06.dat")
!  do i = 1, NO, 1
!    !write(5, *) r(i), theta(i)
!    ! 若是 不是去画图，这些完全没有必要。
!    write(5, *) r(i) * cos(theta(i)), r(i) * sin(theta(i))
!  end do
!  close(5)
    call gettim(ihr,imin,isec,i100th)    
    itime(2,1)=ihr
    itime(2,2)=imin
    itime(2,3)=isec
    itime(2,4)=i100th
    do 50 i=1,4
        if(itime(2,i).ge.itime(1,i))then
            intet(i)=itime(2,i)-itime(1,i)
        else
            intet(i-1)=intet(i-1)-1
            if(i.ne.4)then
                intet(i)=itime(2,i)+60-itime(1,i)
            else
                intet(i)=itime(2,i)+100-itime(1,i)
            endif
        endif
 50    continue 
    write(*,200) intet(1),':',intet(2),':',intet(3),':',intet(4) 

200    format(2x,3(i2,a1),i2)

  stop
end program main
以上是Asym的程序

以下是我的程序
Module Random2d
Implicit None 

    Type Point
      Real(8) x,y
    End Type

Contains
  Subroutine Random_Circle(Num,R,Pnt)
    Integer,Intent(in)::Num
    Real(4),Intent(in)::R
    Type(Point),Dimension(Num),Intent(Out)::Pnt
    Real(4) x,y
    Integer i
      y=0
      Do i=1,Num
        x=5
        Do While(x-x*x+y-y*y<0.25)
          CALL RANDOM(x)
          CALL RANDOM(y)
        EndDo
        Pnt(i).x=R*x
        Pnt(i).y=R*y
      Enddo
      Return
   EndSubroutine
End Module

Program Main
  Use Random2d
  Implicit None
  Integer Num
  Parameter(Num=20000)
  Real(4) R
  Integer i
  Type(Point),Dimension(Num)::Pnt
  Integer(2) itime(2,4),intet(4),ihr,imin,isec,i100th

  call gettim(ihr,imin,isec,i100th)    
        itime(1,1)=ihr
        itime(1,2)=imin
        itime(1,3)=isec
        itime(1,4)=i100th
!    Print*, "Please Input the radius:"
!    Read(*,*) R
     R=100
    Do i=1,Num 
      Call Random_Circle(Num,R,Pnt)
    Enddo
!    Open(1,File='Rd2d_Circle.dat',Status='Unknown')
!      Do i=1,Num
!        Write(1,'(2(f12.6,2x))') Pnt(i).x,Pnt(i).y
!      EndDo
!    Close(1)
    call gettim(ihr,imin,isec,i100th)    
    itime(2,1)=ihr
    itime(2,2)=imin
    itime(2,3)=isec
    itime(2,4)=i100th
    do 50 i=1,4
        if(itime(2,i).ge.itime(1,i))then
            intet(i)=itime(2,i)-itime(1,i)
        else
            intet(i-1)=intet(i-1)-1
            if(i.ne.4)then
                intet(i)=itime(2,i)+60-itime(1,i)
            else
                intet(i)=itime(2,i)+100-itime(1,i)
            endif
        endif
 50    continue 
    write(*,200) intet(1),':',intet(2),':',intet(3),':',intet(4) 

200    format(2x,3(i2,a1),i2)

End Program Main

测试结果是
我的43.71秒,
Asym的44.69秒。





比如我说，产生 100000 事例，看看在 r ＝ 0.5 圆内的几率有多少。
我的程序如下：

! this module is to the date type Operating System independent
module TypeKind
  implicit none
  integer(kind = kind(1)), parameter:: IP = kind(1)
  integer(kind = IP), parameter:: LP = kind(.true.)
  integer(kind = IP), parameter:: SP = kind(1.0) 
  integer(kind = IP), parameter:: DP = kind(1.d0)
end module TypeKind



Module CUMC
  use TypeKind, WP => DP
  implicit none

contains
  ! this subroutine generate Circle-Uniform-Distribution
  ! sample efficient 100%
  subroutine CircleUniform(NO, radii, r)
implicit none
! the number of point you want to simulation
integer(kind = IP), intent(in):: NO
! the radii of the cycle
real(kind = WP), intent(in):: radii   
real(kind = WP), intent(out):: r(NO)
  
integer(kind = IP):: i   ! for loop
real(kind = WP):: r1(NO)
  
call Random_Number( r1 )
call Random_Number( r )

forall( i = 1 : NO : 1 )
  r(i) = max( r(i), r1(i) )
end forall
r = radii * r

return
  end subroutine


  subroutine Probability(a, P)
    implicit none
real(kind = WP), intent(in):: a
real(kind = WP), intent(out):: P

    ! the number of point you want to simulation
    integer(kind = IP), parameter:: NO = 100000
    ! the radii of the cycle
    real(kind = WP):: radii
    real(kind = WP):: r(NO)
    real(kind = WP):: theta(NO)
    integer(kind = IP):: i   ! for loop
    integer(kind = IP):: COUNT
    
    radii = 1._WP
call CircleUniform(NO, radii, r)
COUNT = 0
    do i = 1, NO, 1
  if ( r(i) < a ) then
    COUNT = COUNT + 1
  end if
end do
P = 1._WP * COUNT / NO 
return
  end subroutine

end module
  

program main
  use CUMC
  implicit none
  real(kind = WP):: a
  real(kind = WP):: P
  a = 0.5_WP
  call Probability(a, P)
  write(*, *) "P = ", P

  stop
end program main