C   Test program for chkspk (Check Spike)
C     read data from 'sampledata.txt'
C      data number : num ( = 18000 ) data
C     if datum x(i) seems wrong, chk(i) set to be 1, else 0
C


      program testdespike
      parameter( DATANUM=100000 )

      integer i, num, chk(DATANUM), memnum, sum;
      real x(DATANUM), chkfact

C      num memnum chkfact should be set by user

      num = 18000
      memnum = 1000
      chkfact = 3

      open(unit=10, file='sampledata.txt')
      do i=1, num
       read(10, *) x(i)
      end do
      close(10)

      call chkspk(num, memnum, chkfact, x, chk )

      sum = 0
      do i=1,num
       sum = sum+chk(i)
      end do

      write(*,*) "Error:", sum, "OK:", num-sum, "Total", num

      do i=1, num
      write(*,*) i, chk(i), x(i)
      end do

      stop
      end
       





      subroutine chkspk( num, memnum, chkfact, x, chk )
C
C   From  Jorgen Hojstrup, 'A statistical data screening procedure',
C   Measment Science and Technology, 4, 1530157, 1993, 
C   doi:10.1088/0957-0233/4/2/003
C
C  According to the above article, data quority is checked
C  using the statistics of 'memnum' data before checking data.
C  The criteria of the checking is follows
C
C	diff > chkfact * std  --> chk(i) = 1 : BAD DATA
C	  other               --> chk(i) = 0 : GOOD DATA
C
C   where  std: standard deviation
C          avr: average
C          var: variance
C          cor: correlation between x(i) and x(i-1)
C     
C
      integer num, memnum
      real chkfact, x(num)
      integer chk(num)


      integer i, i1, j
      real fcst, diff
      double precision avr, cor, var, pcor

C at first, set first memnum data are OK

      do i=1, memnum
       chk(i) = 0
      end do

      do i=2, num

        if( i .GT. memnum ) then
	  call statis( i-memnum-1, memnum, num,  x, chk, avr, cor, var)
	else
	  call statis( 0, memnum, num, x, chk, avr, cor, var)
	end if

	if( chk(i-1) .EQ. 0 ) then
	   fcst = cor * x(i-1) + (1.0 - cor) * avr
	else
	   j = i-1
	   do while( (j.GT.0) .AND. (chk(j).NE.0 ) )
	     j=j-1
	   end do
	     i1 = i - j
	     pcor = cor**i1
	     fcst = pcor * x(j) + (1.0-pcor)*avr
        end if
	diff = abs( fcst - x(i) )
	if( diff .GT. (chkfact*sqrt(var) ) ) then
		chk(i) = 1
	else
		chk(i) = 0
        end if

      end do
      end


      subroutine statis( i1, n, tnum, x, chk, avr, cor, var )
C
C   calc. avr cor variance  i=i1 --  i<i1+n .
C        cor :   correlation betweem x[i] and x[i+1] 
C	        calculation are made only chk[i] == 0 */
C
C
C
      integer i1, n, tnum
      real x(tnum)
      integer chk(tnum)
      double precision avr, cor, var

      integer i, num
      real t, u
      double precision s1, s2, s3

      s1 = 0
      s2 = 0
      s3 = 0
      num = 0

      do i=i1, i1+n-1
         if( chk(i) .EQ. 0 ) then
	   s1 = s1 + x(i)
	   num = num + 1
	 end if
      end do
      avr = s1/num

      num = 0
      do i=i1, i1+n-1
         if( chk(i) .EQ. 0 ) then
	   t = x(i) - avr
	   s2 = s2 + t*t
	   num = num +1
	 end if
      end do
      
      var = s2/(num-1)

      num = 0
      do i=i1+1, i1+n-1
         if( (chk(i) .EQ. 0).AND.(chk(i-1) .EQ. 0) ) then
	   t = x(i) - avr
	   u = x(i-1) - avr
	   s3 = s3 + t*u
	   num = num +1
	 end if
      end do

      cor = s3/(num-2)
      cor = cor/var

      end
      
