* Last processed by NICE on 12-Jun-2000 15:53:00
* Customized for :  IEEE, LINUX, UNIX, MOTIF, F77
c     $Header: /scisoft/gildas/sou/nic/main/nic_spike.for,v 1.43 1998/02/03 10:08:26 neri Rel $
c     $Log: nic_spike.for,v $
c     Revision 1.43  1998/02/03 10:08:26  neri
c     *** empty log message ***
c
C Revision 1.42  97/12/18  15:34:45  15:34:45  broguier ()
C add CALPOW variables
C
C Revision 1.32  96/07/03  18:54:06  18:54:06  neri (Roberto Neri)
C *** empty log message ***
C
C Revision 1.31  96/03/25  13:00:06  13:00:06  broguier (Dominique Broguiere)
C remove messages
C
C Revision 1.30  96/03/25  11:46:16  11:46:16  broguier (Dominique Broguiere)
C don't restore data after despiking
C
C Revision 1.29  96/03/19  15:00:42  15:00:42  broguier (Dominique Broguiere)
C remove COUT_SPIKES
C
C Revision 1.27  96/02/19  17:02:55  17:02:55  neri (Roberto Neri)
C *** empty log message ***
C
C Revision 1.21  95/09/27  09:55:48  09:55:48  gildas (Gildas manager)
C Suppress a forgotten call to SIC_LONG
C (S.Guilloteau)
C
C Revision 1.20  95/05/30  18:37:12  18:37:12  broguier (Dominique Broguiere)
C use new variables ip_lst, ip_el, ip_signal, ip_scan_coord,...
C
C Revision 1.17  95/02/27  13:23:11  13:23:11  broguier (Dominique Broguiere)
C use GI4_TRIE and GR4_SORT in SORT_SPIKE
C
C Revision 1.16  95/02/17  08:51:39  08:51:39  gildas (Gildas manager)
C idem
C
C Revision 1.15  1995/02/17  08:50:30  gildas
C Replaced GR4_TRI with GR4_TRIE
C TF
C
C Revision 1.14  1995/02/13  18:28:04  gildas
C Use library version of GR4_TRI
C
C Revision 1.12  1995/01/31  11:13:32  broguier
C keep record number and data associated for each spike
C
C Revision 1.10  95/01/21  18:27:15  18:27:15  neri (Roberto Neri)
C despike_flag updated for toon (-1)
C
C Revision 1.9  95/01/12  12:48:46  12:48:46  broguier (Dominique Broguiere)
C removal of variable IADDRSUBSCAN
C
C Revision 1.8  94/12/24  02:56:48  02:56:48  gildas (Gildas manager)
C Minor syntax correction: replaced a .EQ. by a .EQV.
C TF
C
C Revision 1.6  94/12/22  09:34:25  09:34:25  broguier (Dominique Broguiere)
C .eq. .false. replaced by .not.
C
*
      SUBROUTINE DESPIKE(LINE,SELECTED_CHANNEL,ERROR)
C ----------------------------------------------------------------------
C 	Support routine for command
C	DESPIKE   /CHANNEL /FUNCTION /XRMS /XWINDOW /ITERATION
C Arguments :
C	LINE	C*(*)	Command line			Input
C	ERROR	L	Logical error flag		Output
C
C ----------------------------------------------------------------------
C
C
C  Spike memory
C  ------------
C
C  Memory is allocated to store :
C    - number of spikes by channel,
C    - for each spike, number of record and data value
C  This is used to restore data each time command DESPIKE is called without
C  interpolation.
C
C  IPSPIKE is the pointer to this area (relative to MEMORY area)
C  LENGTH_SPIKE is the number of spikes which can be stored in the area. By
C  default LENGTH_SPIKE is 100 but becomes larger if necessary.
C
C  Structure of Spike Memory :
C      number_of_spikes[NCHAN]
C      index_of_spike[LENGTH_SPIKE] for receiver 1
C      spike_data[LENGTH_SPIKE]     for receiver 1
C      index_of_spike[LENGTH_SPIKE] for receiver 2
C      spike_data[LENGTH_SPIKE]     for receiver 2
C      ...
C  Length of Spike Memory is : NCHAN + 2*LENGTH_SPIKE*NCHAN
C
C
      CHARACTER*(*) LINE
      INTEGER*4 SELECTED_CHANNEL(*)
      LOGICAL ERROR
