sample_standard_normal Subroutine

public subroutine sample_standard_normal(z, sigma_cutoff)

Box–Muller法で標準正規乱数を生成し、任意形状配列へ詰める。

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(out) :: z(:,:)
real(kind=dp), intent(in), optional :: sigma_cutoff

Called by

proc~~sample_standard_normal~~CalledByGraph proc~sample_standard_normal sample_standard_normal proc~sample_photo_raycast_particles sample_photo_raycast_particles proc~sample_photo_raycast_particles->proc~sample_standard_normal proc~sample_shifted_maxwell_velocities sample_shifted_maxwell_velocities proc~sample_shifted_maxwell_velocities->proc~sample_standard_normal proc~init_random_beam_particles init_random_beam_particles proc~init_random_beam_particles->proc~sample_shifted_maxwell_velocities proc~sample_reservoir_face_particles sample_reservoir_face_particles proc~sample_reservoir_face_particles->proc~sample_shifted_maxwell_velocities

Source Code

  subroutine sample_standard_normal(z, sigma_cutoff)
    real(dp), intent(out) :: z(:, :)
    real(dp), intent(in), optional :: sigma_cutoff
    integer :: n_total, i, row, column
    real(dp) :: r, theta, pi, u1, u2, z1, z2, cutoff

    n_total = size(z)
    pi = acos(-1.0_dp)
    cutoff = default_velocity_sigma_cutoff
    if (present(sigma_cutoff)) cutoff = sigma_cutoff
    if (cutoff <= 0.0_dp) error stop "sigma_cutoff must be > 0"

    i = 1
    row = 1
    column = 1
    do while (i <= n_total)
      call random_number(u1)
      call random_number(u2)
      if (u1 <= tiny(1.0_dp)) u1 = tiny(1.0_dp)
      r = sqrt(-2.0_dp*log(u1))
      theta = 2.0_dp*pi*u2
      z1 = r*cos(theta)
      z2 = r*sin(theta)
      if (abs(z1) <= cutoff) then
        z(row, column) = z1
        row = row + 1
        if (row > size(z, 1)) then
          row = 1
          column = column + 1
        end if
        i = i + 1
      end if
      if (i <= n_total .and. abs(z2) <= cutoff) then
        z(row, column) = z2
        row = row + 1
        if (row > size(z, 1)) then
          row = 1
          column = column + 1
        end if
        i = i + 1
      end if
    end do
  end subroutine sample_standard_normal