## KOMPUTASI DISTRIBUSI NEUTRON DALAM
STATISTIK MAXWELL BOLTZMANN
Tuti Purwoningsih (tuti@ut.ac.id) Program Studi Pendidikan Fisika FKIP-UT

## ABSTRAK
Perpindahan neutron diatur dengan beberapa peluang distribusi, seperti peluang distribusi sudut hamburan, peluang distribusi jarak perpindahan, peluang distribusi energi transfer, serta peluang distribusi fluks neutron. Salah satu aplikasi dari pola distribusi ini adalah pemodelan reaksi antara neutron dengan elemen-elemen penyusun jaringan yang terkait dengan serapan neutron dan dosis yang terserap oleh jaringan tumor otak pada terapi BNCT (Boron Neutron Capture Therapy). Dalam penelitian ini dibahas analisis komputasi tentang pola distribusi fluks neutron dalam suatu sistem reaktor. Variabel dalam penelitian ini adalah banyaknya neutron yang disimulasikan, serta kedalaman sistem reaktor yang dalam penelitian ini menggunakan sistem reaktor berbentuk silinder. Hasil simulasi menunjukkan bahwa dengan neutron sebanyak 100.000 diperlukan waktu eksekusi sekitar 20 – 120 menit untuk menghasilkan pola distribusi fluks neutron yang bergantung pada kedalaman sistem yang digunakan. Untuk semua kedalaman, puncak pola distribusi fluks neutron berada pada bin ke tiga. Dari perbandingan antara hasil simulasi dengan eksperimen dari literatur, dapat disimpulkan bahwa melalui analisis simulasi untuk distribusi fluks neutron diperoleh suatu distribusi Poisson yang mengikuti statistik Maxwell-Boltzmann.
Kata kunci: komputasi, Maxwell-Boltzmann, neutron

## ABSTRACT
The migration of neutron is arranged by some probability distributions such as probability of spread distribution, probability of distance distribution, probability of energy distribution and probability of flux distribution. One application of these pattern distributions is modelling the reaction between neutron and elements which compose the tissue related to the absorption of neutron in brain cancer tissues. This article explores computation analysis of pattern of distribution of neutron flux in a reactor system. Variables were the amount of neutron simulated and the depth of cylindrical reactor system. Simulations showed that 20-120 minutes was needed in executing 100,000 neutrons to build the distribution pattern of neutrons flux. This pattern was also depended on the depth of the system. In all depths, the peak of neutron flux distribution pattern was in the 3rd bin. Comparison between this simulations and experiment results in literatures showed that by analyzing the simulation of the distribution of neutron flux, a Poisson distribution which follows the Maxwell-Boltzmann was resulted.
Key words: computation, Maxwell-Boltzmann, neutron
Dalam ilmu mekanika statistik dibahas sistem-sistem makroskopik yang terdiri atas banyak partikel penyusun sistem tersebut. Sifat-sifat makroskopik sistem seperti suhu, energi dan sebagainya diturunkan dari kelakuan dan sifat-sifat partikel pendukungnya yang sangat banyak itu. Sifat dan kelakuan partikel-partikel penyusun sistem makroskopik secara individual disebut dengan keadaan mikroskopik dari sistem makroskopik tersebut. Banyaknya partikel dalam sistem yang bergerak bebas menyebabkan sulitnya mengetahui besaran fisis sifat-sifat makroskopik sistem. Namun, besaran-besaran tersebut dapat diramalkan keadaaan fisisnya berdasarkan peluangnya (probabilitas).
Salah satu prinsip yang dapat diterapkan dalam probabilitas adalah prinsip acak (random), yaitu semua partikel memiliki peluang yang sama untuk berada di titik manapun dalam sistem dan bisa bergerak ke arah manapun dengan momentum manapun sesuai dengan batasan-batasan yang berlaku. Di alam, gerak partikel yang mengalami stokastik adalah neutron. Perpindahan neutron diatur dengan peluang distribusi seperti peluang distribusi  sudut hamburan, peluang distribusi jarak perpindahan, peluang distribusi energi transfer, serta peluang distribusi fluks.
Banyak artikel ilmiah yang telah diterbitkan membahas metode Monte Carlo dalam bidang radiasi nuklir dari tahun 90an sampai sekarang (misalnya: Hendricks, 1994; Kumada, et.al, 2007), namun belum ada yang membahas komputasi neutron dalam mekanika statistik berdasarkan distribusi statistik Maxwell-Boltzmann. Pembahasan ini dianggap penting karena pembahasan tersebut merupakan dasar aplikasi mekanika statistik dalam radiasi nuklir, seperti dalam terapi pengobatan Boron Neutron Capture Therapy (BNCT).
Secara umum tujuan dari penulisan artikel ini adalah menentukan komputasi distribusi fluks neutron dalam statistik Maxwell-Boltzmann. Diharapkan dengan lebih memaknai konsep distribusi fluks neutron dalam statistik Maxwell-Boltzmann, dapat dikembangkan bidang-bidang fisika lain yang mengaplikasikan karakteristik fluks neutron.
Dalam mekanika statistik, banyak perumusan matematika tingkat tinggi, sehingga perlu pemahaman yang lebih mendalam tentang matematikanya sebelum memahami mekanika statistiknya. Mekanika statistik membahas sistem-sistem makroskopik yang terdiri atas sangat banyak partikel penyusun sistem tersebut, sehingga sangat sulit untuk menentukan besaran fisisnya secara pasti satu per satu partikel (Gentle, 2005).
Terdapat tiga jenis fungsi distribusi yang telah dikenal dalam mekanika statistik, yaitu: Maxwell-Boltzmann, Bose-Einstein, dan Fermi Dirac. Jika ‘partikel’nya masing-masing dapat dibedakan (bukan partikel identik), maka fungsi distribusinya mengikuti Statistik Maxwell-Boltzmann. Namun, jika yang dihitung adalah kumpulan partikel identik, maka ada dua kemungkinan yaitu jenis fungsi Boson jika distribusinya mengikuti Statistik Bose-Einstein dan jenis fungsi Fermion jika distribusinya mengikuti Statistik Fermi Dirac  (Cottingham & Greenwood, 2004).

