medfilt := function(x, wid)
{
  mm := array(0,wid)
  wh := as_integer(wid/2)

  if(len(x::shape) != 1)
  {  print "Input array to medfilt must be one dimension only!"
     return mm ;
  }
  nelts := x::shape

  mm := array(0,nelts)

  for( ii in 1:nelts)
  {
    j1 := ii - wh;
    if(j1 < 1) j1 := 1
    j2 := ii + wh;
    if(j2 > nelts) j2 := nelts

    nn :=  x[j1:j2]

    mm[ii] := sort(nn)[wh+1]

  }
  mm::shape := nelts

  return mm;
}

# function to calculate the rms of a function x
calcrms := function(x)
{
  nx     := len(x)
  xmean  := sum(x)/nx
  x2mean := sum( x*x ) / nx
  x2rms  :=  x2mean -(xmean*xmean)

  if(x2rms > 0.0) xrms := sqrt(x2rms)
  else            xrms := 0.0;

  return xrms;
}

# clear accumulator
clracc := function(xx, ref nnacc, ref vv)
{

  val nnacc := 0;
  ss := xx::shape

  val vv := array(0,ss)
  val vv::shape := ss
}

vaccum := function(xx, ref nn, ref vv)
{
   val nn +:=1

   val vv := vv + xx
}

vaverg := function(xx, nn, vv)
{
  print "averaging, n=", nn
  return vv/nn;
}

# average median filtered series of scans
afiltscans := function(scan1, scan2, ref avx, ref avy, nwidth)
{
   avxx := array(0,1023)
   avxx::shape := 1023
   nxx  := 0
   avyy := array(0,1023)
   avyy::shape := 1023
   nyy  := 0

   for(ii in scan1:scan2)
   {
     print "doing scan", ii

     signal(ii)
     vsig := uniget('vsig')
     vx := vsig.data.arr[1,]
     vy := vsig.data.arr[2,]

     vxm := medfilt(vx,nwidth)
     vaccum(vxm, nxx, avxx)

     vym := medfilt(vy,nwidth)
     vaccum(vym, nyy, avyy)
   }
   val avx := vaverg( vx, nxx, avxx)
   val avy := vaverg( vy, nyy, avyy)

}

normdd := function(scanno, avx, avy)
{
  
  signal(scanno)
  vsig := uniget('vsig')
  vx := vsig.data.arr[1,]
  vy := vsig.data.arr[2,]

  vsig.data.arr[1,] := vx/avx
  vsig.data.arr[2,] := vy/avy

  uniput('globalscan1', vsig)
}


# function to paste together many scans
# case for scans running 1 to 1023j
spectpaste := function(scan1,scan2, omit=0)
{
  ss2 := scan2;
  if(omit>0) ss2 -:= 1
  npts := (ss2+1-scan1)*896
  newarr := [=]
  newff  := array(0,npts)
  newxx  := array(0,npts)
  newyy  := array(0,npts)
  newarr := [freq=newff, x=newxx, y=newyy]
  jj := npts-895
  
  for(ii in scan1:scan2)
  {
    if(ii == omit) continue;  # skip bad scan

    signal(ii)
    vsig := uniget('vsig')

    # get array of frequencies: convert to MHZ
    ff := vsig.data.desc.chan_freq.value / 1.0e6 ;

    xx := vsig.data.arr[1,]
    yy := vsig.data.arr[2,]

    # edit out central freq peak
    xx[512] := 0.5*(xx[511] + xx[513])
    yy[512] := 0.5*(yy[511] + yy[513])

    # drop 64 pixels from the front, 63 from the rear.
    # this is for spacing 40MHz spectra at 35 MHz intervals
    # cut spectra each have 896 points

    ffcut := ff[65:960]
    xxcut := xx[65:960]
    yycut := yy[65:960]

   # insert points from rear of array
    newarr.freq[jj:(jj+895)] := ffcut
    newarr.x[jj:(jj+895)]    := xxcut
    newarr.y[jj:(jj+895)]    := yycut

    jj -:= 896;
    
  }
  return newarr;
}


# purge function: works only on 1-D arrays
# removes really big RFIs
cpurge := function(x, x1, x2)
{
  xout := x;

  for(i in 1:len(x))
  {  if( x[i] > x2) xout[i] := 0;
     if( x[i] < x1) xout[i] := 0;
  }
  return xout;
}


