!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 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 @spikes3. 
! 3.- This wonderful program will give you the
!     Local oscillator and IF frequencies of all
!     the spikes in the loaded spectrum, and a plot
!     with the location of the spikes.
!
! We start by calculating the rms of the spectrum
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

! Now we define all the variables used in the macro

define double fspike            !!Frequency of the spike (MHz)
define double fcentral          !!Central sky frequency (MHz)
define double fimage            !!Image frequency (MHz)
define double flo               !!Frequency of the local oscillator (MHz)
define double fif               !!IF frequency (+-6250 or +-9430) (MHz)
define double fact              !!Doppler correction
define double ifreq             !!IF frequency of the spike (GHz)
define integer Nchannel         !!Number of channels of the spectrum
define integer j                !!Counters
define integer k
define integer count            !!Counter on the number of spikes
define double rms               !!rms (K)
define double pos_spike[100]    !!Vector with the position of the spikes (MHz)
define double if_spike[100]     !!Vector with the IF of the spikes (GHz)
define double freq_spike[100]   !!Vector with the frequencies of the spikes(GHz)
define real yrange ypos

!!We assignate the values to the things which are common to all the spectrum,
!!i.e. the Number of channels, the rms, ....

let Nchannel R%HEAD%SPE%NCHAN
let rms R%HEAD%BAS%SIGFI
let count 1

!!We indicate the limits of the window for plot purposes...
let yrange user_ymax-user_ymin
if (pro%narg.eq.1) then
   let ypos 'pro%arg[1]'
else
   let ypos user_ymin+yrange*.7
endif

!Characterization of the system (frequencies of the local oscilator
! image, doppler correction... All the frequencies are in MHz !!!)

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

!We make a loop over all the channels of the spectrum. If the intensity
!of the channel is high (i.e. > 5rms) and if the intensity of the continous
!channels is small (<4 times the intensity of the channel) we would identify
!it as a spike. The maximum accepted width to identify it as a spike
!is 2 channels, so take care when using this macro with narrow spectral lines.

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 pos_spike[count] rx[i]
     let freq_spike[count] fspike
     let if_spike[count] ifreq
     let count count+1

  endif
next

if (count.eq.1)
  say "Congratulations ! No spikes on the spectrum"
else
  plot
  say "------------------------"
  say "f_LO (GHz)" 'flo/1000.0' /format a10 f9.3
  say "------------------------"
  say "Spike Frequency (GHz)" "IF_spike (GHz)"
  say "---------------------" "--------------"

  for i 1 to count-1
     pen 1
     g\draw relo 'pos_spike[i]' ypos  /user
     g\draw line 'pos_spike[i]' ypos+yrange/10 /user
     g\draw text 'pos_spike[i]' ypos+2*yrange/10 'nint(if_spike[i]*100.0)|100.0' 6 90 /user
     pen 0
     say 'freq_spike[i]/1000' "        " 'if_spike[i]' /format f9.3 a8 f9.2
  next
say "---------------------" "--------------"
endif