## Distribusi Maxwell-Boltzmann
Dalam mekanika statistik, gerak partikel merupakan suatu individu yang bergerak tidak teratur (brownian).  Jika suatu sistem partikel dengan massa masing-masing adalah m 1 , m 2 , m 3 , ...., m N dan kecepatannya adalah v 1 , v 2 , v 3 , ..., v N . maka energi rata-rata suatu partikel sistem tersebut adalah sebagai berikut.
, = 1 ∑ 1 2 2 (1)
dengan N merupakan jumlah partikel total. Jika semua partikel memiliki massa yang sama maka
, = 1 1 2 2 = 1 2 1 2 = 1 2 2 (2)
di mana v rms disebut root-mean-square velocity dan didefinisikan sebagai berikut.
2 = 1 ( ∑ 2 ) = 1 ( 1 2 + 2 2 + 3 2 +   … ) (3)
Sistem yang sesuai untuk memeriksa hubungan energi internal total pada sebuah sistem dari N molekul pada suhu mutlak T adalah sampel gas ideal. Teori kinetik elementer dari gas
menunjukkan bahwa hanya jika energi kinetik molekuler rata-ratanya 3 2 kT hukum gas ideal memiliki bentuk PV = NkT sehingga energi molekuler total haruslah 3 2 E NkT  . Oleh karena energi kinetik rata-rata sebuah partikel gas adalah 2 , 1 2 K ave rms E mv  , maka dari persamaan gas ideal dapat disimpulkan bahwa energi rata-rata sebuah partikel gas berhubungan dengan suhu mutlak, hal ini ditunjukkan oleh persamaan berikut ini.
− = 3 2 (4)
dengan k adalah konstanta Boltzman dan T adalah suhu mutlak.
Persamaan (4) menunjukkan bahwa suhu mutlak suatu gas ideal berbanding lurus dengan energi kinetik rata-rata partikel-partikelnya sehingga satu Kelvin setara dengan 2,07  10-23 J atau 1,29  10-4 eV (per partikel). Jika energi kinetik rata-rata yang diberikan dalam persamaan (4) disubstitusikan ke persamaan (2) maka akan didapat v rms sebagai berikut.
= 3 (5)
Kecepatan dalam persamaan (5) bersesuaian dengan energi kinetik rata-rata partikel gas. Bentuk distribusi kecepatan dalam kasus ini adalah sama di segala arah (isotropik), sehingga komponen sudut (anguler) dapat diintegrasikan dengan batas-batasnya. Bentuk distribusi kecepatan tersebut sesuai dengan distribusi Maxwell-Boltzmann.
Dalam distribusi Maxwell-Boltzman, gejala-gejala termodinamika dipahami dengan pendekatan sebagai peluang terbesar dari banyak partikel. Hal ini jelas berbeda dengan teori kinetik yang mempersoalkan gerak sebenarnya dari suatu partikel. Ilmu mekanika statistik berusaha menentukan peluang terbesar bagaimana sejumlah energi tertentu terdistribusi di antara anggota sistem partikel dengan kesetimbangan termal pada suhu mutlak T. Partikel-partikelnya dianggap berinteraksi satu sama lain dan dengan dinding wadahnya sedemikian sehingga cukup untuk menimbulkan kesetimbangan termal, tetapi tidak terlalu besar sehingga geraknya terkorelasi secara kuat.
Molekul gas merupakan partikel identik yang  jaraknya cukup berjauhan sehingga satu sama lainnya dapat dibedakan, dan fungsi distribusi Maxwell-Boltzmann yang berlaku untuknya adalah
( ) = − / (6)
Harga A bergantung dari jumlah partikel dalam sistem, dan  yang memegang peranan serupa dengan konstanta normalisasi suatu fungsi gelombang adalah k yaitu konstanta Boltzmann yang besarnya sebagai berikut.
k = 1,381  10 -23 J/K = 8,616  10 -5 eV/K
Dalam teori kinetik elementer, kecepatan sebuah molekul dengan energi rata-rata 3/2 kT adalah sebagai berikut.
= ̅ 2 = 3 (7)
karena 2 3 1 2 2 mv kT  . Kecepatan ini diberi lambang
karena akar rata-rata kecepatan
molekuler kuadrat (the root-mean-square speed) tidak sama dengan rata-rata aritmatik sederhana  . Hubungan antara   dan bergantung dari hukum distribusi yang berlaku untuk kecepatan molekuler dalam sistem tertentu. Untuk distribusi Maxwell-Boltzmann,
sehingga kecepatan rms adalah sekitar 9 persen lebih besar dari kecepatan rata-rata aritmatik.
Gambar 1. Distribusi kecepatan Maxwell-Boltzmann
Distribusi kecepatan persamaan (6) tidak simetris, sehingga kecepatan berpeluang terbesar v p lebih kecil dari   atau . Untuk mendapatkan v p , samakan turunan n(v) terhadap v dengan nol, dan pecahan persamaan yang timbul. Sehingga besarnya kecepatan berpeluang terbesar dapat dituliskan menjadi persamaan berikut.
(8)
v