# function to paste together many scans
# after normalizing each by dividing by 
# median-filtered version of itself
# case for scans running 1 to 1023
filtpaste := function(scan1,scan2, wfilter, omit=0)
{
  ss2 := scan2;
  if(omit>0) ss2 -:= 1
  npts := (ss2+1-scan1)*896
  newarr := [=]
  newff  := array(0,npts)
  newxx  := array(0,npts)
  newyy  := array(0,npts)

  newarr := [freq=newff, x=newxx, y=newyy]

  jj := npts-895
  
  for(ii in scan1:scan2)
  {
    if(ii == omit) continue;

    signal(ii)
    vsig := uniget('vsig')

    # get array of frequencies: convert to MHZ
    ff := vsig.data.desc.chan_freq.value / 1.0e6 ;

    xx := vsig.data.arr[1,]
    xx::shape := len(xx)

    yy := vsig.data.arr[2,]
    yy::shape := len(yy)

    # edit out central freq peak
    xx[512] := 0.5*(xx[511] + xx[513])
    yy[512] := 0.5*(yy[511] + yy[513])

    # edit out peaks at 256 and 768
    xx[256] := 0.5*(xx[255] + xx[257])
    yy[256] := 0.5*(yy[255] + yy[257])
    xx[768] := 0.5*(xx[767] + xx[769])
    yy[768] := 0.5*(yy[767] + yy[769])

    # do filters

    xxfilt := medfilt(xx, wfilter)
    yyfilt := medfilt(yy, wfilter)

    xx /:= xxfilt
    yy /:= yyfilt

    # drop 64 pixels from the front, 63 from the rear.
    # this is for spacing 40MHz spectra at 35 MHz intervals
    # cut spectra each have 896 points

    ffcut := ff[65:960]
    xxcut := xx[65:960]
    yycut := yy[65:960]

    xxfiltcut := xxfilt[65:960]
    yyfiltcut := yyfilt[65:960]

   # insert points from rear of array
    newarr.freq[jj:(jj+895)] := ffcut
    newarr.x[jj:(jj+895)]    := xxcut
    newarr.y[jj:(jj+895)]    := yycut

    jj -:= 896;
    
  }

  return newarr;
}


# function to average together a list of 
# pasted spectra that were saved with write_value
rfave := function( rflist )
{
   nn := len(rflist)

   for( ii in 1:nn)
   {
     fname := spaste( rflist[ii], '.data')
     print "Adding", ii, fname
     xx := read_value( fname )

     if(ii==1)    # put first value into accumulator
     {  newxx := xx;
     }
     else
     {  newxx.x +:= xx.x;
	newxx.y +:= xx.y;
     }
   }
   newxx.x /:= nn
   newxx.y /:= nn

   return newxx;
}


# clip function: works only on 1-D arrays
cclip := function(x, x1, x2)
{
  xout := x;

  for(i in 1:len(x))
  {  if( x[i] > x2) xout[i] := x2;
     if( x[i] < x1) xout[i] := x1;
  }
  return xout;
}

# make rfilist from array x, using xt as the threshold
# if rfi exceeds limit in successive bins, the width (w)
# is incremented
rfilist := function(x, xt, nneg=1)
{
  nx := len(x.x)
  print "rfilist:", nx, xt

  rf := [=]
  rf := [ f=[0], s=[0], w=[0] ]
  jj := 0
  lastii := -10

  for(ii in 1:nx)
  {
    px := nneg * x.x[ii];
    py := nneg * x.y[ii];

    if(( px >xt)||( py >xt))
    {   
	if((ii-lastii) < 2)
	{ rf.w[jj] +:= 1   # increment width here
	}
	else
	{ jj +:= 1
	  rf.w[jj] := 1;
	  rf.s[jj] := 0
	  rf.f[jj] := 0
	}
	rf.f[jj] := x.freq[ii]
	rf.s[jj] := max(px, rf.s[jj])
	rf.s[jj] := max(py, rf.s[jj])
	lastii := ii
    
    }
  }
  return rf;
}

# will print to file if fname=filename
listrfilist := function(rr, fname='F', skale=0)
{
  print "RFI list", len(rr.f)
  print " Freq(MHz) P(milli) W(chan) "

  if(fname != 'F') 
  {  outf := open(paste(">", fname))
     write(outf,"RFI list")
     write(outf, " Freq(MHz)  P(milli)  W(chan) ")
  }

  nx := len(rr.f)

  for(ii in 1:nx)
  {
    P := 1000.0*(rr.s[ii] - skale)

    st := sprintf(" %8.2f   %7.1f   %5d  ",
	rr.f[ii], P, rr.w[ii]);

    if(fname != 'F') write(outf,st)
    else print st
  }
  outf := F

}

