! Procedure of weighted mean : wmean.nic
! Can be used for the signal subscans signal[ &3 : (&3+npts-1) , nchan]
!			 a ttempo[nrecord,1]
!			 a flux[npairs,nchan]
! &1 is the input variable to mean, its dimension is [&4,nch]
! &2 is the input variable for the weights, its dimension is [&4,nch]
! &3 is the index of the first data to consider on the variable &1
! &4 is the number of data to consider on the variable &1
! &5 is the estimated mean, its dimension is [nch]
! &6 is the standard deviation of this mean,its dimension is [nch]
! &7 is the second dimension of the variable &1 (can be 1 or npairs or nchan)
!
! Extraction of data to mean
define integer i1 npts ix nch
let nch &7
let i1 &3
let npts &4
define real x[npts,nch] xwei[npts,nch]
for k 1 to nch
   for i 1 to npts
      let ix i1-1+i
      let x[i,k] &1[ix,k]
      let xwei[i,k] &2[ix,k]
   next
next
del /var i1 ix
!
! Calculate the mean
if npts.le.1 then
   say "No meaning to do"
   let &5 &1[&3]
   let &6 0
else
   define real xmean[nch] xrms[nch] weisum[nch]
   for k 1 to nch
      for i 1 to npts
         let xmean[k] xmean[k]+xwei[i,k]*x[i,k]
         let weisum[k] weisum[k]+xwei[i,k]
      next
      if weisum[k].ne.0 then
         let xmean[k] xmean[k]|weisum[k]
      else
         say "All the weights are zero !"
         let xmean[k] 0
      endif
      for i 1 to npts
         let xrms[k] xrms[k]+xwei[i,k]*(x[i,k]-xmean[k])**2
      next
      if weisum[k].ne.0 then
         let xrms[k] sqrt(xrms[k]|((npts-1)*weisum[k]))
      else
         let xrms[k] 0
      endif
   next
   let &5 xmean
   let &6 xrms
   del /var weisum xmean xrms
endif
del /var x xwei npts nch