! Procedure flux.nic
! Calculate the flux associated with the on-off scan
! &1 flux du ON-OFF
! &2 the estimated rms on the flux (measures dispersion)
!
! Determine the begining index of each subscan
define integer index[NSUBSCAN] /global
define integer nph1
let index[1] 1
for nph 2 to NSUBSCAN
   let nph1 nph-1
   let index[nph] index[nph1]+nsubrec[nph1]
next
del /var nph1
!
! Determine the number of subscans pairs
define integer npairs /global
let npairs NSUBSCAN|2
!
! Calculate the mean of the signal subscan by subscan
define real ysubs[nchan,NSUBSCAN] ysubs_sig[nchan,NSUBSCAN] /global
define real sigwei /like signal /global
define integer n /global
let sigwei 1
for nph 1 to NSUBSCAN
   let n nph
   @ wmean signal sigwei index[n] nsubrec[n] ysubs[n] ysubs_sig[n] nchan
next
!
! Calculate and remove a first order baseline on signal and y_subs
var data on write
define real x0[npairs] y0[npairs] w0[npairs] A B /global
define integer nph1 nph2
for k 1 to nchan
  if (.not.chan_flag[k]) then
   for np 1 to npairs
      let nph1 2*np-1
      let nph2 2*np
      let x0[np] 2*np-0.5
      let y0[np] (ysubs[k,nph1]+ysubs[k,nph2])|2
      let w0[np] 1|(ysubs_sig[k,nph1]**2+ysubs_sig[k,nph2]**2)
   next
   @ fitline x0 y0 w0 npairs A B
   for nph 1 to NSUBSCAN
      let ysubs[k,nph] ysubs[k,nph]-(A*nph+B)
      for i index[nph] to index[nph]+nsubrec[nph]-1
         let signal[i,k] signal[i,k]-(A*nph+B)
      next
   next
  endif
next
del /var x0 y0 w0 A B nph1 nph2
!
define logical sky_red
let sky_red .TRUE.
If sky_red.and.(.NOT.noise_flag) then
say "!   Reduction of high frequency skynoise"
! Calculate the mean of the positive and negative phases
define real mean_ON[nchan] mean_OFF[nchan]
define integer nph_ON[npairs] nph1 nph_OFF[npairs] nph2
define real wymON[nchan] wymOFF[nchan]
for k 1 to nchan
  if (.not.chan_flag[k]) then
   for np 1 to npairs
      if (mod(np,2).eq.1) then
         let nph_ON[np]  2*np-1
         let nph_OFF[np] 2*np
      else
         let nph_ON[np]  2*np
         let nph_OFF[np] 2*np-1
      endif
      let nph1 nph_ON[np]
      let nph2 nph_OFF[np]
      if ysubs_sig[k,nph1].ne.0 then
         let mean_ON[k] mean_ON[k]+ysubs[k,nph1]|ysubs_sig[k,nph1]**2
         let wymON[k]   wymON[k]+1|ysubs_sig[k,nph1]**2
         !let mean_ON[k] mean_ON[k]+ysubs[k,nph1]*nsubrec[nph1]
         !let wymON[k] wymON[k]+nsubrec[nph1]
      endif
      if ysubs_sig[k,nph2].ne.0 then
         let mean_OFF[k] mean_OFF[k]+ysubs[k,nph2]|ysubs_sig[k,nph2]**2
         let wymOFF[k]   wymOFF[k]+1|ysubs_sig[k,nph2]**2
         !let mean_OFF[k] mean_OFF[k]+ysubs[k,nph2]*nsubrec[nph2]
         !let wymOFF[k] wymOFF[k]+nsubrec[nph_OFF]
      endif
   next
   let mean_ON[k]  mean_ON[k]|wymON[k]
   let mean_OFF[k] mean_OFF[k]|wymOFF[k]
  endif
next
del /var nph1 nph2 wymON wymOFF
!
! High frequency skynoise reduction
define real snoisem[NRECORD] snoisew[NRECORD] moy
for k 1 to nchan
  if (.not.chan_flag[k]) then
   for nph 1 to NSUBSCAN
      if (mod(nph,4).le.1) then
         let moy mean_ON[k]
      else
         let moy mean_OFF[k]
      endif
      for i index[nph] to index[nph]+NSUBREC[nph]-1
         let snoisem[i] snoisem[i]+(signal[i,k]-moy)|ysubs_sig[k,nph]**2 
         let snoisew[i] snoisew[i]+1|ysubs_sig[k,nph]**2
      next
    next
  endif
next
let snoisem snoisem|snoisew
del /var moy snoisew mean_ON mean_OFF
!
! Calculate the mean of the signal after skynoise correction subscan by subscan
let sigwei 1
for k 1 to nchan
  if (.not.chan_flag[k]) then
   for i 1 to NRECORD
      let signal[i,k] signal[i,k]-snoisem[i]
   next
  endif
next
var data on read
var flag wr
let noise_flag .true.
var flag 
@ p_signal
del /var snoisem
say "!   Now starts the flux computation ..."
for nph 1 to NSUBSCAN
   let n nph
   @ wmean signal sigwei index[n] nsubrec[n] ysubs[n] ysubs_sig[n] nchan
next
!
else
   define integer nph_ON[npairs] nph_OFF[npairs]
   for np 1 to npairs
      if (mod(np,2).eq.1) then
         let nph_ON[np]  2*np-1
         let nph_OFF[np] 2*np
      else
         let nph_ON[np]  2*np
         let nph_OFF[np] 2*np-1
      endif
   next
endif
!
!
! Calculate the flux
define real yy[npairs,nchan] ww[npairs,nchan] fluxx[nchan] fluxsig[nchan] /global
define integer nph1 nph2
for np 1 to npairs
   let nph1 nph_ON[np]
   let nph2 nph_OFF[np]
   for k 1 to nchan
      let yy[np,k] (ysubs[k,nph1]-ysubs[k,nph2])|2
      let ww[np,k] 1|(ysubs_sig[k,nph1]**2+ysubs_sig[k,nph2]**2)
   next
next
@ wmean yy ww 1 npairs fluxx fluxsig nchan
let &1 fluxx
let &2 fluxsig
del /var yy ww nph1 nph2 fluxx fluxsig nph_ON nph_OFF
del /var ysubs ysubs_sig sigwei n index npairs