backward monte carlo simulator of one dimensional radiative transfer

bmc_1d_rad_trans.rb

#!/usr/bin/env ruby
# -*- coding: utf-8 -*-

# backward monte carlo simulator of one dimensional radiative transfer
# detector location is the top of atmosphere

np   = 10000 # number of photon
tau  = 1.0   # optical thickness
alb  = 0.5   # surface albedo
omg  = 0.5   # single scattering albedo
zu   = 10000 # upper boundary (m)
srad = 0.0   # summation of radiance
obsh = 10000 # observation height (upward radiance observation)

ext  = tau / zu # extinction coefficient (/m)

for i in 1..np do
  w  = 1.0
  dir = -1.0
  z = obsh
  while z <= zu do
    ftau = - Math.log(rand) # free optical thickness
    znew = z + dir * ftau / ext
    if znew >= zu
      break
    elsif znew <= 0.0
      znew = 0.0
      w = w * alb
      srad = srad + w * Math.exp(- ext * (zu - znew))
      dir = Math.sqrt(rand)
    else
      w = w * omg
      srad = srad + w * Math.exp(- ext * (zu - znew)) / 4.0
      dir = 1.0 - 2.0 * rand
    end
    if (w < 0.25)
      if (w > 0.5 * rand)
        w = 0.5
      else
        w = 0.0
        break
      end
    end
    z = znew
  end
end

puts "Radiance: #{srad / np}"