*
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
*
      LOGICAL   SIC_PRESENT
*
      INTEGER*4 IP_DATAR
      INTEGER*4 IP_NSPIKE,IP_INDEX_SPIKE,IP_DATA_SPIKE
      INTEGER*4 ICHAN,IFUNC
      INTEGER*4 IWINRMS,POINTER
      INTEGER*4 I, J
      INTEGER*4 NSPIKE, NSPIKE_TOTAL
      INTEGER*4 ITER,ITEMAX,IER
      INTEGER*4 SIC_GETVM
*
      REAL*4    FUNWID,XRMS
      LOGICAL   INTERPOLATE_ON_SPIKE
      CHARACTER DO_BLANKING*30
*
* Test if a file has been opened
*
      IF (IADDRRECORD.EQ.0) THEN
         CALL MESSAGE(2,1,'DESPIKE','No input file opened')
         ERROR = .TRUE.
         RETURN
      ENDIF
*
* Option /CHANNEL
*
      IF (SIC_PRESENT(1,0)) THEN
         I = 1
         DO WHILE (SIC_PRESENT(1,I).AND.I.LE.NCHAN)
            CALL SIC_I4 (LINE,1,I,SELECTED_CHANNEL(I),.TRUE.,ERROR)
            IF (ERROR) RETURN
            IF (SELECTED_CHANNEL(I).GT.NCHAN) THEN
               CALL MESSAGE (2,1,'DESPIKE','Bad channel number ')
               ERROR = .TRUE.
               RETURN
            ENDIF
            I = I+1
         ENDDO
         DO J=I,NCHAN
            SELECTED_CHANNEL(J) = 0
         ENDDO
         DO J=1,I-1
            SPIKE_REMOVAL_FLAG(SELECTED_CHANNEL(J)) = 0
         ENDDO
      ELSE
         DO J=1,NCHAN
            SELECTED_CHANNEL(J) = J
            SPIKE_REMOVAL_FLAG(J) = 0
         ENDDO
      ENDIF
*
* Option /FUNCTION
*
      IF (SIC_PRESENT(2,0)) THEN
         CALL SIC_I4 (LINE,2,1,IFUNC,.TRUE.,ERROR)
         IF (ERROR) RETURN
         IF (IFUNC.EQ.1) THEN
            IF (.NOT.SIC_PRESENT(2,2)) THEN
               FUNWID = 2.68
            ELSE
               CALL SIC_R4 (LINE,2,2,FUNWID,.TRUE.,ERROR)
               IF (ERROR) RETURN
            ENDIF
         ENDIF
      ELSE
         IFUNC = 1                       ! default sin(x)/x
         FUNWID = 2.68
      ENDIF
      IF (IFUNC.GT.3) THEN
         CALL MESSAGE (2,1,'DESPIKE','Function not yet implemented')
         ERROR = .TRUE.
         RETURN
      ENDIF
*
* Option /XRMS
*
      IF (SIC_PRESENT(3,0)) THEN
         CALL SIC_R4 (LINE,3,1,XRMS,.TRUE.,ERROR)
         IF (ERROR) RETURN
      ELSE
         XRMS = 5.                       ! Five sigma should be ok
      ENDIF
*
* Option /XWINDOW
*
      IF (SIC_PRESENT(4,0)) THEN
         CALL SIC_I4 (LINE,4,1,IWINRMS,.TRUE.,ERROR)
         IF (ERROR) RETURN
      ELSE
         IWINRMS = 0                     ! default: all the scan
      ENDIF
*
* Option /ITERATION
*
      ITEMAX = 1
      IF (SIC_PRESENT(5,0)) THEN
         CALL SIC_I4 (LINE,5,1,ITEMAX,.TRUE.,ERROR)
         IF (ERROR) RETURN
      ENDIF
