monte carlo simulation of one dimensional radiative transfer

mc_1d_rad_trans.rb

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

# monte carlo simulation of one dimensional radiative transfer
#
# The original code written in Fortran 77
# "Modeling of radiative transfer in cloudy atmospheres
#  and plant canopies using Monte Carlo methods", Tech. rep FTR, 7


np = 10000 # number of photon
tau = 1.0  # optical thickness
sr = 0.5   # surface reflectance
omg = 1.0  # single scattering albedo
zu = tau   # upper boundary
zb = 0.0   # lower boundary
sumr = 0.0 # summation of reflectance
sumt = 0.0 # summation of transmittance

for i in 1..np do
  w = 1.0
  z = zu
  dir = -1.0
  while z <= zu do
    path = - Math.log(rand)
    z = zu + dir * path
    if z <= zb then
      z = zb
      dir = 1.0
      sumt = sumt + w
      w = w * sr
    elsif z <= zu then
      if rand <= 0.5 then
        dir = 1.0 * dir
      else
        dir = -1.0 * dir
      end
      w = w * omg
    else
      sumr = sumr + w
    end
  end
end

puts "Reflectandce top:     #{sumr / np}"
puts "Transmittance bottom: #{sumt / np}"