# calibration scripts
#
# in config file (.stick)
#
# makenv calib -driver zcalib
#
# see also: http://linse-c.psi.ch:8080/sample_environment/83

namespace eval zcalib {
  variable buf
}

proc stdConfig::zcalib {} {
  controller syncedprot 30
  prop write stdSct::complete

  obj ZCALIB out -int
  default 3
  prop check zcalib::command
  prop label "calibration"
  prop enum running,pause,continue,stop
  prop buttons start/0/0-2:pause/1/1-3:continue/2/0,2:stop/3/3

  kids calib {
    node tlist out -text
    default "1.38 logstep0.04 310"
    prop check zcalib::set_tlist
    prop label "T list"
    prop help "a list of temperatures, logstep0.04 creates 25 logarithmic steps per decade in between values"
    prop width 80
    
    node tindex out -int
    default 0
    prop label "T index"
    prop write zcalib::set_tindex
    
    node n_temp upd -int
    prop label "number of T's"
    
    node mode out -int
    default 0
    prop enum waiting,measuring
    prop write zcalib::set_mode
    
    node phase out -int
    default 0
    prop write zcalib::set_phase
    
    node scanpos out -int
    default 0
    prop write zcalib::set_scanpos
    prop help "index within one measuring cycle: 0: reference, 1: 1st chan, 2: ref, 3: 2nd chan, ..."
    
    node channels par -text "3 1 2 4 5 6"
    prop help "first channel is calibrated one"
    prop width 20
    
    node basename par -text "/home/l_samenv/sea/calib_scripts/calib_data/calib%s_p%d_c%d.dat"
    prop width 40
    
    node startdate par -text ""
    
    node tolerance par -float 0.005
    prop help "do not measure when (T_Vti - T_set) > tolerance * T_set"
    
    node slopelimit par -float 0.005
    prop help "do not measure before maxslope < slopelimit"
    
    node maxslope upd
    prop help "max of dR/R/window over all channels"
    
    node window par -float 600
    prop help "window for slope calculation"

    node nmean par -int 10
    
  }
}

proc zcalib::get_temp {index} {
  set T 0
  set logN 0
  set mode single
  set tlist [sctval /calib/tlist]
  foreach item $tlist {
    if {[string match logstep* $item]} {
      if {$T == 0} {
        error "logstep must not appear at the beginning"
      }
      set logstep 0
      scan $item logstep%f logstep
      set mode logstep
    } else {
      if {$mode eq "logstep"} {
        if {$logstep == 0} {
          set logstep $item
          continue
        }
        set n [expr abs(round((log10($item) - log10($T)) / $logstep))]
        if {$index < $n - 1} {
          if {$item > $T} {
            set i [expr 1+$index]
          } else {
            set i [expr -1-$index]
          }
          set dig [expr round(2 - log10($logstep))]
          return [format %.${dig}g [expr $T * pow(10, $i * $logstep)]]
        } else {
          set mode single
          if {$n != 0} {
            incr index [expr 1-$n]
          }
        }
      }
      if {$index == 0} {
        return $item
      }
      set T $item
      incr index -1
    }
  }
  return -$index
}

proc zcalib::clear_buf {} {
  variable buf

  foreach c [hval /calib/channels] {
    set buf($c) [list]
    set buf(dif$c) 1
  }
}

proc zcalib::set_tlist {} {
  set n [get_temp 1000]
  incr n 1000
  updateval [sct parent]/n_temp $n
  sct update [sct target]
}

proc zcalib::activate_nv {} {
  if {[result tt tm] < 12 && [result tt set] < 3 && [result nv] == 0} {
    # do not switch control on when already fixed at low T
  } elseif {[result tt set] > [result tt tm] - 10} {
    nv 1
  } else {
    nv 2
  }
}

proc zcalib::set_tindex {} {
  set i [sct target]
  set T [get_temp $i]
  if {$T <= 0} {
    clientput "close valve and stop pump"
    hepump valve 1
    hepump running 0
    calib 3
    return idle
  }
  clear_buf
  sct update $i
  clientput "next_temp $T"
  tt set $T
  res autoscan 1
  activate_nv
  return idle
}