*
* Option /INTERPOLATE
*
      INTERPOLATE_ON_SPIKE = .FALSE.
      IF (SIC_PRESENT(6,0)) THEN
         INTERPOLATE_ON_SPIKE = .TRUE.
      ENDIF
*
* Allocate temporary memory to store filtered data
*
      IF (IADDR_AA.NE.0) THEN
         CALL FREE_VM(LENGTH_AA,IADDR_AA)
      ENDIF
      LENGTH_AA = NRECORD
      IER = SIC_GETVM (LENGTH_AA,IADDR_AA)
      IP_AA = POINTER(IADDR_AA,MEMORY)
*
*
* For each channel selected, find spikes
*
      I = 1
      DO WHILE (SELECTED_CHANNEL(I).NE.0.AND.I.LE.NCHAN)
         IF (.NOT.RECEIVER_FLAG(SELECTED_CHANNEL(I))) THEN
            ICHAN    = SELECTED_CHANNEL(I)
            IP_DATAR = IP_SIGNAL + (ICHAN-1)*NRECORD
*
*          Restore interpolated values
            IP_NSPIKE      = IPSPIKE+ICHAN-1
            IP_INDEX_SPIKE = IPSPIKE+NCHAN+
     $      2*(ICHAN-1)*LENGTH_SPIKE
            IP_DATA_SPIKE  = IP_INDEX_SPIKE + LENGTH_SPIKE
*DB:03-96    CALL RESTORE_DATA(MEMORY(IP_NSPIKE),MEMORY(IP_INDEX_SPIKE),
*     $      MEMORY(IP_DATA_SPIKE),MEMORY(IP_DATAR))
            NSPIKE = -1
            ITER = 1
            NSPIKE_TOTAL = 0
            DO WHILE (NSPIKE.NE.0 .AND. ITER.LE.ITEMAX)
               NSPIKE = 0
               IF (IFUNC.EQ.1) THEN
                  CALL DESPIKE_SINC (MEMORY(IP_DATAR),MEMORY(IP_AA),
     $            FUNWID)
               ELSEIF (IFUNC.EQ.2) THEN
                  CALL DESPIKE_DIFF2 (MEMORY(IP_DATAR),MEMORY(IP_AA))
               ELSEIF (IFUNC.EQ.3) THEN
                  CALL DESPIKE_LEVEL (MEMORY(IP_DATAR),MEMORY(IP_AA))
               ENDIF
*
* Compare with xrms * mean noise
*
               CALL COMPARE_TO_XRMS (MEMORY(IP_DATAR),MEMORY(IP_AA),
     $         MEMORY(IP_LST),XRMS,IWINRMS,NSPIKE,ICHAN)
*
               IF (NSPIKE.NE.0) WRITE(6,*)'    -> iteration ',ITER,
     $         ': number of spikes found =',NSPIKE
               ITER = ITER+1
               NSPIKE_TOTAL = NSPIKE_TOTAL + NSPIKE
            ENDDO
            IF (NSPIKE_TOTAL.NE.0) THEN
               WRITE(6,'(a,i2,a,i4/)') ' Channel ',ICHAN,
     $         ' ->  Total number of spikes found = ',NSPIKE_TOTAL
               SPIKE_REMOVAL_FLAG(ICHAN) = NSPIKE
               WRITE (DO_BLANKING,'(a,e15.9,2h 1)') 'set blanking ',
     $         BLANKING
               CALL GR_EXEC1 (DO_BLANKING)
*
*            Sort spikes
               CALL SORT_SPIKE(MEMORY(IP_DATA_SPIKE),
     $         MEMORY(IP_INDEX_SPIKE),MEMORY(IP_NSPIKE))
            ENDIF
*
*          Interpolate on spikes
            IF (INTERPOLATE_ON_SPIKE) THEN
               CALL INTERPOLATE_SPIKE(MEMORY(IP_NSPIKE),
     $         MEMORY(IP_INDEX_SPIKE),MEMORY(IP_DATA_SPIKE),
     $         MEMORY(IP_DATAR))