## Interaksi Neutron
Secara sederhana dapat didefinisikan, neutron merupakan partikel yang memiliki suatu beban muatan nol dan massa yang lebih besar sedikit dari proton, atau sering dikenal dengan partikel tidak bermuatan. Sebagai partikel yang netral, lintasan neutron tidak akan terpengaruh oleh elektron orbital yang negatif di sekitar inti maupun medan listrik yang disebabkan oleh inti yang bermuatan positif. Pada saat neutron bertumbukan dengan inti, neutron akan terserap atau terhambur ke arah lain. Jadi hanya neutron yang tidak sempat bertumbukan tetap melaju searah sumbu x. Hal tersebut menyebabkan intensitas neutron tak berhamburan semakin berkurang ketika berkas tersebut semakin dalam menembus. Pergerakan berkas neutron yang searah sumbu x sesuai dengan Gambar 2.
Gambar 2. Pergerakan berkas neutron searah sumbu x
Dalam inti diasumsikan bahwa terdapat N inti/cm 3 sehingga pada tebal yang sangat kecil (mendekati nol) yaitu dx terdapat N.dx inti per cm 2 . Jika tiap inti memiliki luas penampang σ cm 2 , maka bagian luas yang terhalangi oleh inti adalah N. σ.dx sehingga hubungan untuk intensitas berkas neutron dapat ditulis sebagai berikut.
I(x+dx) = (1 – N.σ.dx) I(x) (9)
Selanjutnya didefinisikan penampang lintang makroskopik sebagai berikut.
Σ = Nσ (10)
Dalam persamaan (10), σ memiliki satuan cm 2 /inti yang dikenal dengan penampang lintang mikroskopik sehingga penampang lintang makroskopik memiliki satuan cm -1 . Persamaan ini merupakan probabilitas sebuah neutron akan berhamburan pada ketebalan dx selanjutnya, setelah sebelumnya lolos tanpa hamburan hingga jarak x.
Persamaan (10) menjelaskan tentang peluang atau probabilitas. Probabilitas neutron yang dituliskan sebagai p(x)dx akan mengalami hamburan pada area dx. Probabilitas tersebut bermakna
x d x I(0) I(x) I(x+dx)
dua hal, pertama adalah probabilitas neutron lolos hingga mencapai dx, sedangkan yang kedua probabilitas neutron mengalami hamburan pada area dx. Oleh karena probabilitas neutron bertumbukan pada area dx tidak bergantung pada kejadian sebelumnya, maka probabilitas total p(x)dx diperoleh dengan mengalikan kedua probabilitas berdasarkan persamaan (10) yang dapat dituliskan dalam persamaan berikut.
p(x) dx = Σ exp ( - Σx)dx (11)
Dengan mengetahui probabilitas neutron tersebut, jarak rerata antara dua hamburan yang ditempuh neutron dapat dihitung. Jarak ini biasa dikenal dengan jarak bebas rerata (mean free path) dengan notasi λ yang secara matematis ditulis sebagai berikut.
= ( ) = ( − ) = 1 ∞ 0 ∞ 0
(12)
Neutron-neutron yang terdapat dalam berkas I(x) adalah neutron yang belum mengalami hamburan (fluks tak terhambur/ uncollided flux) yang semuanya melaju dengan arah x positif yang sama. Di pihak lain, neutron-neutron yang telah mengalami hamburan dapat memiliki sebarang arah. Berkas neutron I(x) dapat dituliskan sebagai hasil kali dari laju neutron v dengan satuan cm/s, dan rapat neutron tak terhambur dengan satuan neutron/cm 3 . Secara matematis dapat diperoleh.
( ) = . ′′′ ( ) (13)
Persamaan (12) merupakan bentuk lazim dari penulisan fluks dengan notasi  , sehingga dapat dituliskan.
( ) = . ′′′ ( ) (14)
Secara umum, penampang lintang yang memberikan probabilitas terjadinya tumbukan disebut dengan penampang lintang total yang dinotasikan dengan σ t . Ketika neutron menumbuk inti, maka ada dua kemungkinan kejadian yang dialami neutron yaitu terhambur (scattered) atau terserap (absorbed). Kecenderungan untuk terjadinya hamburan atau serapan diperoleh dengan membagi penampang lintang total dalam penampang lintang hamburan dan penampang lintang serapan yang dihubungkan dalam persamaan matematis berikut.
= + (15)
Dalam peristiwa tumbukan, adalah probabilitas neutron untuk terhambur oleh inti sedangkan adalah probabilitas neutron akan terserap oleh inti.
Reaksi hamburan antara neutron dan inti dapat terjadi melalui dua proses yaitu hamburan elastik dan inelastik, yang dinyatakan sebagai berikut.
= + ′ (16)
dengan menyatakan penampang lintang hamburan elastik, sedangkan menyatakan penampang lintang hamburan inelastik.
Pada hamburan elastik nilai momentum dan energi kinetik neutron sebelum dan setelah hamburan bernilai sama. Pada hamburan inelastik, neutron memberikan sebagian energi kinetiknya ke inti dan meninggalkan inti dalam keadaan tereksitasi sehingga terpenuhi kekekalan momentum (Smith, 2000). Namun nilai energi kinetik neutron sebelum dan setelah hamburan tidak sama. Setelah mendapat energi dari neutron, inti mengeluarkan energi eksitasinya dengan memancarkan satu atau beberapa sinar gamma bersamaan dengan terpancarnya kembali neutron dari inti.

