Files
sea/tcl/drivers/zcalib.tcl
2025-02-28 11:27:07 +01:00

443 lines
11 KiB
Tcl

# 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
}
}