proc zcalib::set_mode {} {
  updateerror /calib/maxslope undefined
  if {[sct target] == 0} {
    res autoscan 1
    activate_nv
  } else {
    res autoscan 0
    nv 0
  }
  clear_buf
  sct update [sct target]
  return idle
}

proc zcalib::set_phase {} {
  sct update [sct target]
  hset [sct parent]/scanpos 0
  return idle
}

proc zcalib::set_scanpos {} {
  sct update [sct target]
  if {[sct target] % 2 != 0} {
    set i [expr [sct target] / 2 + 1]
    set c [lindex [hval /calib/channels] $i]
  } else {
    set c [lindex [hval /calib/channels] 0]
  }
  res s$c/active 1
  return idle
}

proc zcalib::update_chan {value} {
  variable buf

 catch {
  if {[hval /calib] == 1} {
    # how to treat a real pause?
    return
  }
  set time [format %.3f [expr [DoubleTime] - [hgetpropval /calib @basetime]]]
  set c [sct @channel]
  if {$c != [hval /res]} {
    return
  }
  set channels [hval /calib/channels]
  if {[hval /calib/mode] == 0} {
    lappend buf($c) $time
    lappend buf($c) $value
    set i 0
    set vmin 1e30
    set vmax 0
    set idx 0
    set full 0
      
    # determine average
    set sumtime 0
    set sumres 0
    set n 0
    set w [hval /calib/window]
    foreach {t v} $buf($c) {
      if {$t < $time - $w} {
        set idx $i
      } else {
        set sumtime [expr $sumtime + $t]
        set sumres [expr $sumres + $v]
        incr n
      }
      incr i 2
    }
    if {$idx > 0} {
      # cut down buf to window size
      set buf($c) [lrange $buf($c) $idx end]
    }
    if {$n < 2} {
      set dt 0
    } else {
      set meantime [expr $sumtime / double($n)]
      set meanres [expr $sumres / double($n)]
      
      # determine slope
      set sumdt 0
      set sumdr 0
      foreach {t v} $buf($c) {
        set dt [expr ($t - $meantime)]
        set sumdt [expr $sumdt + $dt * $dt]
        set sumdr [expr $sumdr + ($v - $meanres) * $dt]
      }
      # slope is relative R change within window (Ohm/Ohm/10 min.)
      set buf(slope$c) [expr $sumdr / $sumdt / $meanres * [hval /calib/window]]
      set maxslope 0
      foreach ch $channels {
        if {abs($buf(slope$ch)) > abs($maxslope)} {
          set maxslope $buf(slope$ch)
        }
      }
      set dt [expr round($time - [lindex $buf($c) 0])]
    }
    if {$dt > $w / 3} {
      # do not update maxslope before we have a third of the window
      updateval /calib/maxslope [expr abs($maxslope)]
    }
    if {$dt < $w || abs($maxslope) > [hval /calib/slopelimit]} {
      # do not finish waiting mode before we have the full window and the slope limit is reached
      return
    }
    calib phase [hval /calib/phase]
    # switch to measuring mode
    calib mode 1
  } else {
    if {abs([hval /tt] - [hval /tt/set]) / [hval /tt/set] > [hval /calib/tolerance] && [hval /tt] > 1.5} {
      # go to back to waiting mode when out of relative tolerance and not at base
      clientput "tt out of tolerance - WAIT"
      calib mode 0
      return
    }
    set phase [hval /calib/phase]
    if {$phase >= 3} {
      # after 3 measuring cycles go to next point
      calib phase 0
      calib tindex [expr [hval /calib/tindex] + 1]
      calib mode 0
      return
    }
    set scanpos [hval /calib/scanpos]
    set lim [hval /calib/nmean]
    if {$scanpos % 2 == 0} {
      if {$c != [lindex $channels 0]} {
        return
      }
    } else {
      if {$c != [lindex $channels [expr ($scanpos + 1)/ 2]]} {
        return
      }
    }
    lappend buf($c) $value
    if {$scanpos == 0} {
      set lim [expr $lim / 2]
    }
    set len [llength $buf($c)]
    if {$len < $lim} {
      return
    }
    set mean [expr [::tcl::mathop::+ {*}$buf($c)] / double($len)]
    set sum2 0
    set diflog [list]
    foreach r $buf($c) {
      set sum2 [expr $sum2 + pow($r - $mean, 2)]
      lappend diflog [format %.7g [expr $r - $mean]]
    }
    set sigma [expr sqrt($sum2) / double($len)]
    if {$scanpos == 0} {
    } elseif {$scanpos % 2 == 0} {
      set cf [lindex $channels [expr $scanpos/ 2]]
      set buf(ref$cf) $mean
      set buf(ref_sigma$cf) $sigma
      # keep last half of buffer for next time
      set buf($c) [lrange $buf($c) [expr $len / 2] end]
    } else {
      set buf(value$c) $mean
      set buf(sigma$c) $sigma
      set buf(time$c) $time
      set buf($c) [list]
    }
    incr scanpos
    if {$scanpos <= 2 * [llength $channels] - 2} {
      calib scanpos $scanpos
      return
    }
    catch {
        foreach channel $channels {
          if {$channel != [lindex $channels 0]} {
            set filnam [format [hval /calib/basename] [hval /calib/startdate] $phase $channel]
            set exists [file exists $filnam]
            set fil [open $filnam a]
            if {!$exists} {
              puts $fil "#zdat (zcalib by M.Zolliker 2018)"
              puts $fil "#time since midnight, R ref, R meas, sigma, Tset(nominal)"
            }
            puts $fil [format "%.2f	%.9g	%.9g	%.9g" $buf(time$channel) $buf(ref$channel) $buf(value$channel) $buf(sigma$channel) [result tt set]]
            close $fil
            set filnam "[lindex [split [hval /calib/basename] "%"] 0]_[hval /calib/startdate].txt"
            set exists [file exists $filnam]
            set fil [open $filnam a]
            if {!$exists} {
              puts $fil "#zcalib (format 2021 M.Zolliker)"
              puts $fil "#chan, time since midnight, R ref, R meas, sigma, Tset(nominal)"
            }
            puts $fil [format "%d	%.9g	%.9g	%.9g	%.2f" $channel $buf(ref$channel) $buf(value$channel) $buf(sigma$channel) [result tt set] $buf(time$channel)]
            close $fil
          }
        }
    } msg
    if {$msg ne ""} {
      clientput "save error: $msg"
    }
    clear_buf
    incr phase
    calib phase $phase
  }
 } msg
 if {$msg ne ""} {
   clientput "ERROR: $msg"
 }
}