## METODE
Program dikembangkan dalam bahasa MATLAB 7.7.0 (R2008B) menggunakan notebook Atom (TM ) CPU N280 dengan kecepatan 1,66 GHz dan RAM 1,99 GB. Langkah-langkah pemrograman yang digunakan dalam penelitian ini terdiri dari 4 tahap.
Tahap pertama adalah menentukan Keadaan Awal Variabel. Program yang digunakan dalam penelitian ini menggunakan geometri berbentuk silinder dengan jari-jari tetap, dan kedalaman berbeda-beda. Setiap neutron yang disimulasikan memiliki posisi yang random, sehingga posisinya adalah posisi rata-rata. Kecepatannya merupakan kecepatan rata-rata dari konfigurasi neutron yang ada. Kecepatan awal neutron menurut sumbu y dan sumbu z ditetapkan nol, sedangkan pada sumbu x adalah v init . Lebar dari sel dalam orde mean free path sehingga sel dapat dianggap homogen. Tahap kedua adalah menentukan posisi pertama neutron. Secara fisis, neutron yang bergerak bebas amat sulit diketahui posisinya satu per satu secara pasti. Sehingga didefinisikanlah sebuah jalan bebas rata-rata (mean free path) yaitu  panjang lintasan sebelum neutron bertumbukan, yang biasa disimbolkan dengan  . Untuk  yang kontinyu, probabilitas neutron secara komulatif akan berdistribusi secara seragam dalam interval [0,1] yang diwakili oleh suatu bilangan random R yang berdistribusi seragam [0,1]. Sehingga posisi neutron secara random dapat dituliskan dalam persamaan berikut.
1 x LxR  (17)
dengan L adalah panjang lintasan dan R 1 adalah bilangan random yang terdistribusi secara seragam pada [0,1].
Tahap pemrograman ketiga adalah menyeleksi partikel ke dalam sel. Dalam sebuah sel, setiap neutron mempunyai peluang untuk bertumbukan yang berbanding langsung dengan kelajuan relatifnya.
( , ) = | − | ∑ ∑ | − | − 1 =1 =1 (18)
dengan N c adalah jumlah partikel dalam suatu sel.
Tahap terakhir pemrograman adalah menghitung besaran–besaran fisis termodinamika jika terjadi tumbukan. Saat neutron mengalami hamburan elastik, energinya tetap namun mengubah arah gerak, sehingga kecepatan pusat massanya tidak berubah akibat terjadinya tumbukan. Sudut azimut  terdistribusi secara seragam antara 0 sampai 2  , sehingga rapat probabilitasnya menjadi:
φ = 2 π R (19)
dengan R adalah bilangan random seragam [0,1]
Kecepatan setelah tumbukan ditentukan menggunakan persamaan berikut.
∗ = ∗ + 1 2 ∗ (20)
∗ = ∗ − 1 2 ∗ (21)
Fluks neutron ditentukan menggunakan hubungan antara jumlah neutron n dengan distribusi kecepatan v dalam persamaan berikut.
= (22)