*             Remove spikes of spike area
               CALL RAZ_MEMORY(MEMORY(IP_NSPIKE),1)
               CALL RAZ_MEMORY(MEMORY(IP_INDEX_SPIKE),LENGTH_SPIKE)
               CALL RAZ_MEMORY(MEMORY(IP_DATA_SPIKE), LENGTH_SPIKE)
            ENDIF
         ENDIF
         I = I+1
      ENDDO
*
* Free memory
*
      IF (IADDR_AA.NE.0) THEN
         CALL FREE_VM(LENGTH_AA,IADDR_AA)
         IADDR_AA=0
      ENDIF
      RETURN
      END
*
C --------------------------------------------------------------------
C
C DESPIKE_SINC : find and remove (if required) spikes
C
C Algorythm: CLIP SPIKES FROM THE DATA AT LEVEL OF XRMS TIMES THE NOISE.
C
C Input:
C Output:
C
C --------------------------------------------------------------------
*
      SUBROUTINE DESPIKE_SINC (DATAR,DATANEW,FUNWID)
*
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
      INCLUDE 'const.inc'
*
      REAL*4    DATAR(*),DATANEW(*),FUNWID
*
      REAL*4    CONVD,CONVW
      REAL*4    TOTWID
      REAL*4    X,ZQ
      LOGICAL FLAG
*
      INTEGER*4 I,J,ICONO,N
*
      REAL*4    MEAN
*
* compute mean
*
      MEAN = 0
      N = 0
      DO I=1,NRECORD
         IF (ABS(DATAR(I)-BLANKING).GT.0.5) THEN
            MEAN = MEAN+DATAR(I)
            N=N+1
         ENDIF
      ENDDO
      MEAN = MEAN/N
*
C Compute FUNWID and TOTWID (in arcseconds)
c         FUNWID = VLIGHT / (FREQUENCY*1.D9) / (1.67*telecope_diameter)
c     $            /PI*180.*3600.
c
c         FACTOR = VAZM * 60. * 0.5 ! FACTOR is the number of arcsecond
c                               ! represented by one measured point
c                               ! 0.5 second is the integration time
c
c         write(6,'(a,a,f5.2,a,$)')
c     $        'Convolution with sinc(x/width).',
c     $        ' Width (default =',funwid/factor, '): '
c         read(5,*) funwid
*
      TOTWID = 4.D0 * FUNWID
