!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! THE WONDERFUL CLASS PROGRAM TO CATCH SPIKES ! 
! (c) Manuel Gonzalez and Christof Buchbender !
!                 June 2012                   !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!
! HELP
! 1.- Load in CLASS your spiky spectrum
! 2.- Tape @detection_spikes. 
! 3.- This wonderful program will give you the
!     Local oscillator and IF frequencies of all
!     the spikes in the loaded spectrum
!
! We start by calculating the rms of the spectrum
find

define integer j                !!Counters
define integer k
define integer count            !!Counter on the number of spikes
define integer number_spec[100]
define double spike_chn[100]
define double spike_freq[100]
define double spike_if[100]
let count 1

for m1 1 to found

  get n
  set format long
  set plot histo
  set unit f v
  set win 0 10
  sic mess class A-I
  base
  plot
  sic mess class A+I

  define integer Nchannel         !!Number of channels of the spectrum
  define double rms               !!rms (K)
  define double fcentral
  define double fimage
  define double fact
  define double fif
  define double flo
  define double ifreq
  define double fspike
  let Nchannel R%HEAD%SPE%NCHAN
  let rms R%HEAD%BAS%SIGFI
  let fcentral R%HEAD%SPE%RESTF               !!Sky central frequency
  let fimage R%HEAD%SPE%IMAGE                 !!Image frequency
  let fact 1.0/(1.0+R%HEAD%SPE%DOPPLER)       !!Doppler correction
  let fif (fimage-fcentral)/2.0*fact          !!if value (+-6.25, +-9.43)
  let flo fcentral+fif

  for i 2 to Nchannel-2
    let j i-1
    let k i+2 
    if (abs(ry[i]/rms).gt.5.and.abs(ry[i]/ry[k]).gt.4.and.abs(ry[i]/ry[j]).gt.4)
       
       let fspike fcentral+rx[i]           !!Spike frequency (MHz)
       let ifreq abs(fspike-flo)/1000.0    !!IF of the spike GHz()
       if (ifreq.gt.8) then                !!If the IF > 8GHz we fold it into
          let ifreq 15.68-ifreq            !!the good range.
       endif
       let number_spec[count] R%HEAD%GEN%NUM
       let spike_freq[count] R%HEAD%SPE%RESTF+rx[i]
       let spike_chn[count] rx[i]
       let spike_if[count] ifreq
       let count count+1

    endif
  next i
  sic\delete /variable Nchannel rms
  sic\delete /variable fcentral fimage
  sic\delete /variable fact fif
  sic\delete /variable flo ifreq
  sic\delete /variable fspike
next m1

if (count.ne.1) then
  say "============================================"
  say "           LIST OF SPIKY SPECTRA            "
  say "============================================"
  for i 1 to count-1
     say 'i' 'number_spec[i]' 'Spike_chn[i]' 'Spike_freq[i]' 'spike_if[i]' /format  i5  i5 f12.2 f12.2 f9.2
  next
  say "============================================"
  say "                    END                     "
  say "============================================"

endif