## HASIL DAN PEMBAHASAN
Penelitian awal terapi biasanya dilakukan dengan cara in vitro ataupun simulasi dengan bantuan sistem komputasi. Dosimetri radiasi berhubungan dengan penentuan jumlah, ruang, dan distribusi sementara dari energi yang diberikan dalam bahan yang mengalami radiasi ionisasi. Pemodelan geometri dari jaringan yang digunakan dalam penelitian ini adalah bentuk silinder, dengan hasil keluaran program berupa pola distribusi fluks neutron pada kedalaman silinder yang berbeda-beda, yakni: 10 mm, 30 mm, dan 50 mm dengan jari-jari silinder 50 mm. Kemudian silinder dibagi menjadi daerah-daerah yang dalam penelitian ini dinamakan jumlah sel, sedangkan jumlah sel dalam penelitian ini adalah 15 sel.
Jumlah neutron yang disimulasikan untuk setiap kondisi dalam penelitian ini adalah 100.000 neutron dengan  jumlah langkah sebanyak 50 langkah. Grafik yang dihasilkan dari program komputasi adalah dalam bentuk histogram berdasarkan pengembangan dari program komputasi yang dikembangkan oleh Suarga (2007) dalam Program Metode Stokastik.
1 2 3 4 5 6 7 8 9 10 0 0.5 1 1.5 2 2.5 3 3.5 x 10 4 Plot Fluks Partikel Fluks Partikel J umlah partik el
Gambar 3. Distribusi Fluks Neutron dengan kedalaman silinder 10 mm
1 2 3 4 5 6 7 8 9 10 0 0.5 1 1.5 2 2.5 3 x 10 4 Plot Fluks Partikel Fluks Partikel J u m la h p a rt ik e l Gambar 4 Distribusi Fluks Neutron dengan kedalaman silinder 30 mm 1 2 3 4 5 6 7 8 9 10 0 0.5 1 1.5 2 2.5 3 x 10 4 Plot Fluks Partikel Fluks Partikel J u m la h p a rt ik e l
Gambar 5. Distribusi Fluks Neutron dengan kedalaman silinder 50 mm
Waktu yang diperlukan untuk memperoleh pola distribusi tiap-tiap kedalaman berbeda-beda, makin dalam silinder yang disimulasikan maka waktu perhitungannya juga akan semakin lama. Waktu simulasi yang diperlukan untuk setiap kedalaman disajikan dalam Tabel 1.
Tabel 1. Waktu Simulasi Fluks Neutron berdasarkan Kedalaman
Kedalaman
(mm)
Fluks Neutron Max.  (10 9 neutron/s cm) Waktu simulasi (menit) 10 0,5 28,36414 30 1,5 79,18081 50 3,0 153,4194
Dari Gambar 3 sampai Gambar 5 terlihat bahwa pola kurva distribusi fluks neutron yang dihasilkan tidak mengikuti pola distribusi normal atau Gaussian tetapi distribusi Poisson. Hal ini tampak dari pola distribusi fluks neutron untuk setiap kedalaman tertinggi pada bin ke 3, lalu mulai menurun secara tajam di bin ke 6 menuju ke nol. Namun besarnya fluks neutron tidak pernah mencapai angka nol karena energi kinetik neutron tidak akan menjadi nol akibat masih adanya kecepatan yang diberikan oleh inti atom.
Pola distribusi fluks neutron pada Gambar 3 sampai dengan Gambar 5 telah signifikan dengan pola distribusi fluks neutron yang diteliti oleh Kumada et.al (2007) dari hasil eksperimen, seperti terlihat pada Gambar 6.
.
Gambar 6. Pola distribusi termal fluks Neutron yang dilakukan oleh Kumada et al. (2007)
Penyerapan (absorpsi) menggeser spektrum termal neutron naik dari distribusi Maxwell- Boltzmann, hal ini terlihat dari pola distribusi fluks neutron yang dihasilkan bertambah sesuai dengan ukuran penampang lintang mikroskopik absorpsinya. Dalam keadaan interaksi neutron yang sesungguhnya, kesetimbangan yang lengkap tidak pernah dicapai sebelum absorpsi neutron berlangsung.

