      program sp830                        ! 10分間毎の飛来数の推移
      implicit none                        ! ポアソン分布
      integer i, j

      character dummy*10
      integer   Ns, hh1, mm1, ss1, ms1,  hh2, mm2, ss2, ms2
      parameter(Ns=20000)
      real*8    s1, s2, ds12

      real*8 sum, shiftT, deltaT           ! sum;時間差を加算する変数
      parameter( shiftT=120, deltaT=300 )  ! 5分後から解析，10分間毎の飛来数
      data   sum/0d0/
      integer iCNT,CNTperHOUR(10000)       ! 1時間毎の数
      data         CNTperHOUR/10000*0/

      integer day, hou, min, sec           ! 全観測時間の日、時、分、秒
      real*8 average                       ! 飛来時間差の平均値(秒)
      
      integer NPmax, NPmax1                   ! 時間差分布を求めるために
      parameter( NPmax=1000, NPmax1=NPmax+1 )
      real*8 Pxmin, Pxmax, dxP, boxP(NPmax1)  ! Pxmin～Pxmax(個) NPmax分割
      parameter( Pxmin=0, Pxmax=1000 )        ! 0～300secまで150分割
      integer cntP(NPmax)
      data    cntP/NPmax*0/
      real*8  cnt
      
      real*8 poisson, lambda
      integer m

      real*8 x_for_plot, y_for_plot

      dxP      = ( Pxmax - Pxmin ) / NPmax
      boxP(1) = Pxmin
      do j = 1, NPmax
        boxP(j+1) = boxP(j) + dxP
      enddo

      open(1, file='070109.DAT')
      read(1,*) dummy
      read(1,*) dummy
      read(1,*) dummy

      read(1,'(3i2,i7)') hh1, mm1, ss1, ms1
        s1  = hh1*60*60 + mm1*60 + ss1 + ms1*1d-6

      do i = 1, Ns-1
        read(1,'(3i2,i7)') hh2, mm2, ss2, ms2
        s2   = hh2*60*60 + mm2*60 + ss2 + ms2*1d-6
        ds12 = s2 - s1
        if( ds12.lt.0 ) then
            ds12 = ds12 + 24*60*60                    ! 日付が変わる時の処理
        endif

        sum  =   sum + ds12                           ! 経過時間(sec)
        iCNT = ( sum - shiftT ) / deltaT + 1
        if( iCNT.ge.1 ) then
            CNTperHOUR(iCNT) = CNTperHOUR(iCNT) + 1   ! 時間毎の飛来数の推移
        endif

        hh1 = hh2
        mm1 = mm2
        ss1 = ss2
        ms1 = ms2
        s1  = s2
      enddo
      close(1)

      write(*,'(4x,a5,f17.6,a4)')'sum =',sum,' sec'  ! 全観測時間(sec)

      day =   sum/(24*60*60)                          ! 日
      hou =  (sum-day*24*60*60)/(60*60)               ! 時
      min = ((sum-day*24*60*60)-hou*60*60)/60         ! 分
      sec = ((sum-day*24*60*60)-hou*60*60)-min*60     ! 秒
      write(*,'(4x,4(i3,a1))') day,'d',hou,'h',min,'m',sec,'s'

      average = sum / (Ns-1)                          ! 平均時間差
      write(*,'(a9,1pe12.5,a4,i10)') 'average =',average,' sec', iCNT

      open(2, file='sp830.dat')
      do j = 1, iCNT
        x_for_plot = (j+0.5)*deltaT/3600
        y_for_plot = CNTperHOUR(j)
        write(2,*) x_for_plot, y_for_plot            ! 時間毎の飛来数の推移
      enddo
      close(2)

      do i=1, iCNT-1
        cnt = CNTperHOUR(i)
        do j = 1, NPmax                               ! 飛来数の度数分布
            if( boxP(j).le.cnt .and. cnt.lt.boxP(j+1) ) then
                cntP(j) = cntP(j) + 1
                exit
            endif
        enddo
      enddo

      sum = 0
      do i= 1, iCNT-1
        sum = sum + CNTperHOUR(i)
      enddo
      lambda = sum/(iCNT-1)                           ! 飛来数の平均値
      
      open(3, file='sp835.dat')
      do j = 1, NPmax
        x_for_plot = boxP(j) + dxP/2d0
        y_for_plot = float(cntP(j))/dxP/(iCNT-1)      ! 時間差の確率密度の分布
        m          = boxP(j)
        write(3,*) x_for_plot,   y_for_plot,   poisson(m,lambda)
      enddo
      close(3)

      stop
      end

      real*8 function poisson(m,lambda)               ! ポアソン分布関数
      implicit none                                   ! lambda 平均値
      integer m, i
      real*8 lambda, kaijyo
      if(m.gt.0) then
        kaijyo = 1d0
        do i=1, m
         kaijyo = kaijyo * i
        enddo
      else
        kaijyo = 1d0
      endif
      poisson = lambda**m / kaijyo * exp(-lambda)
      return
      end