c234567
      program spillstat3
      common/pawc/hmemor(2000000)
      dimension icc(1),cont(5000),aamp(2000)
      character*6 aa

      call hlimit(2000000) 
      do ii =1,12 
      call hbook1(110+ii,' spill 5ms ',1000,1000.,6000.,0.) ! full range
      enddo
      call hbook1(101,' spill 1ms ',4000,1500.,5500.,0.)
      call hbook1(102,' spill 5ms ',800,1500.,5500.,0.)
      call hbook1(103,' spill 25ms ',160,1500.,5500.,0.)

      call hbook1(200,' mean 1 ms ',50,50.,150.,0.)
      call hbook1(201,' sd/sqrt(mean) 1 ms ',60,1.,4.,0.)
      call hbook1(202,' skewness 1 ms ',50,-1.,1.,0.)
      call hbook1(203,' kurtosis  1 ms ',40,1.,5.,0.)
    
      call hbook1(204,' lf noise ',25,0.,0.5,0.)
      call hbook1(205,' 100 Hz noise ',25,0.,0.25,0.)
      call hbook1(206,' total  noise ',30,0.,0.6,0.)


 
      call hbook2(301,' sd/sqrt(mean) vs skewness ',
     1 60,1.,4.,50,-1.,1.,0.)
      call hbook2(302,' sd/sqrt(mean) vs kurtosis ',
     1 60,1.,4.,40,1.,5.,0.)
      call hbook2(303,' skewness  vs kurtosis ',
     1 50,-1.,1.,40,1.,5.,0.)

      call hbook2(304,' sd/sqrt(mean) vs lf noise ',
     1 60,1.,4.,25,0.,0.5,0.)
      call hbook2(305,' sd/sqrt(mean) vs 100 Hz',
     1 60,1.,4.,25,0.,0.25,0.)
      call hbook2(306,' sd/sqrt(mean) vs total noise',
     1 60,1.,4.,30,0.,0.6,0.)


      do ii = 1,12
      call hbook1(8000+ii,' freq ',250,0.,250.,0.)
      call hbook1(9000+ii,' RMS ',50,0.,1000.,0.)
      call hbook1(9100+ii,' First Difference 1 ms ',50,-100.,150.,0.)
      call hbook2(10000+ii,' corrected folded spill, spill ',
     1 231,0.,924.,500,1100.,6100.,0.)
      call hbook1(11000+ii,' Corrected folded time',
     1 231,0.,924.,0.)
      call hbook1(12000+ii,' Uncorrected folded time',
     1 231,0.,924.,0.)
      enddo 

      do ii = 1,12
      icc(1) = ii
       if(icc(1).eq.1)
     1 open(unit=19,file='elis-12567-0336.txt',status='unknown') ! golden file 
       if(icc(1).eq.2)
     1 open(unit=19,file='elis-12567-0337.txt',status='unknown') ! golden file 
       if(icc(1).eq.3)
     1 open(unit=19,file='elis-12567-0338.txt',status='unknown') ! golden file 
       if(icc(1).eq.4)
     1 open(unit=19,file='elis-12567-0339.txt',status='unknown') ! golden file 
       if(icc(1).eq.5)
     1 open(unit=19,file='elis-12567-0340.txt',status='unknown') ! golden file 
       if(icc(1).eq.6)
     1 open(unit=19,file='elis-12567-1146.txt',status='unknown') ! golden file 
       if(icc(1).eq.7)
     1 open(unit=19,file='elis-2465-0170.txt',status='unknown') ! dataplot31.kumac
       if(icc(1).eq.8)
     1 open(unit=19,file='elis-2566-1230.txt',status='unknown') ! dataplot32.kumac
       if(icc(1).eq.9)
     1 open(unit=19,file='elis-2066-1129.txt',status='unknown') ! dataplot33.kumac
      if(icc(1).eq.10)
     1 open(unit=19,file='elis-12288-1360.txt',status='unknown') ! gold 22
      if(icc(1).eq.11)
     1 open(unit=19,file='elis-13019-0133.txt',status='unknown') !  new 23
      if(icc(1).eq.12)
     1 open(unit=19,file='elis-13046-0011.txt',status='unknown') ! new 23

      read(19,5510,END=999)aa,nrun,nbur
 5510 format(a6,2i8)
      write(*,*)' '
      write(*,*)' run burst ',nrun,nbur



      nnn =0
      nn = 500000 !  max input - see arrays
      do k =1,nn
      read(19,5513,END=999)aa,kk,kid,kcode,kstamp
      nnn = nnn+1
 5513 format(a6,2i8,2z10)
      fm = float(kstamp)/40000d0 ! ms
      call hfill(110+ii,fm,0.,1.) ! full spill 5 ms
      call hfill(101,fm,0.,1.)
      call hfill(102,fm,0.,1.)
      call hfill(103,fm,0.,1.)
