#!/bin/sh
# the next line restarts using wish \
exec wish $0 $@

set Home [file join $env(TCLTOOLS_HOME)]
source [ file join $Home TclToolInits.tcl ]
TclToolInits GUI GPH PLT UTILS UDF UDFDB UDFUTILS FFT

# READ a set of satellite potential data. Fit it using a Savitsky-Golay
#   filter.  Plot the data, the fit, and difference between the data 
#   and the fit.

# THIS procedure acquires the data.  The input is the UDF Key, the PIDF
#   sensor and unit number, and the beginning and ending year, day, hour, 
#   minute, and second.  Return the time and data arrays and the time of
#   one measurement.

proc ReadData { KeY sN uN bY bD bH bM bS eY eD eH eM eS tD pD dT } {
   global Pidf0 ExDa

   upvar $tD tM
   upvar $pD dV
   upvar $dT mT

# OPEN the data file

   FileOpen $KeY 0 $bY $bD $bH $bM $bS $eY $eD $eH $eM $eS

# GET to ths start time 

   FilePos $KeY 0 $bY $bD $bH $bM $bS FINE 

# READ the PIDF file

   ReadPIDF $KeY 0 

# GET the converstion algorithm

   set uI [UUunitAlg $KeY $uN $sN]

# GET the VIDF number of the sensor

   set vN [lindex $Pidf0(SenVidfN) $sN]

# THIS is the beginning and ending time in milliseconds

   set bMs [expr $bH * 3600000 + $bM * 60000 + $bS]
   set eMs [expr $eH * 3600000 + $eM * 60000 + $eS]

# THIS is the end time in the required time compare format

   set Tdone [format "%4d%03d%08d" $eY $eD $eMs]

# INITIALIZE looping flag 

   set GoOn 2

# INITIALIZE the index counter 

   set nE 0

# LOOP over the data until we have reached the end time

   while { $GoOn > 0 } { 

# READ the data

      set S [ReadUDF $KeY 0 0 0 0 1] 

# CHECK for out of data

      if { [lindex [split $S "|"] 0] == -1 } { break }
      
# CONVERT the data to units and save.  Save beginning time as well 

      set dV($nE) [ConvertToUnits $KeY 0 0 $uI rV]
      set tM($nE) [list $ExDa(0,BYR) $ExDa(0,BDY) $ExDa(0,BMSEC) 0]

# CHECK to see if the current time is between the requested start and stop
#   time.  If after the stop time we are done.

      set Ts [format "%4d%03d%08d" $ExDa(0,BYR) $ExDa(0,BDY) $ExDa(0,BMSEC)]
      set Te [format "%4d%03d%08d" $ExDa(0,EYR) $ExDa(0,EDY) $ExDa(0,EMSEC)]
      set tList [lsort -ascii -increasing [list $Ts $Tdone $Te]]
      set GoOn [lsearch $tList $Tdone]

# INCREMENT the index counter

      incr nE
   }

# THIS is the interval of time in seconds spaned by the measurement

   set mT [expr ($ExDa(0,EMSEC) - $ExDa(0,BMSEC)) / 1000.0]

# DONE - return the number of data points read

   return $nE
}

# THIS is the time range we are interested in

set bY 2005 ; set bD 33 ; set bH 14 ; set bM 0  ; set bS 0;
set eY 2005 ; set eD 33 ; set eH 14 ; set eM 45 ; set eS 0;

# THIS is the plot time range 

set BT [list $bY $bD $bH $bM $bS 0]
set ET [list $eY $eD $eH $eM $eS 0]

# EXPAND the input time by 5 minutes to take into account the effect of the
#   filter length.  If will corrupt data at the beginning and ending of the 
#   time period  

set sY $bY ; set sD $bD ; set sH $bH ; set sM [expr $bM - 5] ; set sS $bS
set fY $eY ; set fD $eD ; set fH $eH ; set fM [expr $eM + 5] ; set fS $eS

if { $sM < 0 } {
   incr sM 60 ; incr sH -1
   if { $sH < 0 } { incr sH 24 ; incr sD -1 }
   if { $sD < 1 } {
      incr sY -1
      incr sD [TUyearLeap $sY]
   }
}
if { $fM > 59 } {
   incr fM -60 ; incr fH 
   if { $fH > 23 } { incr fH -24 ; incr fD }
   set DiY [TUyearLeap $fY]
   if { $fD > $DiY } {
      incr fY 
      incr fD -$DiY
   }
}

# INITIALIZE the UDF common variables

UDFInits

# ESTABLISH the UDF Key

set KeY [GetDataKey CLUSTERII CLUSTER-3 EFW POTENTIAL SCPOT_SP]

# READ the data 

set nE [ReadData $KeY 0 1 $sY $sD $sH $sM $sS $fY $fD $fH $fM $fS tM pV mT] 

# CONVERT the data times into offsets in seconds from the beginning of the
#   time at which the plot starts and centered on the measurement.

# THIS is the start and stop time the plot and the same for the times
#   over which data is acquired

set bMs [expr $bH * 3600000 + $bM * 60000 + $bS]
set eMs [expr $eH * 3600000 + $eM * 60000 + $eS]
set sMs [expr $sH * 3600000 + $sM * 60000 + $sS]
set fMs [expr $fH * 3600000 + $fM * 60000 + $fS]

