program sp820 ! AS data for 20000 events implicit none ! 指数乱数を使ったシミュレーション integer i, j integer Ntmax, Ntmax1 parameter( Ntmax=150, Ntmax1=Ntmax+1 ) real*8 tmin, tmax, dt, box(Ntmax1) parameter( tmin=0, tmax=300d0 ) real*8 lambda ! 平均自由行程(飛来平均時間差) integer Ns_simu ! サンプリング数 parameter( Ns_simu = 200000 ) ! 20000回 integer cnt_simu( Ntmax ) data cnt_simu / Ntmax*0 / real*8 ds12_simu ! 時間差 real*8 sum_simu, ave_simu ! 時間差の積算値, 時間差の平均値 data sum_simu/ 0d0 / integer day, hou, min, sec ! 観測時間の日、時、分、秒 real*8 deltaT ! 時間幅 parameter( deltaT=3600 ) ! 3600=1時間 integer iCNT, CNTperHOUR(10000) ! 1時間毎の飛来数 data CNTperHOUR/10000*0/ real*8 r ! 一様乱数 integer iy data iy/1234567/ real*8 x, y ! グラフ描画のための出力変数 real*8 f, t ! 確率密度関数の文関数 f(t) = 1d0/lambda * exp(-t/lambda) ! dt = ( tmax - tmin ) / Ntmax box(1) = tmin do j=1,Ntmax box(j+1) = box(j) + dt enddo open(1, file='sp810_lambda.dat') ! <-- lambda read(1,*) lambda close(1) ! lambda = ave ! 指数乱数によるシミュレーション do i=1,Ns_simu-1 call ransu(iy,r) ! 一様乱数 ds12_simu = - lambda * log(r) ! 指数乱数 -> 飛来時間差 do j=1,Ntmax ! 飛来時間差の度数分布 if( box(j).le.ds12_simu .and. ds12_simu.lt.box(j+1) ) then cnt_simu(j) = cnt_simu(j) + 1 go to 2000 endif enddo 2000 continue sum_simu = sum_simu + ds12_simu ! 観測経過時間(sec) iCNT = sum_simu / deltaT + 1 ! 1時間毎の飛来数 CNTperHOUR(iCNT) = CNTperHOUR(iCNT) + 1 enddo write(*,'(/a9,f12.1,a4)')'sum_simu=',sum_simu,' sec' day = sum_simu/(24*60*60) ! 日 hou = (sum_simu-day*24*60*60)/(60*60) ! 時 min = ((sum_simu-day*24*60*60)-hou*60*60)/60 ! 分 sec = ((sum_simu-day*24*60*60)-hou*60*60)-min*60 ! 秒 write(*,'(4x,4(i3,a1))') day,'d',hou,'h',min,'m',sec,'s' ! 観測時間 ave_simu= sum_simu / (Ns_simu-1) ! 時間差の平均 write(*,'(a9,1pe12.5,a4,i10)') 'ave_simu=',ave_simu,' sec', iCNT open(4,file='sp820.dat') ! 飛来時間差(sec)の分布 lambda = ave_simu do j=1,Ntmax x = box(j) + dt/2d0 y = float(cnt_simu(j))/dt/(Ns_simu-1) ! 飛来時間差の確率密度の分布 write(4,*) x, y, f(x) ! 期待曲線 enddo close(4) open(5, file='sp825.dat') ! 1時間毎の飛来数 do j=1, iCNT x = j + 0.5 write(5,*) x, CNTperHOUR(j) enddo close(5) stop end