      program sp840     ! 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/7654321/

      real*8 x, y                          ! グラフ描画のための出力変数
      real*8 f, t                          !
      f(t) = 1d0/lambda * exp(-t/lambda)   ! 確率密度関数の文関数

      real*8 detection_Time                ! 観測時間(sec)

          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)

      open(1, file='sp810_detection_Time.dat') ! <-- detection_Time
      read(1,*) detection_Time
      close(1)
!     write(*,*) detection_Time

!     指数乱数によるシミュレーション
!     simulationする観測時間 sum_simu が detection_Timeになるまで

!     do i=1,Ns_simu-1
      i = 0
      do while ( sum_simu .le. detection_Time )
        i = i + 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

      Ns_simu = i
      write(*,'(/a9,i10)')'N_sampl =',Ns_simu

      write(*,'(a9,f17.6,a4)')'sum_sim2=',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_sim2=',ave_simu,' sec', iCNT

      open(4,file='sp840.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='sp845.dat')                ! 1時間毎の飛来数
      do j=1, iCNT
        x = j + 0.5
        write(5,*) x, CNTperHOUR(j)
      enddo
      close(5)
      
      stop
      end