!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!	this subroutine compute the mean and the rms 
!	of mesurements with different rms.
!
!	R.MORENO 18.02.95
!
define integer ntab ntab2 bid1 iscan
define integer nchannel
define character dummy*30 filename*40
!
let dummy 'scan_list'"None"
if (dummy.ne."None") then
 sic output "scan.tmp"
 say 'scan_list' /format a255
 sic output
 system "cp scan.tmp scan_list.tmp"
 system "rm scan.tmp"
!
 let nmax 0
 def real t1[100,100] t2[100,100]
 def real ti[100]
 def int ii
 def char c*1
 for j 1 to 100
   col x j /file "scan_list.tmp"
   if (nxy.ne.0) then
     let nmax nmax+1
     let iscan nint(x)
     if (iscan.lt.10) then
       let filename "000"'iscan'".onf"
     else if (iscan.lt.100) then
       let filename "00"'iscan'".onf"
     else if (iscan.lt.1000) then
       let filename "0"'iscan'".onf"
     else
       let filename 'iscan'".onf"
     endif
     sic search 'filename'
     if (.not.sic$exist) then
       say "Cannot open "'filename'
       let onoff_error yes    
       return
     endif
     col x 1 y 6 z 7 /file 'filename'
     if (nxy.eq.0) then
       say 'filename'" is empty"
       let onoff_error yes    
       return
     endif
     if (j.eq.1) then    ! find number of channels
       let nchannel = x[1]
       for i 2 to nxy
	 if (x[i].gt.nchannel) then
		let nchannel = x[i]
	 end if
       next i
       if (nchannel.le.0) then
	 say "Incorrect file "'filename'
         let onoff_error yes
         return
       endif
       if (nchannel.gt.100) then
         say "Number of channels found greater than 100"
       endif
     endif
! read nrecord in each scan to compute integration time
     get 'iscan'
     if (olam.ne.0.or.obet.ne.0) then
       say "NIC is not yet ready to combine scans with RA or DEC offsets"
       let  onoff_error yes
       return
     endif
     let ti[j] nrecord
     for i 1 to nchannel
	let t1[i,j] = y[i]
	let t2[i,j] = z[i] 
     next i
! Unit (jy or counts)
     system "touch .tmp"
     col x 1 /file ".tmp"         ! just to close 'filename'
     accept c /col 'filename' /line 6 /format "58x,1a"
     let ii j
     if (c.ne."J".and.c.ne."s") then
      say "Cannot retrieve unit in "'filename'
      say "Unit will be fixed according to jansky_flag"
      let onoff_unit jansky_flag
     else 
      if (ii.eq.1) then    ! first scan
       let onoff_jy no
       if (c.eq."J") then
         let onoff_jy yes
       endif
      else
       if (onoff_jy.and.c.ne."J") then
         say " "
         say " Incompatible units"
         let onoff_error yes
         return
       else if (.not.onoff_jy.and.c.eq."J") then
         say " "
         say " Incompatible units"
         let onoff_error yes
         return
       endif
      endif
     endif
   endif
 next


 say " "
 say " "
 say "Using procedure @ p_onoff-sum to compute mean, rms"
 say "and snr of the ON-OFF data in file:" 'filename'
 on error return
 col x 1 /file 'filename'
!
! test the channels and observations number
!
 if (nchannel.gt.0) then
         say " "
	 say "Number of channels : " 'nchannel'
	 say "Number of observations found : " 'nmax'
	 say " "
!
!	define the adapted size variable
!
  define double val[nchannel,nmax] sigma[nchannel,nmax]
  if exist(mean) then
   del /var mean  sigma_mean  mean_chan sigma_mean_chan onoff_time
  endif
  define double mean[nchannel] sigma_mean[nchannel] /global 
  define double mean_ext sigma_ext sigma2
  define double mean_corr sigma_corr s_to_n_corr
  define double s_to_n[nchannel]
  define integer info
!
!
! compute the mean and rms value for each channels
!
  for j 1 to nmax
        for i 1 to nchannel
	   let val[i,j]   t1[i,j]
           let sigma[i,j] t2[i,j]
        next
  next j
  del /var t1 t2
  def real mean_chan[nmax,nchannel] sigma_mean_chan[nmax,nchannel] /global
  def real onoff_time[nmax] /global
  for i 1 to nchannel
	let mean[i] = 0.0
	let sigma_mean[i] = 0.0
	let info = 0
	for j 1 to nmax
	  if (sigma[i,j].gt.0.0) then
            let sigma2 = sigma[i,j]**2
	    let mean[i] = mean[i]+val[i,j]|sigma2
	    let sigma_mean[i] = sigma_mean[i]+1.0|sigma2
	    let mean_chan[j,i]  mean[i]|sigma_mean[i]
            let sigma_mean_chan[j,i] sqrt(1.0|sigma_mean[i])
	  else
	    if (info.ne.1) then
	      let info = 1
	      say "WARNING channel " 'i' " : sigma = 0, data not used (Flagged ?)"
	    endif
	  endif
	next j
        if (sigma_mean[i].gt.0) then
	  let mean[i] = mean[i]|sigma_mean[i]
	  let sigma_mean[i] = sqrt(1.0|sigma_mean[i])
	  let s_to_n[i] = abs(mean[i]|sigma_mean[i])
        end if
  next i
  let onoff_time[1] ti[1]|2
  for j 2 to nmax
     let ii j-1
     let onoff_time[j] ti[j]|2+onoff_time[ii]
  next
!  -------------------------------------------------------------
! 
  for i nchannel to 1 by -1
!
!  truncation of the data
!
   say "Channel : " i "    Mean : " mean[i] "    Rms : " sigma_mean[i] "   Snr : " s_to_n[i] /format a i2 a f13.5 a f11.5 a f6.2  
  next i
!
!  compute the mean value of all the external channels, and their rms
!
  let mean_ext = 0.0
  let sigma_ext = 0.0
  for i 1 to nchannel
     if (i.ne.ref_chan) then
	if (sigma_mean[i].gt.0) then
		let sigma2 = sigma_mean[i]**2
		let mean_ext = mean_ext+mean[i]|sigma2
		let sigma_ext = sigma_ext+1.0|sigma2
	end if
     end if
  next i
  if (sigma_ext.gt.0) then
	   let mean_ext = mean_ext|sigma_ext
	   let sigma_ext = sqrt(1.0|sigma_ext)
  end if
!
  say " "
  say "Mean value of the adjacent channels    : " mean_ext "   Rms : " sigma_ext /format a f13.5 a f11.5
  say " "
 else
	say " "
	say "Scan list is empty ..."
 endif
endif