set bT [list $bY $bD $bMs 0]
set eT [list $eY $eD $eMs 0]
set sT [list $sY $sD $sMs 0]
set fT [list $fY $fD $fMs 0]

# THIS Base will give the difference in seconds.

set Base 1

# THIS does the conversion 

for  { set I 0 } { $I < $nE } { incr I } {
    set tM($I) [expr [TUtimeConv $bT $tM($I) 0 Base] + $mT / 2.0]
}

# START up the graphics

GraphicsOn TK ETones SWAP

# DEFINE the two plot windows, top one will hold the difference plot and
#   the lower one the data and SV fit

GenWindow 0 0.1 0.1 0.0 0.9 0.49 0.0  0. -8. 0. 7200. -4. 0. X
GenWindow 1 0.1 0.51 0.0 0.9 0.9 0.0  0. -2. 0. 7200.  2. 0. X

# DEFINE the plots as time based in X spanning the requested time

DefinePlot 0 0 -1 -1 X $BT $ET
DefinePlot 1 1 -1 -1 X $BT $ET

# SET up two grid definitions one which spans the plot time and one which
#   which spans the time over which the data was acquired. Both The X axis of the grid
#   spans the time offsets from the start time of the plot

set CT [list $bY $bD $bMs 0]

set gI(0) [TUtimeConv $CT $bT 0 Base]
set gI(1) [TUtimeConv $CT $eT 0 Base]
set gI(2) -1.0
set gI(3) 1.0
set gI(4) [expr int([lindex $PltInfo(0,Pix) 0])]
set gI(5) 1
set gI(6) POINT
set gI(7) POINT
set gI(8) COLUMN
set gI(9) KEEP
set gI(10) NEW:ADD:END
set gI(11) 10000.0
set gI(12) -1.0e31
set gI(13) -1.0e30
set gI(14) 1.0e30
set gI(15) NO
set gI(16) NO

set gIa(0) [TUtimeConv $CT $sT 0 Base]
set gIa(1) [TUtimeConv $CT $fT 0 Base]
set gIa(2) -1.0
set gIa(3) 1.0
set gIa(4) [expr int(($gIa(1) - $gIa(0)) / $mT)]
set gIa(5) 1
set gIa(6) POINT
set gIa(7) POINT
set gIa(8) COLUMN
set gIa(9) KEEP
set gIa(10) NEW:ADD:END
set gIa(11) 10000.0
set gIa(12) -1.0e31
set gIa(13) -1.0e30
set gIa(14) 1.0e30
set gIa(15) NO
set gIa(16) NO

# STICK the input data in the full grid.  Then do a fill on the grid to fill 
#   in any empty grids

set YY(0) 0.0
TUdataGrid $nE 1 tM tM $nE YY YY pV gData Norm gIa
TUgridFill1D gData gIa X HF -1 1

# THIS returns the X coordinates of the grid centers in the full grid

TUgridInfo 32 gIa XX

# STICK the full grid in the plot grid.  Then do a fill on the grid to fill 
#   in any empty grids.  This lops of the extra five minutes on the ends

TUdataGrid $gIa(4) 1 XX XX $gIa(4) YY YY gData pGrid Norm gI
TUgridFill1D pGrid gI X HF -1 1

# PLOT this data

Plot2D 0 Xb Yb Ab pGrid -1 0 LINE [list 250 < -1.0e30]

# Pad the full grid to 2^N, filling with 0.  Need this length in forming
#    the S-G filter

set End 1
while { $End <= $gIa(4) } {
   set End [expr $End * 2]
}
for { set I $gIa(4) } { $I < $End } { incr I } { set gData($I) 0.0 }

# FILTER the data using a S-G filter with 8 points on either side of the
#   center.  Since potential is spin averaged (~4s) this is about 1 min
#   data - less than the 5 mins we added at beginning and end.

set fOps(0) 8
set fOps(1) 8
set fOps(2) 2
TUdataFilter gData $End SAVITZKYGOLAY fOps fDSG

# GRID the filtered data and plot it

TUgridInfo 32 gIa XX
TUdataGrid $gIa(4) 1 XX XX $gIa(4) YY YY fDSG fGrid Norm gI
TUgridFill1D fGrid gI X HF -1 1
set oP [list $GphInfo(Blue) NONE HOLD HOLD 10 *]
Plot2D 0 Xb Yb Ab fGrid -1 0 LINE $oP

# DIFFERENCE filtered and non-filtered potential data and plot it

for { set I 0 } { $I < $gI(4) } { incr I } { 
   set Grid($I) [expr $pGrid($I) - $fGrid($I)]
}

Plot2D 1 Xb Yb Ab Grid -1 0 LINE $oP

# SET plots to have 5 major and minor ticks along Y axis

PLTinfoChg 0 Y TICKS MJNUMBER 5
PLTinfoChg 0 Y TICKS MNNUMBER 5
PLTinfoChg 1 Y TICKS MJNUMBER 5
PLTinfoChg 1 Y TICKS MNNUMBER 5

# DRAW the plot boxes.

PlotOutlines 0 [list X 1.0 POT 1] [list BOT BOT BOTH BOTH]
PlotOutlines 1 [list "" 1.0 "DIF" 1 ] [list NONE BOT BOTH BOTH]
set Lab "START TIME: $bY $bD ${bH}:${bM}:${bS}"
PLTaxisRec 1 Xt NONE AXIS 1.0 0.0 right $Lab