temperature.py

   1 import Melcor
   2 import time
   3 import stripchart
   4
   5 VERSION = "0.2"
   6
   7 # buzzwords: 'integrator windup' for integral term built up during a slow approach.
   8
   9 class error (Exception) :
  10     "Errors with the temperature controller"
  11     pass
  12
  13 class errorMelcor (error) :
  14     pass
  15 class errorOutOfRange (error) :
  16     pass
  17
  18 def _check1(functionCall) :
  19     (err, val) = functionCall
  20     if err != 0 :
  21         raise errorMelcor
  22     return val
  23
  24 def _check0(functionCall) :
  25     err = functionCall
  26     if err != 0 :
  27         raise errorMelcor
  28
  29
  30 def melcor2double(value) :
  31     (err, doub) = Melcor.melcor2double(value)
  32     if err != 0 :
  33         raise errorMelcor, "Error converting melcor to double"
  34     return doub
  35 def double2melcor(doub) :
  36     (err, val) = Melcor.double2melcor(doub)
  37     if err != 0 :
  38         raise errorMelcor, "Error converting double to melcor"
  39     return val
  40
  41 def check_range(raw_output, min, max) :
  42     if raw_output < min :
  43         raise errorOutOfRange, '%g < %g' % (raw_output, min)
  44     if raw_output > max :
  45         raise errorOutOfRange, '%g > %g' % (raw_output, max)
  46
  47 class tempController :
  48     "Pretty wrappers for controlling a Melcor MTCA Temperature Controller"
  49     def __init__(self, controller=1, device='/dev/ttyS0', maxCurrent=0.2) :
  50         """
  51         (controller, device, maxCurrent) -> (tempController instance)
  52         controller : MTCA controller Id
  53         device     : serial port you're using to connect to the controller
  54         maxCurrent : initial maximum allowed current (in Amps)
  55         Set maxCurrent = None if you don't want to adjust from it's prev. value.
  56
  57         0.2 A is the default max current since it seems ok to use without fluid
  58         cooled heatsink.  If you are cooling the heatsink, use 1.0 A, which seems
  59         safely below the peltier's 1.2 A limit.
  60         """
  61         self.verbose = False
  62         self.setpoint = 20.0 # degrees C
  63         self.Tmin = 5.0 # setup some protective bounds for sanity checks
  64         self.Tmax = 50.0
  65         self.specMaxCur = 4.0 # Amps, the rated max current from controller specs
  66         self.T = Melcor.tempController(controller, device)
  67         self.Tstrip = stripchart.stripchart(pipename='Tstrip_pipe',
  68                                             title='Temp strip')
  69         self.Cstrip = stripchart.stripchart(pipename='Cstrip_pipe',
  70                                             title='Current strip')
  71         if maxCurrent != None : # if None, just leave maxCurrent at it's prev. val.
  72             self.setMaxCurrent(maxCurrent) # Amps
  73     def getTemp(self) :
  74         "Returns the current process temperature in degrees Celsius"
  75         val = self.read(Melcor.REG_HIGH_RESOLUTION)
  76         temp = val/100.0
  77         if self.Tstrip.status == 'open' :
  78             self.Tstrip.add_point(temp)
  79         return temp
  80     def getAmbientTemp(self) :
  81         "Returns room temperature in degrees Celsius"
  82         val = self.read(Melcor.REG_AMBIENT_TEMPERATURE)
  83         # convert (Fahrenheit*10) to Celsius
  84         return (val/10.0 - 32)/1.8
  85     def setSetpoint(self, setpoint) :
  86         "Set the temperature setpoint in degrees Celsius"
  87         val = double2melcor(setpoint)
  88         self.write(Melcor.REG_SET_POINT_1, val)
  89     def getSetpoint(self) :
  90         "Get the temperature setpoint in degrees Celsius"
  91         val = self.read(Melcor.REG_SET_POINT_1)
  92         return melcor2double(val)
  93     def setMaxCurrent(self, maxCur) :
  94         """
  95         Set the max current in Amps.
  96         (Note to Melcor enthusiasts: set's both the 'above' and 'below' limits)
  97         """
  98         maxPercent = maxCur / self.specMaxCur * 100
  99         val = double2melcor(maxPercent)
 100         self.write(Melcor.REG_HIGH_POWER_LIMIT_ABOVE, val)
 101         self.write(Melcor.REG_HIGH_POWER_LIMIT_BELOW, val)
 102         self.maxCurrent = maxCur
 103     def getMaxCurrent(self) :
 104         """
 105         () -> (currentLimitAbove (A), currentLimitBelow (A), currentLimitSetpoint (deg C))
 106         """
 107         per = self.read(Melcor.REG_HIGH_POWER_LIMIT_ABOVE)
 108         curLimAbove = melcor2double(per)/100.0 * self.specMaxCur
 109         per = self.read(Melcor.REG_HIGH_POWER_LIMIT_BELOW)
 110         curLimBelow = melcor2double(per)/100.0 * self.specMaxCur
 111         val = self.read(Melcor.REG_POWER_LIMIT_SETPOINT)
 112         curLimSet = melcor2double(val)
 113         return (curLimAbove, curLimBelow, curLimSet)
 114     def getPercentCurrent(self) :
 115         """
 116         Returns the percent of rated max current being output.
 117         See getCurrent()
 118         """
 119         val = int(self.read(Melcor.REG_PERCENT_OUTPUT))
 120         if val > 2**15 :
 121             val -= 2**16
 122         return float(val)/10.0
 123     def getCurrent(self) :
 124         """
 125         The returned current is not the actual current,
 126         but the current that the temperature controller
 127         calculates it should generate.
 128         If the voltage required to generate that current
 129         exceeds the controllers max voltage (15V on mine),
 130         then the physical current will be less than the
 131         value returned here.
 132         """
 133         percentOutput = self.getPercentCurrent()
 134         cur = self.specMaxCur * percentOutput / 100.0
 135         if self.Cstrip.status == 'open' :
 136             self.Cstrip.add_point(cur)
 137         return cur
 138     def setCoolingGains(self, propband=0.1, integral=0, derivative=0) :
 139         """
 140         (propband, integral, derivative, dead_band) -> None
 141         propband   : propotional gain band in degrees C
 142         integral   : integral weight in minutes (0.00 to 99.99)
 143         derivative : derivative weight in minutes (? to ?)
 144         See 5.10 and the pages afterwards in the manual for Melcor's explaination.
 145         Formula (from Cornell BioPhys El Producto Beamline notes)
 146         P_cout = -1/T_prop * [ (T_samp - T_set)
 147                                + 1/t_int * int_-inf^t (T_samp(t')-T_set(t')) dt'
 148                                + t_deriv * dT_samp/dt
 149         Where P_cout is the percent of the rated max current that the controller
 150          would like to output if you weren't limiting it,
 151         T_prop is the propband input to this function,
 152         T_samp is the measured temperature of the sample in deg C,
 153         T_set is the setpoint in deg C,
 154         t_int is the integral input to this function,
 155         the integral with respect to t' is actually only from the time that
 156          T_samp has been with T_prop of T_set (not -inf), and
 157         t_deriv is the derivative input to this function.
 158
 159         Cooling is output 1
 160         """
 161         check_range(propband, 0, 99.9)
 162         check_range(integral, 0, 99.99)
 163         check_range(derivative, 0, 9.99)
 164
 165         val = double2melcor(propband)
 166         self.write(Melcor.REG_PROPBAND_1, val)
 167         val = int(integral * 100)
 168         self.write(Melcor.REG_INTEGRAL_1, val)
 169         val = int(derivative * 100)
 170         self.write(Melcor.REG_DERIVATIVE_1, val)
 171     def getCoolingGains(self) :
 172         "() -> (propband, integral, derivative)"
 173         val = self.read(Melcor.REG_PROPBAND_1)
 174         propband = melcor2double(val)
 175         val = self.read(Melcor.REG_INTEGRAL_1)
 176         integral = val/100.0
 177         val = self.read(Melcor.REG_DERIVATIVE_1)
 178         derivative = val/100.0
 179         return (propband, integral, derivative)
 180     def setHeatingGains(self, propband=0.1, integral=0, derivative=0) :
 181         """
 182         (propband, integral, derivative, dead_band) -> None
 183         propband   : propotional gain band in degrees C
 184         integral   : integral weight in minutes (0.00 to 99.99)
 185         derivative : derivative weight in minutes (? to ?)
 186         Don't use derivative, dead time.
 187         Cycle time?
 188         Histerysis?
 189         Burst?
 190         See 5.10 and the pages afterwards in the manual for Melcor's explaination.
 191         Formula (from Cornell BioPhys El Producto Beamline notes)
 192         P_cout = -1/T_prop * [ (T_samp - T_set)
 193                                + 1/t_int * int_-inf^t (T_samp(t')-T_set(t')) dt'
 194                                + t_deriv * dT_samp/dt
 195         Where P_cout is the percent of the rated max current that the controller
 196          would like to output if you weren't limiting it,
 197         T_prop is the propband input to this function,
 198         T_samp is the measured temperature of the sample in deg C,
 199         T_set is the setpoint in deg C,
 200         t_int is the integral input to this function,
 201         the integral with respect to t' is actually only from the time that
 202          T_samp has been with T_prop of T_set (not -inf), and
 203         t_deriv is the derivative input to this function.
 204
 205         Heating is output 2
 206         """
 207         check_range(propband, 0, 99.9)
 208         check_range(integral, 0, 99.99)
 209         check_range(derivative, 0, 9.99)
 210
 211         val = double2melcor(propband)
 212         self.write(Melcor.REG_PROPBAND_2, val)
 213         val = int(integral * 100)
 214         self.write(Melcor.REG_INTEGRAL_2, val)
 215         val = int(derivative * 100)
 216         self.write(Melcor.REG_DERIVATIVE_2, val)
 217     def getHeatingGains(self) :
 218         "() -> (propband, integral, derivative)"
 219         val = self.read(Melcor.REG_PROPBAND_2)
 220         propband = melcor2double(val)
 221         val = self.read(Melcor.REG_INTEGRAL_2)
 222         integral = val/100.0
 223         val = self.read(Melcor.REG_DERIVATIVE_2)
 224         derivative = val/100.0
 225         return (propband, integral, derivative)
 226     def getFeedbackTerms(self) :
 227         """
 228         Experimental
 229         """
 230         pid = int(self.read(Melcor.REG_PID_POWER_1))
 231         if pid > 2**15 :
 232             pid -= 2**16
 233         prop = int(self.read(Melcor.REG_PROP_TERM_1))
 234         if prop >= 2**15 :
 235             prop -= 2**16
 236         ntgrl = int(self.read(Melcor.REG_INTEGRAL_TERM_1))
 237         print ntgrl
 238         if ntgrl >= 2**15 :
 239             ntgrl -= 2**16
 240         deriv = int(self.read(Melcor.REG_DERIVATIVE_TERM_1))
 241         if deriv >= 2**15 :
 242             deriv -= 2**16
 243         return (pid, prop, ntgrl, deriv)
 244     def checkFeedbackTerms(self) :
 245         pid, prop, ntgrl, deriv = self.getFeedbackTerms()
 246         pout = self.getPercentCurrent()
 247         T = self.getTemp()
 248         Tset = self.getSetpoint()
 249         if T > Tset : # cooling
 250             Tprop, tint, tderiv = self.getCoolingGains()
 251         else : # heating
 252             Tprop, tint, tderiv = self.getHeatingGains()
 253         print "pid(read) %g =? sum(calc from terms) %g =? cur(read) %g" % (pid, prop+ntgrl+deriv, pout)
 254         print "read:     prop %d, integral %d, deriv %d" % (prop, ntgrl, deriv)
 255         print "my calcs: prop %g" % ((Tset-T)/Tprop)
 256     def setTemp(self, setpoint, tolerance=0.3, time=10.0) :
 257         """
 258         Changes setpoint to SETPOINT and waits for stability
 259         """
 260         self.setSetpoint(setpoint)
 261         while self.isStable(setpoint, tolerance, time) != True :
 262             pass
 263     def setTemp_funkygain(self, setpoint, dead_time, heat_rate, cool_rate,
 264                           peltier_efficiency_fn, outside_equilib_rate,
 265                           tolerance=0.3, time=10.0) :
 266         """
 267         Highly experimental, see diffusion.py
 268         """
 269         mode = ""
 270         T = self.getTemp()
 271         # full steam ahead
 272         print "full steam ahead"
 273         self.setSetpoint(setpoint)
 274         self.setHeatingGains(0.1, 0, 0)
 275         self.setCoolingGains(0.1, 0, 0)
 276         if T < setpoint :
 277             mode = "Heating"
 278             self._heat_until_close(setpoint, dead_time, heat_rate)
 279         elif T > setpoint :
 280             mode = "Cooling"
 281             self._cool_until_close(setpoint, dead_time, cool_rate)
 282         # coast
 283         print "coast while temperature equilibrates"
 284         self.setHeatingGains(100, 0, 0)
 285         self.setCoolingGains(100, 0, 0)
 286         time.sleep(dead_time*2)
 287         cool_prop, heat_prop = self.calcPropBands()
 288         print "calculated prop bands: c %g, h %g deg C" % (cool_prop, heat_prop)
 289         print "reset integral gain, and bump to predicted props"
 290         # pop down to reset integral gain, could also jump setpoint...
 291         self.setHeatingGains(0.1, 0, 0)
 292         self.setHeatingGains(heat_prop, 0, 0)
 293         self.setCoolingGains(0.1, 0, 0)
 294         self.setCoolingGains(cool_prop, 0, 0)
 295         time.sleep(dead_time*4)
 296         # now add in some integral to reduce droop
 297         print "set integral gains to %g" % (dead_time*4)
 298         self.setHeatingGains(heat_prop, dead_time*4, 0)
 299         self.setCoolingGains(cool_prop, dead_time*4, 0)
 300         time.sleep(dead_time*8)
 301         print "wait to enter tolerance band"
 302         while (self.getTemp()-setpoint) :
 303             time.sleep(dead_time)
 304         print "should be stable now"
 305         if not self.isStable(setpoint, tolerance, time) :
 306             raise error, "Algorithm broken ;)"
 307     def _heat_until_close(self, setpoint, dead_time, heat_rate) :
 308         while self.getTemp() < setpoint - 0.5*rate*dead_time :
 309             time.sleep(dead_time/10.0)
 310     def calcPropBands_HACK(setpoint, peltier_efficiency_fn, outside_equilib_rate) :
 311         heat_loss = outside_equilib_rate * (setpoint - self.getAmbientTemp())
 312         required_current = heat_loss / peltier_efficiency_fn(setpoint)
 313         if required_current > self.maxCurrent :
 314             raise errorOutOfRange, "Can't source %g Amps", required_current
 315         fraction_current = required_current / self.maxCurrent
 316         droop = 0.5 # expected droop in deg C on only proporitional gain
 317         # droop / T_prop = fraction current
 318         T_prop = droop / fraction_current
 319         if setpoint > self.getAmbientTemp()+5 : # heating
 320             return (T_prop*10, T_prop)
 321         elif setpoint < self.getAmbientTemp()+5 : # cooling
 322             return (T_prop, T_prop*10)
 323         else : # right about room temperature
 324             return (T_prop, T_prop)
 325     def getMode(self) :
 326         mcode = self.read(Melcor.REG_AUTO_MANUAL_OP_MODE)
 327         if mcode == 0 :
 328             return 'auto'
 329         elif mcode == 1 :
 330             return 'manual'
 331         else :
 332             raise error, "Unrecognized mode code %d" % mcode
 333     def setMode(self, mode) :
 334         if mode == 'auto' :
 335             mcode = 0
 336         elif mode == 'manual' :
 337             mcode = 1
 338         else :
 339             raise error, "Unrecognized mode %s" % mode
 340         self.write(Melcor.REG_AUTO_MANUAL_OP_MODE, mcode)
 341     def getManualCurrent(self) :
 342         val = int(self.read(Melcor.REG_MANUAL_SET_POINT))
 343         if val > 2**15 : # convert to signed
 344             val -= 2**16
 345         pct = float(val)/10.0 # stored value is percent * 10
 346         return self.specMaxCur * pct / 100.0
 347     def setManualCurrent(self, amps) :
 348         if amps > self.maxCurrent :
 349             raise error, "Suggested current %g > max %g" % \
 350                 (amps, self.maxCurrent)
 351         pct = amps / self.specMaxCur * 100.0
 352         val = int(pct * 10.0)
 353         if val < 0 : # convert to unsigned
 354             val += 2**16
 355         self.write(Melcor.REG_MANUAL_SET_POINT, val)
 356     def calcPropBands_ZN_get_step_response(self,
 357                                            initial_current=None,
 358                                            initial_wait_time=20.0,
 359                                            current_step=0.1,
 360                                            response_wait_time=40.0,
 361                                            plotVerbose=False) :
 362         """
 363         Ziegler-Nichols tuning, using step response for input.
 364         Process must be stable when calling this function.
 365         """
 366         if initial_current == None :
 367             initial_current = self.getCurrent()
 368         original_mode = self.getMode()
 369         if original_mode == 'manual' :
 370             original_current = self.getManualCurrent()
 371         else :
 372             self.setMode('manual')
 373         self.setManualCurrent(initial_current)
 374         Tarr = []
 375         tarr = []
 376         start = time.time()
 377         tm = start
 378         # get some stability data before stepping
 379         while tm < start + initial_wait_time :
 380             tarr.append(tm-start)
 381             Tarr.append(self.getTemp())
 382             tm = time.time()
 383         # step the output
 384         self.setManualCurrent(initial_current + current_step)
 385         Tlast = Tarr[-1]
 386         while tm < start + initial_wait_time + response_wait_time :
 387             Tnow = self.getTemp()
 388             if Tnow != Tlast : # save us some trouble averaging later
 389                 tarr.append(tm-start)
 390                 Tarr.append(self.getTemp())
 391             Tlast = Tnow
 392             tm = time.time()
 393
 394         self.setMode(original_mode)
 395         if original_mode == 'manual' :
 396             self.setManualCurrent(original_current)
 397         if plotVerbose == True :
 398             from pylab import figure, plot, title, xlabel, ylabel
 399             figure(10)
 400             plot(tarr, Tarr, 'r.-')
 401             xlabel('time (s)')
 402             ylabel('Temp (C)')
 403             title('Plant step response')
 404         return (tarr, Tarr)
 405     def calcPropBands_ZN_analyze_step_response(self, tarr, Tarr,
 406                                                initial_wait_time=20,
 407                                                textVerbose=False) :
 408         """
 409         Analyze the step response to determine dead time td,
 410         and slope at point of inflection spoi.
 411         """
 412         i=0
 413         while tarr[i] < initial_wait_time :
 414             i += 1
 415         # find point of inflection (steepest slope)
 416         istep = i
 417         ipoi = istep
 418         def slope(i) :
 419             return (Tarr[i+1] - Tarr[i])/(tarr[i+1]-tarr[i])
 420         for i in range(istep, len(tarr)-1) :
 421             if slope(i) > slope(ipoi) :
 422                 ipoi = i
 423         print "Max slope at t = %g, T = %g" % (tarr[ipoi], Tarr[ipoi])
 424         spoi = slope(ipoi)
 425         # find the dead time
 426         # find the initial temperature
 427         initialT = 0.0
 428         for i in range(istep) :
 429             initialT += Tarr[i]
 430         initialT /= float(istep)
 431         print "Initial temperature %g" % initialT
 432         deltaT = Tarr[ipoi] - initialT
 433         rise_t_to_poi = spoi/deltaT
 434         td = tarr[ipoi] - initial_wait_time - rise_t_to_poi
 435         return (td, spoi)
 436     def calcPropBands_ZN_compute_terms(self, td, spoi) :
 437         kp = 1.2/spoi
 438         ki = 2*td
 439         kd = td/2
 440         return (kp, ki, kd)
 441     def calcPropBands_ZN(self, initial_current=None,
 442                          initial_wait_time=20.0,
 443                          current_step=0.1,
 444                          response_wait_time=40.0,
 445                          textVerbose=False,
 446                          plotVerbose=False) :
 447         tarr, Tarr = self.calcPropBands_ZN_get_step_response( \
 448             initial_current, initial_wait_time,
 449             current_step, response_wait_time,
 450             plotVerbose)
 451         td, spoi = self.calcPropBands_ZN_analyze_step_response( \
 452             tarr, Tarr, initial_wait_time, textVerbose)
 453         return self.calcPropBands_ZN_compute_terms(td, spoi)
 454     def stripT(self, on) :
 455         if on :
 456             self.stripT = True
 457             pipename = 'Temp_pipe'
 458             self.stripTpipe = os.popen(pipename)
 459             os.system("stripchart -u %s" % pipename)
 460     def isStable(self, setpoint, tolerance=0.3, maxTime=10.0) :
 461         """
 462         Counts how long the temperature stays within
 463         TOLERANCE of SETPOINT.
 464         Returns when temp goes bad, or MAXTIME elapses.
 465         """
 466         stable = False
 467         startTime = time.time()
 468         stopTime = startTime
 469         while abs(self.getTemp() - setpoint) < tolerance :
 470                 stopTime = time.time()
 471                 if (stopTime-startTime) > maxTime :
 472                         print "Stable for long enough"
 473                         break
 474         if stopTime-startTime > maxTime :
 475                 return True
 476         else :
 477                 return False
 478     def setFilterTime(self, seconds) :
 479         """
 480         Positive values to affect only monitored values.
 481         Negative values affect both monitored and control values.
 482         """
 483         decSeconds = int(seconds*10)
 484         if decSeconds < 0 : # convert (unsigned int) -> (2's compliment signed)
 485                 decSeconds += 2**16
 486         self.write(Melcor.REG_INPUT_SOFTWARE_FILTER_1, decSeconds)
 487     def getFilterTime(self) :
 488         """
 489         Positive values to affect only monitored values.
 490         Negative values affect both monitored and control values.
 491         """
 492         val = self.read(Melcor.REG_INPUT_SOFTWARE_FILTER_1)
 493         if val >= 2**15 : # convert (2's complement signed) -> (unsigned int)
 494                 val -= 2**16
 495         return val/10.0
 496     def sanityCheck(self) :
 497         "Check that some key registers have the values we expect"
 498         self._sanityCheck(Melcor.REG_UNITS_TYPE,   2) # SI
 499         self._sanityCheck(Melcor.REG_C_OR_F,       1) # C
 500         self._sanityCheck(Melcor.REG_FAILURE_MODE, 2) # off
 501         self._sanityCheck(Melcor.REG_RAMPING_MODE, 0) # off
 502         self._sanityCheck(Melcor.REG_OUTPUT_1,     1) # cool
 503         self._sanityCheck(Melcor.REG_OUTPUT_2,     1) # heat
 504     def _sanityCheck(self, register, expected_value) :
 505         val = self.read(register)
 506         if val != expected_value :
 507                 print "Register %d, expected %d, was %d" % (register,
 508                                                     expected_value,
 509                                                     val)
 510                 raise error, "Controller settings error"
 511     def read(self, register) :
 512         """
 513         (register) -> (value)
 514         Returns the value of the specified memory register on the controller.
 515         Registers are defined in the Melcor module.
 516         See melcor_registers.h for a pointers on meanings and manual page nums.
 517         """
 518         (err, val) = self.T.read(register)
 519         if err != 0 :
 520                 raise errorMelcor
 521         return val
 522     def write(self, register, value) :
 523         """
 524         (register, value) -> None
 525         Sets the value of the specified memory register on the controller.
 526         Registers are defined in the Melcor module.
 527         See melcor_registers.h for a pointers on meanings and manual page nums.
 528         """
 529         err = self.T.write(register, value)
 530         if err != 0 :
 531                 raise errorMelcor
 532     def getDeadtimeData(self, num_oscillations=10, curHysteresis=0.8) :
 533         orig_heat_gains = self.getHeatingGains()
 534         orig_cool_gains = self.getCoolingGains()
 535         if self.verbose :
 536             print "Measuring dead time"
 537             print " go to bang-bang"
 538         self.setHeatingGains(0.1, 0, 0)
 539         self.setCoolingGains(0.1, 0, 0)
 540         def isHeating(cur) :
 541             if cur > curHysteresis :
 542                 return True
 543             elif cur < -curHysteresis :
 544                 return False
 545             else :
 546                 return None
 547         i=0
 548         timeArr = [0.0]
 549         temp = self.getTemp()
 550         cur = self.getCurrent()
 551         heat_first = isHeating(cur)
 552         start_time = time.time()
 553         tm = 0
 554         if verbose :
 555             print " Wait to exit hysteresis region"
 556         while heat_first == None and tm < 30:
 557                 temp = t.getTemp()
 558                 cur = t.getCurrent()
 559                 heat_first = isHeating(temp, cur)
 560                 tm = time.time()-start_time
 561         if tm > 30 :
 562                 raise error, "after 30 seconds, still inside hysteresis region"
 563         if self.verbose :
 564             print " Read oscillations"
 565         heating = heat_first
 566         start_time = time.time()
 567         tempArr = [temp]
 568         curArr = [cur]
 569         if verbose :
 570                 print "Temp %g\t(%g),\tCur %g,\tTime %d" % (temp, temp-Tset, cur, 0)
 571         while i < numOscillations*2 :
 572                 temp = t.getTemp()
 573                 tm = time.time()-start_time
 574                 cur = t.getCurrent()
 575                 tempArr.append(temp)
 576                 timeArr.append(tm)
 577                 curArr.append(cur)
 578                 check_signs(temp,cur)
 579                 if heating == True and isHeating(temp, cur) == False :
 580                         print "Transition to cooling (i=%d)" % i
 581                         heating = False
 582                         i += 1
 583                 elif heating == False and isHeating(temp, cur) == True :
 584                         print "Transition to heating (i=%d)" % i
 585                         heating = True
 586                         i += 1
 587         if verbose :
 588             print " Restoring gains"
 589         self.setHeatingGains(*orig_heat_gains)
 590         self.setCoolingGains(*orig_cool_gains)
 591     def time_getTemp(self) :
 592         "Rough estimate timeing of getTemp(), takes me about 0.1s"
 593         start = time.time()
 594         for i in range(10) :
 595             self.getTemp()
 596         stop = time.time()
 597         return (stop-start)/10.0
 598
 599 def _test_tempController() :
 600     t = tempController(controller=1, maxCurrent=0.1)
 601
 602     print "Temp     = %g" % t.getTemp()
 603     print "Current  = %g" % t.getCurrent()
 604     print "Setpoint = %g" % t.getSetpoint()
 605
 606     print "Setting setpoint to 5.0 deg C"
 607     t.setSetpoint(5.0)
 608     sp = t.getSetpoint()
 609     print "Setpoint = %g" % sp
 610     if sp != 5.0 :
 611         raise Exception, "Setpoint in %g != setpoint out %g" % (sp, 5.0)
 612     time.sleep(10) # give the controller some time to overcome any integral gain
 613     c = t.getCurrent()
 614     print "Current  = %g" % c
 615     mca, mcb, mct = t.getMaxCurrent()
 616     if t.getTemp() < mct : # we're below the high power limit setpoint, use mcb
 617         if c != mcb :
 618             raise Exception, "Current not at max %g, and we're shooting for a big temp" % mcb
 619     else :
 620         if c != mca :
 621             raise Exception, "Current not at max %g, and we're shooting for a big temp" % mca
 622
 623
 624     print "Setting setpoint to 50.0 deg C"
 625     t.setSetpoint(50.0)
 626     sp = t.getSetpoint()
 627     print "Setpoint = %g" % sp
 628     if sp != 5.0 :
 629         raise Exception, "Setpoint in %g != setpoint out %g" % (sp, 5.0)
 630     time.sleep(10)
 631     c = t.getCurrent()
 632     print "Current  = %g" % c
 633     print "Success"
 634     mca, mcb, mct = t.getMaxCurrent()
 635     if t.getTemp() < mct : # we're below the high power limit setpoint, use mcb
 636         if -c != mcb :
 637             raise Exception, "Current not at min %g, and we're shooting for a big temp" % (-mcb)
 638     else :
 639         if -c != mca :
 640             raise Exception, "Current not at min %g, and we're shooting for a big temp" % (-mca)
 641
 642 def test() :
 643     _test_tempController()
 644
 645 if __name__ == "__main__" :
 646     test()