## KESIMPULAN
Bentuk grafik pola distribusi fluks neutron dalam penelitian ini sesuai dengan pola distribusi Poisson yang sesuai dengan statistik Maxwell Boltzmann. Keterbatasan memori komputer dan banyaknya partikel yang disimulasikan menyebabkan waktu yang cukup lama untuk melakukan simulasi. Perbedaan nilai fluks neutron yang diperoleh dalam penelitian ini dengan literatur akibat adanya ketidakpastian statistik yang digunakan dalam pendekatan numerik.
Penelitian ini telah dapat memperkirakan karakteristik dari proses perjalanan neutron untuk mencapai ketebalan tertentu. Prinsip dari proses perjalanan neutron adalah neutron bergerak dengan lintasan yang lurus, dan keluar dari lintasan tersebut hanya ketika neutron bertumbukan dengan inti yang membuatnya terhambur ke lintasan dengan arah baru atau hilang karena terabsorpsi. Sejarah kehidupan neutron terdiri atas beberapa kali tumbukan hamburan yang akhirnya mengalami reaksi absorpsi yang menyebabkan neutron hilang atau mati.
Untuk menyempurnakan penelitian, maka disarankan: 1) Perlu dilakukan kajian komputasi yang lebih mendalam mengenai solusi dari pola distribusi fluks neutron agar optimalisasi algoritma dan kode lebih baik; 2) Penelitian ini dapat dilakukan dengan menggunakan program komputer yang lebih baik dan canggih; 3) Melakukan variasi geometri reaktor agar dapat diaplikasikan untuk berbagai kasus dengan geometri yang berbeda-beda.

## REFERENSI
Cottingham, W. N., & Greenwood, D.A. (2004). An introduction to nuclear physics, (2 nd ed), United Kingdom: Cambridge University Press. Gentle, J.E. (2005). Random number generation and Monte Carlo method, (2 nd ed), New York: Springer
Hendricks, J.S. (1994). A Monte Carlo code for particle transport. Journal Los Alamos Science, 22. Kumada, H., Yamamoto, K., Matsumura, A., Yamamoto, T., & Nakagawa, Y. (2007). Development of JCDS, a computational dosimetry system at JAEA for boron neutron capture therapy, First European Workshop on Monte Carlo Treatment Planning. Journal of Physics: Conference Series, 74 (2007) 012010.
Smith, F.A. (2000). Primer in applied radiation physics. London: World Scientific. Suarga. (2007). Fisika komputasi: Solusi problema fisika dengan matlab. Yogyakarta: Andi.