# process to repeat processing of all the filtered spectra.
pf30proc := function()
{
#  dopen('rs918_SP')
#  lscans()
#  print "filtering az096a"
#   zz := filtpaste(8,33,7)
#   write_value(zz, "az096afilt.data")

#  print "filtering az186a"
#   zz := filtpaste(36,61,7)
#   write_value(zz, "az186afilt.data")

#  dopen('rf276_SP')
#  lscans()
#  print "filtering az276a"
#   zz := filtpaste(64,89,7)
#   write_value(zz, "az276afilt.data")

#  print "filtering az006a"
#   zz := filtpaste(92,117,7)
#   write_value(zz, "az006afilt.data")

#  dopen('rf918b_SP')
#  lscans()
#  print "filtering az096b"
#   zz := filtpaste(120,145,7)
#   write_value(zz, "az096bfilt.data")

#  print "filtering az186b"
#   zz := filtpaste(148,173,7)
#   write_value(zz, "az186bfilt.data")

  dopen('rf276b_SP')
  lscans()
  print "filtering az276b"
   zz := filtpaste(176,201,7)
   write_value(zz, "az276bfilt.data")

  print "filtering az006b"
   zz := filtpaste(204,229,7)
   write_value(zz, "az006bfilt.data")

 # now make grand average of all spectra
   rlist := ['az096afilt','az186afilt','az276afilt','az006afilt','az096bfilt','az186bfilt','az276bfilt','az006bfilt']

   cc := rfave(rlist)
   write_value(cc, "filtallaz.data")

  # make average of 96 degree data
   slist := ['az096afilt', 'az096bfilt']
   dd96 := rfave(slist)
   write_value(dd96, "filt096az.data")

  # make average of 186 degree data
   slist := ['az186afilt', 'az186bfilt']
   dd186 := rfave(slist)
   write_value(dd186, "filt186az.data")

  # make average of 276 degree data
   slist := ['az276afilt', 'az276bfilt']
   dd276 := rfave(slist)
   write_value(dd276, "filt276az.data")

  # make average of 006 degree data
   slist := ['az006afilt', 'az006bfilt']
   dd006 := rfave(slist)
   write_value(dd006, "filt006az.data")
   
  # difference each direction by the grand average:

   print "differencing 96 degree spectrum"
   diff96 := dd96;
   diff96.x := dd96.x - cc.x
   diff96.y := dd96.y - cc.y
   write_value(diff96, "diff096.data")

   print "differencing 186 degree spectrum"
   diff186 := dd186;
   diff186.x := dd186.x - cc.x
   diff186.y := dd186.y - cc.y
   write_value(diff186, "diff186.data")

   print "differencing 276 degree spectrum"
   diff276 := dd276;
   diff276.x := dd276.x - cc.x
   diff276.y := dd276.y - cc.y
   write_value(diff276, "diff276.data")

   print "differencing 006 degree spectrum"
   diff006 := dd006;
   diff006.x := dd006.x - cc.x
   diff006.y := dd006.y - cc.y
   write_value(diff006, "diff006.data")

}


