! Procedure flux.nic
! Calcule le flux associe au scan considere
! &1 nombre de bolometres a considerer
! &2 flux du ON-OFF
! &3 estimation de l'incertitude de ce flux
! signal
despike/iter 3/interp
signal
!pause
!
! Applique les facteurs de calibration
set gain_amplifier
!set gain_receiver on 1. 1.16 1.13 1. 0.91 0.95 0.96
set gain_chan
set gain_elevation
set opacity
!symbol cal/inquire "Facteur de calibration : "
!var data on write
!let signal signal*cal /where signal.ne.blanking[1]
var *
!
! Determine l'index de depart de chaque 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 le nb de bolometres et le nb de paires 
define integer nbolo npairs/global
let nbolo &1
let npairs NSUBSCAN|2
!
! Calcule la moyenne par subscans
define real ysubs[nbolo,NSUBSCAN] ysubs_sig[nbolo,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] nbolo
next
!say "La moyenne du signal sur chacun des subscans"
!for k 1 to nbolo
!say 'k' 'ysubs[k,1]' 'ysubs_sig[k,1]'
!next
!pause
!
! Soustraction d'une ligne de base d'ordre 1
var data on write
define real x0[npairs] y0[npairs] w0[npairs] A B/global
define integer nph1 nph2
for k 1 to nbolo
   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
!   say "bolometre " 'k' "a*x+b=" '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
!	 if signal[i,k].ne.blanking[1] then 
            let signal[i,k] signal[i,k]-(A*nph+B)
!	 endif
      next
   next
next
!signal
del/var x0 y0 w0 A B nph1 nph2
!var data on read
!
! Calcul de la moyenne des phases positives et negatives
define real mean_ON[nbolo] mean_OFF[nbolo]
define integer nph_ON[npairs] nph1 nph_OFF[npairs] nph2
define real wymON[nbolo] wymOFF[nbolo]
for k 1 to nbolo
   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]
next
!exa mean_ON
!exa mean_OFF
del/var nph1 nph2 wymON wymOFF
!
! Soustraction du bruit de ciel
define real snoisem[NRECORD] snoisew[NRECORD] moy
for k 1 to nbolo
   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
!	 if signal[i,k].ne.blanking[1] then
            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
!         endif
      next
    next
next
let snoisem snoisem|snoisew
!for i 1 to 10
!   say 'snoisem[i]'
!next
del/var moy snoisew mean_ON mean_OFF
!
! Calcul des moyennes par subscans corrigees par le bruit de ciel
let sigwei 1
!var data write
for k 1 to nbolo
   for i 1 to NRECORD
!      if signal[i,k].ne.blanking[1] then
         let signal[i,k] signal[i,k]-snoisem[i]
!      endif
   next
next
var data on read
signal
!pause
for nph 1 to NSUBSCAN
   let n nph
   @wmean signal sigwei index[n] nsubrec[n] ysubs[n] ysubs_sig[n] nbolo
next
!say "La moyenne du signal sur le premier subscan et sur les bolometres"
!for k 1 to nbolo
!say 'k' 'ysubs[k,1]' 'ysubs_sig[k,1]'
!next
!
! Calcul du flux
define real yy[npairs,nbolo] ww[npairs,nbolo] fluxx[nbolo] fluxsig[nbolo]/global
define integer nph1 nph2
for np 1 to npairs
   let nph1 nph_ON[np]
   let nph2 nph_OFF[np]
!   if (mod(np,2).eq.1) then
!      let nph1 2*np-1
!      let nph2 2*np
!   else
!      let nph1 2*np
!      let nph2 2*np-1
!   endif
   for k 1 to nbolo
      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 nbolo
!for k 1 to nbolo
!  say "Flux : " 'flux[k]' "+\-" 'fluxsig[k]'
!next
let &2 fluxx
let &3 fluxsig
del/var yy ww nph1 nph2 fluxx fluxsig nph_ON nph_OFF
del/var snoisem ysubs ysubs_sig sigwei n index npairs nbolo