proc zcalib::set_chan {channel} {
  _res updatescript /res/s$channel/raw zcalib::update_chan
  hsetprop /res/s$channel/raw count 10
  hsetprop /res/s$channel/raw dif 1
}

proc zcalib::set_calchan {channel} {
  _res updatescript /res/s$channel/raw zcalib::update_chan
  hsetprop /res/s$channel/raw count 5
  hsetprop /res/s$channel/raw dif 1
}

proc zcalib::install_scripts {} {
  set channels [hval [sct]/channels]
  foreach channel $channels {
    _res updatescript /res/s$channel/raw zcalib::update_chan
  }
  hsetprop /res @calchan [lindex $channels 0]
}

proc zcalib::command {} {
  switch [sct target] {
    0 { #run
      # update n_temp
      hset [sct]/tlist [hval [sct]/tlist]
      set today [clock format [clock seconds] -format "%Y-%m-%d"]
      hupdate [sct]/startdate $today
      hupdate [sct]/tindex 0
      sct @basetime [clock scan $today]
      install_scripts
      nv autoflow/getTemp zcalib::tmts
      hset [sct]/tindex 0
    }
    1 { #pause
      activate_nv
      res autoscan 1
      calib mode 0
    }
    2 { #continue
      nv 0
      calib tindex [hval /calib/tindex]
      calib mode 0
      sct update 0
      install_scripts
      return idle
    }
    3 { # stop
      nv autoflow/getTemp ccu4flow::tmts
      foreach channel [hval [sct]/channels] {
        _res killupdatescript /res/s$channel/raw zcalib::update_chan
      }
    }
  }
  sct update [sct target]
  return idle
}

proc zcalib::tmts {} {
  set tm [silent 1 result tt tm]
  set ts [silent $tm hval /res/s[hgetpropval /res @calchan]]
  if {$ts < $tm} {
    return $ts
  } else {
    return $tm
  }
}