! Procedure to correct for platforming effects in FTS spectra
! at 200kHz and 50kHz resolution
! written by CK and MG based on ideas of JeP
! 2-Feb-2012
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Delete Global Variables if they exist!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if exist(channelWindow) then
del /var channelWindow
endif

if exist(chan) then
del /var chan
endif

if exist(velo) then
del /var velo
endif

if exist(debug) then
del /var debug
endif


! Debug option: 1 shows plots during correction and pauses
!               0 does not pause                
def int debug /global
let debug 1
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! List of Lines and thei width in GHz 
!!! which are expected to appear in the 
!!! bandwidth of the Spectra that has to be 
!!! corrected.
!!! Different Receivers can be treated at once
!!! Lines that are not in the bandwidth of the 
!!! actual spectrum are ignored. 
!!! You can add as many lines as you want
!!! Or leave it entirly empty 
!!! the latter defaults to 'set window 0 0'.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

begin data gag_scratch:lines.dat

  115271	50
  110201	50
  230542	50
  220402	50
  
end data gag_scratch:lines.dat


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Procedure to get the channel for a specific Frequency
!
! Input &1 Frequency in GHz
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

begin procedure freqToChannel

! read the Input
define real freqInput 
let freqInput &1

define real diffFreqCentralFreq
define int channelDiff

! variable 'chan' stores the corresponding channel
! it has to be global for later use in  prodecure win

if exist(chan) then
del /var chan
endif

define int chan /global

! calculate the difference between central frequency 
! (stored in R%HEAD%SPE%RESTF) and the input frequency 
let diffFreqCentralFreq freqInput-R%HEAD%SPE%RESTF

! get the number of channels corresponding to 
! that difference
let channelDiff diffFreqCentralFreq/freq_step

! add the difference in channels to the 
! central channel 
let chan R%HEAD%SPE%RCHAN+channelDiff
exa chan

end procedure freqToChannel

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
! Get the Velocity for a specific Frequency
!
! Input &1 Frequency in GHz
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


begin procedure freqToVelocity

! read the Input
define real freqInput 
let freqInput &1

define real diffFreqCentralFreq
define int veloDiff

! variable 'velo' stores the corresponding velocity
! it has to be global for later use in  prodecure win

if exist(velo) then
del /var velo
endif

define int velo /global

! calculate the difference between central frequency 
! (stored in R%HEAD%SPE%RESTF) and the input frequency 
let diffFreqCentralFreq freqInput-R%HEAD%SPE%RESTF

! get the difference in velocity 
! to do so; get the number of Channels corresponding 
! to diffFreqCentralFreq and multiply by the 
! channel width in velocity
let veloDiff velo_step*(diffFreqCentralFreq/freq_step)
 
! add the velocity difference to the central velocity
let velo R%HEAD%SPE%VOFF+veloDiff
exa velo

end procedure freqToVelocity

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
! Construct the line windows from 
!
! Input &1 windowSize in MHz
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
begin procedure win
  
  ! read in the Line Frequencies 
  ! for that windows with width 'windowSize'
  ! are to be created
 
  greg1\column x 1 y 2 /file gag_scratch:lines

  ! lenght of the array storing the window 
  ! number of line (nxy) times 2 for the two window borders
  ! we have to add +1 since "set wind /var " expects 
  ! the first entry in the array to indicate the number of 
  ! windows
  def int lenght
  let lenght nxy*2+1
  if (lenght.gt.1) then
 	 ! Windows is the array that stores all line Windows 
 	 if exist(Windows) then
 		  del /var Windows
  	endif
	def real Windows[lenght] /global
	
	! first entry has to be the window number in the array
	let Windows[1] nxy
	exa Windows

	! Now fill the array with the corresponding window 
	!borders in velocity
	for i 1 to nxy by 1
		say 'x[i]-y[i]'
		@freqToVelocity 'x[i]-y[i]'
		let Windows['2*i'] velo
		say 'x[i]+y[i]'
		@freqToVelocity 'x[i]+y[i]'
		let Windows['2*i+1'] velo
	next i

	set wind /var Windows

  else if (lenght.eq.1) then
	set wind 0 0
  endif

