sample_standard_normal Subroutine

public subroutine sample_standard_normal(z, sigma_cutoff, rng_seed, rng_counter)

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

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(out) :: z(:,:)
real(kind=dp), intent(in), optional :: sigma_cutoff
integer(kind=i32), intent(in), optional :: rng_seed
integer(kind=i64), intent(in), optional :: rng_counter(:)

Calls

proc~~sample_standard_normal~~CallsGraph proc~sample_standard_normal sample_standard_normal interface~random_uniform random_uniform proc~sample_standard_normal->interface~random_uniform none~next_uniform counter_stream_type%next_uniform proc~sample_standard_normal->none~next_uniform pcg64_random_double53 pcg64_random_double53 none~next_uniform->pcg64_random_double53 proc~open_uniform_from_double53 open_uniform_from_double53 none~next_uniform->proc~open_uniform_from_double53

Called by

proc~~sample_standard_normal~~CalledByGraph proc~sample_standard_normal sample_standard_normal proc~sample_photo_emission_velocity sample_photo_emission_velocity proc~sample_photo_emission_velocity->proc~sample_standard_normal proc~sample_reservoir_face_particle_counter sample_reservoir_face_particle_counter proc~sample_reservoir_face_particle_counter->proc~sample_standard_normal proc~sample_shifted_maxwell_velocities sample_shifted_maxwell_velocities proc~sample_shifted_maxwell_velocities->proc~sample_standard_normal proc~sample_photo_ray_slot sample_photo_ray_slot proc~sample_photo_ray_slot->proc~sample_photo_emission_velocity

Source Code

  subroutine sample_standard_normal(z, sigma_cutoff, rng_seed, rng_counter)
    real(dp), intent(out) :: z(:, :)
    real(dp), intent(in), optional :: sigma_cutoff
    integer(i32), intent(in), optional :: rng_seed
    integer(i64), intent(in), optional :: rng_counter(:)
    type(counter_stream_type) :: stream
    integer :: n_total, i, row, column
    real(dp) :: r, theta, pi, u1, u2, z1, z2, cutoff

    ! Counter draws 100, 101, ... feed the Box-Muller pairs in order.
    if (present(rng_counter)) call stream%open(rng_seed, rng_counter, 100_i64)
    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)
      if (present(rng_counter)) then
        u1 = stream%next_uniform()
        u2 = stream%next_uniform()
      else
        call random_uniform(u1)
        call random_uniform(u2)
      end if
      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
    call stream%close()
  end subroutine sample_standard_normal