# process to do processing of all the Feb 7-8 spectra.
# for the Feb7-8 run.
pf07proc := function()
{
#  dopen('r918a_SP')
#  lscans()
#   print "pasting az096a"
#   zz := spectpaste(22,47)
#   write_value(zz, "raw096a.data")
#   ff := filtpaste(22,47,7)
#   write_value(ff, "filt096a.data")

#  print "pasting az186a"
#   zz := spectpaste(50,76, omit=66)
#   write_value(zz, "raw186a.data")
#   ff := filtpaste(50,76,7, omit=66)
#   write_value(ff, "filt186a.data")

#   dopen('r206a_SP')
#   lscans()
#   print "pasting az276a"
#   zz := spectpaste(79,104)
#   write_value(zz, "raw276a.data")
#   print "filtering az276a"
#   ff := filtpaste(79,104,7)
#   write_value(ff, "filt276a.data")

#   print "pasting az006a"
#    zz := spectpaste(107,132)
#    write_value(zz, "raw006a.data")
#   print "filtering az006a"
#    ff := filtpaste(107,132,7)
#    write_value(ff, "filt006a.data")

#   dopen('r918b_SP')
#   lscans()
#   print "pasting az096b"
#   zz := spectpaste(135,160)
#   write_value(zz, "raw096b.data")
#   print "filtering az096b"
#   ff := filtpaste(135,160,7)
#   write_value(ff, "filt096b.data")
#
#   print "pasting az186b"
#    zz := spectpaste(163,188)
#    write_value(zz, "raw186b.data")
#   print "filtering az186b"
#    ff := filtpaste(163,188,7)
#    write_value(ff, "filt186b.data")

#  dopen('r206b_SP')
#  lscans()
#  print "pasting az276b"
#  zz := spectpaste(191,216)
#  write_value(zz, "raw276b.data")
#  print "filtering az276b"
#  ff := filtpaste(191,216,7)
#  write_value(ff, "filt276b.data")

#  print "pasting az006b"
#   zz := spectpaste(219,244)
#   write_value(zz, "raw006b.data")
#  print "filtering az006b"
#   ff := filtpaste(219,244,7)
#   write_value(ff, "filt006b.data")

#  print "pasting az096c"
#   zz := spectpaste(247,272)
#   write_value(zz, "raw096c.data")
#  print "filtering az096c"
#   ff := filtpaste(247,272,7)
#   write_value(ff, "filt096c.data")

 # now make grand average of all raw spectra
 # leave out 096c

#  rwlist := ['raw096a', 'raw186a', 'raw276a', 'raw006a', 'raw096b',
#'raw186b', 'raw276b', 'raw006b']
#  ccw := rfave(rwlist)
#  write_value(ccw, "rawallaz.data")

 # now make grand average of all filtered spectra

#  rlist := ['filt096a','filt186a','filt276a','filt006a','filt096b','filt186b','filt276b','filt006b']

#  cc := rfave(rlist)
#  write_value(cc, "filtallaz.data")

  # make average of 96 degree data
#  slist := ['filt096a', 'filt096b']
#  dd96 := rfave(slist)
#  write_value(dd96, "avfilt096az.data")

  # make average of 186 degree data
#  slist := ['filt186a', 'filt186b']
#  dd186 := rfave(slist)
#  write_value(dd186, "avfilt186az.data")

  # make average of 276 degree data
#  slist := ['filt276a', 'filt276b']
#  dd276 := rfave(slist)
#  write_value(dd276, "afvilt276az.data")

  # make average of 006 degree data
#  slist := ['filt006a', 'filt006b']
#  dd006 := rfave(slist)
#  write_value(dd006, "avfilt006az.data")
   
  # difference each direction by the grand average:

#  print "differencing 96 degree spectrum"
#  diff96 := dd96;
#  diff96.x := dd96.x - cc.x
#  diff96.y := dd96.y - cc.y
#  write_value(diff96, "diff096.data")

#  print "differencing 186 degree spectrum"
#  diff186 := dd186;
#  diff186.x := dd186.x - cc.x
#  diff186.y := dd186.y - cc.y
#  write_value(diff186, "diff186.data")

#  print "differencing 276 degree spectrum"
#  diff276 := dd276;
#  diff276.x := dd276.x - cc.x
#  diff276.y := dd276.y - cc.y
#  write_value(diff276, "diff276.data")

#  print "differencing 006 degree spectrum"
#  diff006 := dd006;
#  diff006.x := dd006.x - cc.x
#  diff006.y := dd006.y - cc.y
#  write_value(diff006, "diff006.data")
   
  # get difference: lasers on-off: filtered data
# print "laser on-off: filtered data"
# onfilt := read_value("onfilt.data")
# offilt := read_value("filtallaz.data")
# diffonoff := offilt;
# diffonoff.x := onfilt.x - offilt.x
# diffonoff.y := onfilt.y - offilt.y

# write_value(diffonoff, "diffonoffilt.data")

# diffonoff.x := (onfilt.x - offilt.x)/(onfilt.x + offilt.x)
# diffonoff.y := (onfilt.y - offilt.y)/(onfilt.y + offilt.y)
# write_value(diffonoff, "diffonoffiltnorm.data")
  
  # get difference: lasers on-off: raw data
  print "laser on-off: raw data"
  onfilt := read_value("onraw.data")
  offilt := read_value("rawallaz.data")
  diffonoff := offilt;
  diffonoff.x := onfilt.x - offilt.x
  diffonoff.y := onfilt.y - offilt.y
  write_value(diffonoff, "diffonofraw.data")
 
  diffonoff.x := (onfilt.x - offilt.x)/(onfilt.x + offilt.x)
  diffonoff.y := (onfilt.y - offilt.y)/(onfilt.y + offilt.y)
  write_value(diffonoff, "diffonofrawnorm.data")
}

pplt := function(xx)
{
  pg.clear()
  pg.plotxy1(xx.freq, xx.x)
  pg.plotxy1(xx.freq, xx.y)
}

plot_thing := function(name)
{
  xx := read_value(name)

  pplt(xx)
}