*
* Compute FUNC() if the function has not yet been initialized
*
      IF (.NOT.INIT_FUNC) THEN
         CALL FUNGEN(FUNC)               ! func(i)=0.01*sin(pi*i/50)/(pi*i/50) (
         INIT_FUNC = .TRUE.
      ENDIF
*
* Find spikes
*
      DO  I = 1,NRECORD
         IF (ABS(DATAR(I)-BLANKING).GT.0.5) THEN
*
            CONVD = DATAR(I) * FUNC(1)
            CONVW = FUNC(1)
*
            IF (I.GT.1) THEN             ! Add to CONVD  datar(j)*sinc() with j<
               FLAG = .TRUE.
               J = I-1
               DO WHILE (J.GE.1 .AND. FLAG)
                  X = ABS(I-J)
                  IF (X.LE.TOTWID) THEN
                     X = X / FUNWID
                     ICONO = INT (50.D0 * X + 0.5D0) + 1
                     ZQ = FUNC(ICONO)
                     IF (ABS(DATAR(J)-BLANKING).GT.0.5) THEN
                        CONVD = CONVD + DATAR(J) * ZQ
                     ELSE
*                       convd = convd+(mean+(drand48()-0.5d0)*2.*rms)*zq
*                       or convd = mean + rangau(rms) (cf util/gag/sysfor.f)
                        CONVD = CONVD + MEAN*ZQ
                     ENDIF
                     CONVW = CONVW + ZQ
                  ELSE
                     FLAG = .FALSE.
                  ENDIF
                  J = J-1
               ENDDO
            ENDIF
*
            IF (I.LT.NRECORD) THEN       ! Add to CONVD  datar(j)*sinc() with j>
               FLAG = .TRUE.
               J = I+1
               DO WHILE (J.LE.NRECORD .AND. FLAG)
                  X = ABS(J-I)
                  IF(X.LE.TOTWID) THEN
                     X = X / FUNWID
                     ICONO = INT(50.D0 * X + 0.5D0) + 1
                     ZQ = FUNC(ICONO)
                     IF (ABS(DATAR(J)-BLANKING).GT.0.5) THEN
                        CONVD = CONVD + DATAR(J) * ZQ
                     ELSE
*                    convd = convd+(mean+(drand48()-0.5d0)*2.*rms)*zq
                        CONVD = CONVD + MEAN*ZQ
                     ENDIF
                     CONVW = CONVW + ZQ
                  ELSE
                     FLAG = .FALSE.
                  ENDIF
                  J = J+1
               ENDDO
            ENDIF
*
            CONVD = CONVD / CONVW
*
            DATANEW(I) = DATAR(I)-CONVD
         ELSE
            DATANEW(I) = 0
         ENDIF
*
      ENDDO
      RETURN
      END
*
C --------------------------------------------------------------------
C
C DESPIKE_DIFF2 :
C
C Algorythm: compute  datar(i)-(datar(i-1)+datar(i+1))/2.
C
C Input:
C
C Output:
C
C --------------------------------------------------------------------
*
      SUBROUTINE DESPIKE_DIFF2 (DATAR,DATANEW)
*
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
*
      REAL*4    DATAR(*),DATANEW(*)
*
      INTEGER*4 I,J
      REAL*4 MEAN
*
* Compute mean
*
      MEAN = 0
      J = 0
      DO I=1,NRECORD
         IF (ABS(DATAR(I)-BLANKING).GT.0.5) THEN
            MEAN = MEAN+DATAR(I)
            J=J+1
         ENDIF
      ENDDO
      MEAN = MEAN/J
*
* Compute DATANEW(I)=DATAR(I)-(DATAR(I-1)+DATAR(I+1))/2.
*
      DATANEW(1) = DATAR(1)-MEAN
      IF (DATAR(1).EQ.BLANKING) DATANEW(1) = 0
      DATANEW(NRECORD) = DATAR(NRECORD)-MEAN
      IF (DATAR(NRECORD).EQ.BLANKING) DATANEW(NRECORD) = 0
*
      DO  I = 2,NRECORD-1
         IF (DATAR(I).GT.BLANKING) THEN
            IF (ABS(DATAR(I-1)-BLANKING).GT.0.5) THEN
               DATANEW(I) = DATAR(I) - DATAR(I-1)/2.
            ELSE
               DATANEW(I) = DATAR(I)-MEAN/2.
            ENDIF
            IF (ABS(DATAR(I+1)-BLANKING).GT.0.5) THEN
               DATANEW(I) = DATANEW(I)-DATAR(I+1)/2.
            ELSE
               DATANEW(I) = DATANEW(I)-MEAN/2.
            ENDIF
         ELSE
            DATANEW(I) = 0
         ENDIF
      ENDDO
*
      RETURN
      END
*
C --------------------------------------------------------------------
C
C DESPIKE_LEVEL :
C
C --------------------------------------------------------------------
*
      SUBROUTINE DESPIKE_LEVEL (DATAR,DATANEW)
*
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
*
      REAL*4    DATAR(*),DATANEW(*)
*
      INTEGER*4 I,J
      REAL*4 MEAN
*
* Compute mean
*
      MEAN = 0
      J = 0
      DO I=1,NRECORD
         IF (ABS(DATAR(I)-BLANKING).GT.0.5) THEN
            MEAN = MEAN+DATAR(I)
            J=J+1
         ENDIF
      ENDDO
      MEAN = MEAN/J
*
* If condition
*
      DO  I = 1,NRECORD
         IF (ABS(DATAR(I)-BLANKING).GT.0.5) THEN
            DATANEW(I) = DATAR(I)-MEAN
         ELSE
            DATANEW(I) = 0
         ENDIF
      ENDDO
*
      RETURN
      END
*
* -----------------------------------------------------------------
*   FUNGEN :
*       this subroutine sets up the convolving function used to         *
*       filter and regrid the scan data.                                *
* -----------------------------------------------------------------
      SUBROUTINE FUNGEN(FUNC)
*
      REAL*8  U,XF,XF1,FUNC(*)
      INTEGER I
      U = 1.0D0
      XF = DACOS(-1.D0)/50.0D0
      XF1 = 0.0D0
      FUNC(1) = 0.01D0
      DO I = 2,201
         XF1 = XF1 + XF
         IF (I.GE.151) THEN
            U = ( 201.D0 - I ) / 50.D0
         ENDIF
         FUNC(I) = 0.01D0 * U * DSIN(XF1) / XF1
      ENDDO
      RETURN
      END
*
* --------------------------------------------------------------------
*
* COMPARE_TO_XRMS
*
* --------------------------------------------------------------------
*
C+HPUX
c$OPTIMIZE LEVEL1
C-HPUX
      SUBROUTINE COMPARE_TO_XRMS (DATAR,DATANEW,LST,XRMS,IWINRMS,
     $NSPIKE,ICHAN)
*
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
      REAL*4    XRMS,DATAR(*),DATANEW(*),LST(*)
      INTEGER*4 IWINRMS,NSPIKE,ICHAN
*
      REAL*4    SD,SD1,AM,ABS_DATA
      INTEGER*4 HWINRMS,I
*
* Compare with xrms * mean noise
*
      AM = 0
      IF (IWINRMS.EQ.0) THEN
         SD = 0.
         DO  I = 1,NRECORD
            AM = MAX (ABS(DATANEW(I)),AM)
            SD = SD + DATANEW(I)*DATANEW(I)
         ENDDO
         SD = SD / NRECORD
         SD = SQRT(SD)
         SD1 = SD*XRMS
         HWINRMS = 0
      ELSE
         SD = 0.
         HWINRMS = IWINRMS/2
         DO I=1,IWINRMS+1
            SD = SD + DATANEW(I)*DATANEW(I)
         ENDDO
         DO I=1,HWINRMS+2,NRECORD-HWINRMS-1
            AM = MAX (ABS(DATANEW(I)),AM)
         ENDDO
      ENDIF
      AM = AM*0.999
*
*
* Find spikes
      IF (IWINRMS.NE.0) THEN             ! compute rms on window
         DO I=1,NRECORD
            IF (I.GT.(HWINRMS+1) .AND. I.LE.(NRECORD-HWINRMS)) THEN
               SD = SD - DATANEW(I-(HWINRMS+1))**2
     $         + DATANEW(I+HWINRMS)**2
            ENDIF
            SD1 = SD/(IWINRMS+1)
            SD1 = SQRT(SD1)*XRMS
*
            ABS_DATA = ABS(DATANEW(I))
            IF (ABS_DATA.GT.SD1 .AND. ABS_DATA.GE.AM) THEN
               WRITE(6,'(a,i5,a,f8.2,a,f8.0)')
     $         '     Record=',I,' LST=',LST(I),'  SIGNAL=',DATAR(I)
               NSPIKE = NSPIKE+1
               CALL ADD_SPIKE(I,DATAR(I),ICHAN)
               DATAR(I) = BLANKING
            ENDIF
         ENDDO
      ELSE
         DO I=1,NRECORD
            IF (ABS(DATANEW(I)).GT.SD1.AND.ABS(DATANEW(I)).GE.AM) THEN
               WRITE(6,'(a,i5,a,f8.2,a,f8.0)')
     $         '     Record=',I,' LST=',LST(I),'  SIGNAL=',DATAR(I)
               NSPIKE = NSPIKE+1
               CALL ADD_SPIKE(I,DATAR(I),ICHAN)
               DATAR(I) = BLANKING
            ENDIF
         ENDDO
      ENDIF
      RETURN
      END
*
*------------------------------------------------------------------------
*
* RESTORE_DATA : restore original values instead of interpolated values
*
*------------------------------------------------------------------------
      SUBROUTINE RESTORE_DATA (NSPIKE_CHAN,SPIKE_RECORD,SPIKE_DATA,
     $DATAR)
*
      INTEGER NSPIKE_CHAN                ! number of spikes for this channel
      INTEGER SPIKE_RECORD(*)            ! number of record associated
      REAL    SPIKE_DATA(*)              ! data associated to this record
      REAL    DATAR(*)
      INTEGER I
*
      DO I =1,NSPIKE_CHAN
         DATAR(SPIKE_RECORD(I)) = SPIKE_DATA(I)
      ENDDO
      NSPIKE_CHAN = 0                    ! raz number of spikes
      END
*
*------------------------------------------------------------------------
*
* ADD_SPIKE : add a spike in spike arrays
*
* Spike Memory Organization (Total = nchan+2*nchan*length_spike)
*
* ipspike 				= number of spikes for each receiver
*     +nrec [1..length_spike]		= record number of spikes (receiver 1)
*     +length_spike [1..length_spike]	= signal of spikes (receiver 1)
*
* ... and so on
*
* nrec = number of receivers
* length_spike allocated in steps of 100 spikes
*
*------------------------------------------------------------------------
*
      SUBROUTINE ADD_SPIKE (IRECORD,VAL,ICHAN)
*
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
*
      INTEGER IRECORD                    ! index of record (input)
      REAL    VAL                        ! = signal[irecord,ichan] (input)
*
      INTEGER*4 IP,IADDR,POINTER,IP1,IP2
      INTEGER SIC_GETVM
      INTEGER ICHAN
      INTEGER NSPIKE                     ! number of spikes for this channel
      INTEGER OLD_LENGTH_SPIKE,LENGTH,IER,I
*
* -- Test if enough memory to add a spike
      CALL R4TOR4(MEMORY(IPSPIKE+ICHAN-1),NSPIKE,1)  !number of spikes before
      NSPIKE = NSPIKE + 1                ! new number of spikes
      IF (NSPIKE.GT.LENGTH_SPIKE) THEN
*
*       Allocate more memory for spikes data
         OLD_LENGTH_SPIKE = LENGTH_SPIKE
         LENGTH_SPIKE = LENGTH_SPIKE + 100
         LENGTH = NCHAN + 2*NCHAN*LENGTH_SPIKE
         IER = SIC_GETVM(LENGTH,IADDR)
         IP =  POINTER(IADDR,MEMORY)
*
*      Copy information
         CALL R4TOR4(MEMORY(IPSPIKE),MEMORY(IP),NCHAN)
         DO I=1,NCHAN
            IP1 = IPSPIKE + NCHAN+ 2*(I-1)*OLD_LENGTH_SPIKE
            IP2 = IP      + NCHAN+ 2*(I-1)*LENGTH_SPIKE
            CALL R4TOR4(MEMORY(IP1),MEMORY(IP2),MEMORY(IPSPIKE+I-1))
            IP1 = IP1 + OLD_LENGTH_SPIKE
            IP2 = IP2 + LENGTH_SPIKE
            CALL R4TOR4(MEMORY(IP1),MEMORY(IP2),MEMORY(IPSPIKE+I-1))
         ENDDO
*
*       Free "old" memory
         LENGTH = NCHAN + 2*NCHAN*OLD_LENGTH_SPIKE
         CALL FREE_VM(LENGTH,IADDRSPIKE)
*
*       Update IPSPIKE and IADDRSPIKE
         IPSPIKE = IP
         IADDRSPIKE = IADDR
      ENDIF
*
*
* -- Add a spike in spike array
*
*
      CALL R4TOR4(NSPIKE,MEMORY(IPSPIKE+ICHAN-1),1)  ! number of spikes
      IP = IPSPIKE + NCHAN+ 2*(ICHAN-1)*LENGTH_SPIKE + NSPIKE - 1
      CALL R4TOR4(IRECORD,MEMORY(IP),1)  ! index
      IP = IP+LENGTH_SPIKE
      CALL R4TOR4(VAL,MEMORY(IP),1)      ! data
      RETURN
      END
*
*------------------------------------------------------------------------
*
* SORT_SPIKE : sort spike area
*
*------------------------------------------------------------------------
      SUBROUTINE SORT_SPIKE (DA,RA,N)
      INTEGER N, RA(N), RRA, I, J, L, IR
      REAL DA(N), DDA
*
      IF (N.EQ.1) RETURN
      L = N/2+1
      IR = N
      DO WHILE (N.NE.(N+1))
         IF (L.GT.1) THEN
            L = L-1
            RRA = RA(L)
            DDA = DA(L)
         ELSE
            RRA = RA(IR)
            DDA = DA(IR)
            RA(IR) = RA(1)
            DA(IR) = DA(1)
            IR = IR-1
            IF (IR.EQ.1) THEN
               RA(1) = RRA
               DA(1) = DDA
               RETURN
            ENDIF
         ENDIF
         I = L
         J = L+L
         DO WHILE (J.LE.IR)
            IF (J.LT.IR) THEN
               IF (RA(J).LT.RA(J+1)) J=J+1
            ENDIF
            IF (RRA.LT.RA(J)) THEN
               RA(I) = RA(J)
               DA(I) = DA(J)
               I = J
               J = J+J
            ELSE
               J = IR+1
            ENDIF
         ENDDO
         RA(I) = RRA
         DA(I) = DDA
      ENDDO
      END
*
*
*------------------------------------------------------------------------
*
* INTERPOLATE_SPIKE : interpolate when spike
*
*------------------------------------------------------------------------
      SUBROUTINE INTERPOLATE_SPIKE (NSPIKE_CHAN,SPIKE_RECORD,SPIKE_DATA,
     $DATAR)
*
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
*
      INTEGER NSPIKE_CHAN                ! number of spikes for this channel
      INTEGER SPIKE_RECORD(*)            ! number of record associated
      REAL    SPIKE_DATA(*)              ! data associated to this record
      REAL    DATAR(*)
*
      INTEGER INDEX,ISPIKE,J,DIFF,LSPIKE
      REAL    DATA1,DATA2
*
      ISPIKE = 1
      DO WHILE (ISPIKE.LE.NSPIKE_CHAN)
*
* ---- Find "length of spike"
         LSPIKE = 1
         DIFF = 1
         J = ISPIKE
         DO WHILE (DIFF.EQ.1 .AND. (J+1).LE.NSPIKE_CHAN)
            DIFF = SPIKE_RECORD(J+1) - SPIKE_RECORD(J)
            IF (DIFF.EQ.1) THEN
               LSPIKE = LSPIKE+1
               J      = J +1
            ENDIF
         ENDDO
*
* ---- Interpolate
         INDEX = SPIKE_RECORD(ISPIKE)
         IF (INDEX.GT.1) THEN
            DATA1 = DATAR(INDEX-1)
         ELSE
            DATA1 = DATAR(INDEX+LSPIKE)
         ENDIF
         IF (INDEX.GT.1) THEN
            DATA2 = DATAR(INDEX+LSPIKE)
         ELSE
            DATA2 = DATAR(INDEX-1)
         ENDIF
         IF (ABS(DATA1-BLANKING).LT.1. .AND.
     $   ABS(DATA2-BLANKING).GT.05) DATA1 = DATA2    ! data1=blanking
         IF (ABS(DATA2-BLANKING).LT.0.5 .AND.
     $   ABS(DATA1-BLANKING).GT.1.) DATA2 = DATA1    ! data2=blanking
         DO J=1,LSPIKE
            DATAR(INDEX+J-1) = (DATA1+DATA2)/2.
         ENDDO
         ISPIKE = ISPIKE + LSPIKE
      ENDDO
      RETURN
      END