end procedure win


begin procedure subband 

  set variable spectro read
  define real jump1 jump2 nfts dum
  define real ty ty2 ty3 /like ry
  define integer band reso
  
  !define int reso 
  let reso 200
  say "Resolution="'reso'"kHz

  if ((freq_step.ge.0.194).and.(freq_step.le.(0.200))) then  !reso 1 = 200kHz
    say "Here"
    exa nchan ! 20727 = number of channels from header
    exa rchan !  8961 = reference channel
    exa fres  ! 0.1953
    exa vres
    !if ((fres.gt.0.2).or.(fres.lt.0.19)) then
    !  pause " ERROR: fres.ne.0.195 MHz "
    !endif
    let nfts = (nchan+2)|3 
    let jump1 = nfts   ! edge in number of channels
    let jump2 = 2*nfts 
    let jump1 = vres*(jump1-rchan)+voff ! edges in velocity
    let jump2 = vres*(jump2-rchan)+voff

  else if ((freq_step.ge.0.045).and.(freq_step.le.(0.055))) then

    exa nchan 
    exa rchan 
    exa fres  
    exa vres
    if ((fres.gt.0.06).or.(fres.lt.0.04)) then
      pause " ERROR: fres.ne.50kHz "
    endif
    let jump1 = 12494 
    let jump2 = 12494+12289 ! information from GP, 2-Feb-2012
    let jump1 = vres*(jump1-rchan)+voff ! edges in velocity
    let jump2 = vres*(jump2-rchan)+voff

  else 
    pause " ERROR: wrong resolution "
  end if

  if (jump1.lt.jump2) then
    continue
  else 
    let dum = jump1
    let jump1 = jump2
    let jump2 = dum
  endif
  say " Band edges = "'jump1'" "'jump2'" km/s " ! -1227 517 km/s

  ! this gets the line windows for the current spectrum
  @win 
  
  if (debug.eq.1) then
 	 plot
	 pause "original spectrum"  
  endif
  ! subtract baselines from each of the three FTS units in the 4GHz band:
 
  let ry -1000.00 /where rx.gt.jump1
  if (debug.eq.1) then
  	pl
  	pause
  endif
  if (debug.ne.1) then
	  base 1 
  else if (debug.eq.1) then
	base 1 /plot; draw window; 
	pause "1st FTS unit handled "
  endif
  let ty = ry /where rx.lt.jump1
  

  get
  let ry -1000.00 /where (rx.lt.jump1)
  let ry -1000.00 /where (rx.gt.jump2)

  if (debug.eq.1) then
  	pl
  	pause
  endif 
  if (debug.ne.1) then
	  base 1 
  else if (debug.eq.1) then
 	 base 1 /plot; draw window; 
  	pause
  endif 
  let ty2 = ry /where ((rx.gt.jump1).and.(rx.lt.jump2))
 
  get 
  let ry -1000.00 /where rx.lt.jump2

  if (debug.eq.1) then
  	pl
  	pause
  endif 

  if (debug.ne.1) then
	  base 1 
  else if (debug.eq.1) then
 	 base 1 /plot; draw window; 
  	pause "3rd FTS unit handled"
  endif 
  let ty3 = ry /where rx.gt.jump2
   

  let ry = ty /where rx.lt.jump1
  let ry = ty2 /where (rx.gt.jump1).and.(rx.lt.jump2)
  let ry = ty3 /where rx.gt.jump2


  if (debug.eq.1) then
  	plot
        draw wind
	pause " new spectrum with FTS edges "
  endif 
  !smooth box 10; plot; 
  !set window jump1 jump2
  !draw window
  

  delete /var jump1 jump2 nfts dum ty ty2 ty3 band reso

end procedure subband
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


pen 0 /weight 3
define real test
set plot histo
set format long

file in FTSOdp20120215.30m
file out fts-red.30m single /overwrite

set tel *E0*UO*
set line *12CO*
set source NGC604
find


for i 1 to found
 get n
 @ subband 
 write
next