c--------------------------------
c     folded data

      ffb = 923.9926  !  fold time
      tempx = mod(float(kstamp),ffb)  !  folded spill
      tempy = float(kstamp)/40000d0  ! ms spill
      xxx = (tempy -1500.)/1000.
      corr = 2.4*xxx**3
      ctempx = tempx - corr ! corrected folded spill
      if(ctempx.lt.0.)ctempx = ctempx +924.0
      call hfill(10000+ii,ctempx,tempy,1.)
      call hfill(11000+ii,ctempx,0.,1.)
      if(tempy.gt.1200..and.tempy.lt.2200.)then
c     uncorrected folded
      call hfill(12000+ii,tempx,0.,1.)
      endif


c--------------------------------
      enddo
 999  continue

      write(*,*)' count  ', nnn
      write(*,*)' '
      close (19)

c     RMS plots
      call hunpak(102,cont,' ',0) ! 5 ms interval 1.5 - 5.5
      do jj =1,800
      call hfill(9000+ii,cont(jj),0.,1.)
      enddo

c     1.5 - 5.5 sec
      call hunpak(101,cont,' ',0) ! 1 ms interval

c     first difference plots
      do i = 1,3999
      temp = cont(i+1) -cont(i) ! first diference of spill
      call hfill(9100+ii,temp,0.,1.)
      enddo

c     spill statistics fot 1 ms
      fnn = 4000
      smean = 0.
      do i = 1,fnn
      smean = smean +cont(i)
      enddo
      smean = smean/fnn
      sd = 0.
      skew = 0.
      skurt = 0.
      do i = 1,fnn
      sd = sd + (cont(i) - smean)**2
      skew = skew + (cont(i) - smean)**3
      skurt = skurt + (cont(i) - smean)**4
      enddo

      sd =sqrt(sd/fnn)
      skew = skew/sd**3/fnn
      skurt = skurt/sd**4/fnn
c     try excess RMS <<<<<<<<<<<<<<<<<<
      sd = sqrt(sd**2-n)  !  mod
      sdn = sd/sqrt(smean)

      write(*,*)' mean,sd/sqrt(mean),skewness,kurtosis ',
     1  smean,sdn,skew,skurt    
 
      call hfill(200,smean,0.,1.)
      call hfill(201,sdn,0.,1.)
      call hfill(202,skew,0.,1.)
      call hfill(203,skurt,0.,1.)

      call hfill(301,sdn,skew,1.)
      call hfill(302,sdn,skurt,1.)
      call hfill(303,skew,skurt,1.)
c-------------------------------------------


*     periodogram    
 
      ssm = 0.
      do i =1,1000
      ssm = ssm + cont(i)
      enddo
      ssm = ssm/1000.
*     ssm = 0. !    <<<<<<<<< note <<<<<<<<<<<<<
      do i = 1,1000
      cont(i) = cont(i) - ssm
      enddo
      pi2 = 6.283185
      pi = 3.141592654
c     form periodogram
      do iw = 0,999 ! freq
      xs = 0.
      xc =0.
      do i = 1,1000
      t = i
      t = t/1000.  ! ms
      fq = iw
      xs = xs + cont(i)*sin(pi2*fq*t)
      xc = xc + cont(i)*cos(pi2*fq*t)
      enddo
      aamp(iw+1) = (xs**2 + xc**2)
      aamp(iw+1) = sqrt(aamp(iw+1))
      temp =iw
      call hfill(8000+ii,temp,0.,aamp(iw+1))
      enddo


c     noise levels
c     background ~ 400.
      sn1 = 0.
      do j = 1,40
      sn1 = sn1 + aamp(j) - 400.
      enddo
      sn2 = 0.
      do j = 95,105
      sn2 = sn2 + aamp(j) - 400.
      enddo
      sn2 = abs(sn2)  ! avoid small negative numbers
      sn3 = 0.
      do j = 1,250
      sn3 = sn3 + aamp(j) !  total signal
      enddo
      sn4 = 0.
      do j = 1,250
      sn4 =  sn4 + aamp(j) -400.  ! total noise 
      enddo
     
c     for consistency - less sensitive to 400 subtraction
c     but can omit 50 Hz signal

c     sn4 = sn1 + sn2        

 
      write(*,*)' noise:  low ,100 Hz, total ',sn1,sn2,sn3        
      write(*,*)' noise fractions low, 100 Hz, total noise ',
     1  sn1/sn3, sn2/sn3, sn4/sn3

      call hfill(204,sn1/sn3,0.,1.)
      call hfill(205,sn2/sn3,0.,1.)
      call hfill(206,sn4/sn3,0.,1.)

      call hfill(304,sdn,sn1/sn3,1.)
      call hfill(305,sdn,sn2/sn3,1.)
      call hfill(306,sdn,sn4/sn3,1.)

c234567
      
c--------------------------------------------
*     call hreset(100,' ')
      call hreset(101,' ')
      call hreset(102,' ')
      call hreset(103,' ')

      enddo  ! loop over files

      call hprint(200)
      call hprint(201)
      call hprint(202)
      call hprint(203)
      call hprint(204)
      call hprint(205)
      call hprint(206)

      call hprint(301)
      call hprint(302)
      call hprint(303)
      call hprint(304)
      call hprint(305)
      call hprint(306)

*     periodogram 
      do ii =1,12
*     call hprint(8000+ii)
      enddo


      call hrput(0,'spillstat3.dat','n')
      stop
      end



