Berikut adalah source code simulasi Monte Carlo pada kasus peluruhan :
Source code dalam bentuk .cpp dapat didownload disini.
//********************************************************************************
#include <iostream>
#include <ctime>
#include <cstdlib>
#include <cmath>
#include <fstream>
using namespace std;
int main()
{
double N0 = 10000000; //inti induk awal;
double N0_an =N0; //untuk analitik
double lambda=0.2; //hitung konstanta peluruhan (konversi)
double p = 1-(exp(-lambda)); //peluang meluruh
double hl_step=0;
double rat_nN0=0; //peluang meluruh
double sum_hl=0; //half life wadah kosong
double hl_an=log(2)/lambda; //analitical half-life
int t=1; //untuk t=0 tanpa iterasi
srand( time(0) ); //avoid same generated number
int n = 0; // counter untuk menghitung inti yang meluruh
ofstream file_; //buat file
file_.open ("decay_out10000000.txt"); //buka, lalu tulis file
cout <<"#t(sec)"<<"\t"<<"#Undecayed (MC)"<<"\t"<<"#Undecayed (analitic)"<< "\n"; //saat t masih 0
file_<<t-1<<"\t"<<N0<< "\t"<< N0_an<< "\n"; //saat t masih 0
cout <<t-1<<"\t"<<N0<< "\t"<< N0_an<< "\n"; //saat t masih 0
Source code dalam bentuk .cpp dapat didownload disini.
while (N0>0) // hanaya run ketika : kondisi berhenti (according to both N0 or time)
{
for (int i = 0; i < N0; ++i)
{
double random=static_cast <double> (rand()) /( static_cast <double> (RAND_MAX)); //random keluaran
if (random<p){ //marking decayed isotop ("1" dice)
n = ++n; }
}
if (n>0){ //jika tak ada yang meluruh (n=0) random tidak dicount
rat_nN0=(n/N0); //rasio n/N0 untuk menghitung half life (estimasi MC)
hl_step=(log(0.5)/(log(1-(rat_nN0))));
sum_hl=sum_hl+hl_step;
N0=N0-n; //sisa inti induk MC
double N0_a=N0_an*(exp(-lambda*t)); //sisa inti induk (analitik)
cout <<t<<"\t"<<N0<<"\t"<<N0_a<<"\t"<<"\t"<<hl_step<<"\t"<<hl_an<<"\n";
file_<<t<<"\t"<<N0<<"\t"<<N0_a<<"\t"<<"\t"<<hl_step<<"\t"<<hl_an<<"\n";
if (N0<1){
break;
}
n=0; // set ulang counter
t=++t;
}
}
file_.close ();
double hl_MC=sum_hl/t;
cout<<"\n \n" ;
cout<<"--------------------------------\n";
cout<<"half-life (analitic)"<<"\t"<<"="<<hl_an<<"(s)\n";
cout<<"half-life (MC)"<<"\t \t"<<"="<<hl_MC<<"(s)\n";
file_.open ("hl_out10000000.txt");
file_<<hl_MC<<"\t"<<hl_an<<"\n";
file_.close ();
return 0;
}
//*******************************************************************************
Source code dalam bentuk .cpp dapat didownload disini.
Contoh hasilnya adalah :
Penjelasan lengkap dalam bentuk slide tersedia pada link youtube di sini.
Membutuhkan jasa pembuatan simulasi Matlab, C++ , dll ? Klik di sini
Simulasi Peluruhan Menggunakan Metode Acak Monte Carlo C++
May 12, 2017
1
ini C++ biasa mba, bukan code MCNP. Thanks
ReplyDelete