diff --git a/controls/__init__.py b/controls/__init__.py new file mode 100644 index 0000000..b7ebec2 --- /dev/null +++ b/controls/__init__.py @@ -0,0 +1,18 @@ +# -*- coding: utf-8 -*- +'''Chemical Engineering Design Library (ChEDL). Utilities for process modeling. +Copyright (C) 2016, Caleb Bell + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program. If not, see .''' + +from __future__ import division diff --git a/controls/controller_tuning.py b/controls/controller_tuning.py new file mode 100644 index 0000000..1198d56 --- /dev/null +++ b/controls/controller_tuning.py @@ -0,0 +1,391 @@ +# -*- coding: utf-8 -*- +'''Chemical Engineering Design Library (ChEDL). Utilities for process modeling. +Copyright (C) 2016, Caleb Bell + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program. If not, see .''' + +from __future__ import division + + + +def IAE(ys, ysp): + r'''Calculates integral absolute error (IAE) for a variable. Also + called the Sum of Absolute Errors (SAE) . + + .. math:: + IAE = \int_0^\infty |y_{sp}(t) - y_s(t)|dt + + Parameters + ---------- + ys : array-like + Varying values of y + ysp : float + Desired value for y ("set point") + + Returns + ------- + IAE : float + Integral squared error + + Examples + -------- + >>> IAE([9.9, 9.92, 9.94, 9.96], 10) + 0.27999999999999936 + ''' + iae = sum([abs(yi-ysp) for yi in ys]) + return iae + + +def ISE(ys, ysp): + r'''Calculates integral squared error (ISE) for a variable. Also + called the Sum of Squared Errors (SSE) . + + .. math:: + ISE = \int_0^\infty |y_{sp}(t) - y_s(t)|^2dt + + Parameters + ---------- + ys : array-like + Varying values of y + ysp : float + Desired value for y ("set point") + + Returns + ------- + IAE : float + Integral squared error + + Examples + -------- + >>> IAE([9.9, 9.92, 9.94, 9.96], 10) + 0.021599999999999935 + ''' + ise = sum([abs(yi-ysp)**2 for yi in ys]) + return ise + + +def ITAE(ys, ts, ysp): + r'''Calculates Integral Time Absolute Error (ITAE) for a variable. + + .. math:: + ITAE = \int_0^\infty t|y_{sp}(t) - y_s(t)|dt + + Parameters + ---------- + ys : array-like + Varying values of y + ts : array-like + Varying values of time, [time] + ysp : float + Desired value for y ("set point") + + Returns + ------- + ITAE : float + Integral absolute error + + Examples + -------- + >>> ITAE([9.9, 9.92, 9.94, 9.96], [0, 1, 2, 3], 10) + 0.3199999999999985 + ''' + itae = sum([ts[i]*abs(ys[i]-ysp) for i in range(len(ys))]) + return itae + + +def ITSE(ys, ts, ysp): + r'''Calculates Integral Time Squared Error (ITSE) for a variable. + + .. math:: + ITSE = \int_0^\infty t|y_{sp}(t) - y_s(t)|^2 dt + + Parameters + ---------- + ys : array-like + Varying values of y + ts : array-like + Varying values of time, [time] + ysp : float + Desired value for y ("set point") + + Returns + ------- + ITSE : float + Integral squared error + + Examples + -------- + >>> ITAE([9.9, 9.92, 9.94, 9.96], [0, 1, 2, 3], 10) + 0.018399999999999927 + ''' + itse = sum([ts[i]*abs(ys[i]-ysp)**2 for i in range(len(ys))]) + return itse + +### Controller form conversions + +def interactive_to_noninteractive(Kc, tauI, tauD): + r'''Converts `Kc`, `tauI`, and `tauD` from the interactive (older) + controller form to the ideal or standard form (newer). Primes designate + interactive terms. + + .. math:: + K_c = K_c'\left(\frac{\tau_I' + \tau_D'}{\tau_I'}\right) + + \tau_I = \tau_I' + \tau_D' + + \tau_D = \frac{\tau_I' \tau_D'}{\tau_i' + \tau_D'} + + Parameters + ---------- + Kc : float + Controller gain, [] + tauI : float + Integral time constant, [time] + tauD : float + Derivative time constant, [time] + + Returns + ------- + Kc : float + Controller gain, [] + tauI : float + Integral time constant, [time] + tauD : float + Derivative time constant, [time] + + Examples + -------- + >>> interactive_to_noninteractive(0.07236, 7.236, 2.764) + (0.09999999999999999, 10.0, 2.0000303999999995) + ''' + Kc_non = Kc*(tauI + tauD)/tauI + tauI_non = tauI + tauD + tauD_non = (tauI*tauD)/(tauI + tauD) + return Kc_non, tauI_non, tauD_non + +#print interactive_to_noninteractive(0.07236, 7.236, 2.764) + +def noninteractive_to_interactive(Kc, tauI, tauD): + r'''Converts `Kc`, `tauI`, and `tauD` from the ideal or standard controller + form (newer) to the interactive (older) form. Primes designate interactive + terms. + + .. math:: + K_c' = K_c\left(0.5 + \sqrt{0.25 - \tau_D/\tau_I}\right) + + \tau_I' = \tau_I\left(0.5 + \sqrt{0.25 - \tau_D/\tau_I}\right) + + \tau_D' = \frac{\tau_D}{0.5 + \sqrt{0.25 - \tau_D/\tau_I}} + + Parameters + ---------- + Kc : float + Controller gain, [] + tauI : float + Integral time constant, [time] + tauD : float + Derivative time constant, [time] + + Returns + ------- + Kc : float + Controller gain, [] + tauI : float + Integral time constant, [time] + tauD : float + Derivative time constant, [time] + + Notes + ----- + If tauD/tauI < 0.25, the conversion is not possible - the square root of a + negative number is being taken. + + Examples + -------- + >>> noninteractive_to_interactive(0.1, 10., 2.) + (0.0723606797749979, 7.23606797749979, 2.7639320225002106) + ''' + if tauD/tauI > 0.25: + raise Exception(ValueError) + Kc_int = Kc*(0.5 + (0.25 - tauD/tauI)**0.5) + tauI_int = tauI*(0.5 + (0.25 - tauD/tauI)**0.5) + tauD_int = tauD/(0.5 + (0.25 - tauD/tauI)**0.5) + return Kc_int, tauI_int, tauD_int + + +def noninteractive_to_parallel(Kc, tauI, tauD): + r'''Converts `Kc`, `tauI`, and `tauD` from the ideal or standard + controller form to `Kp`, `KI`, and `KD` in the parallel form. + + .. math:: + K_p = K_c + + K_I = \frac{K_c}{\tau_I} + + K_D = K_c\tau_D + + Parameters + ---------- + Kc : float + Controller gain, [] + tauI : float + Integral time constant, [time] + tauD : float + Derivative time constant, [time] + + Returns + ------- + Kp : float + Proportional gain, [] + KI : float + Integral time constant, [1/time] + KD : float + Derivative time constant, [time] + + Examples + -------- + >>> noninteractive_to_parallel(0.1, 10., 2.) + (0.1, 0.01, 0.2) + ''' + Kp = Kc + KI = Kc/tauI + KD = Kc*tauD + return Kp, KI, KD + +#print noninteractive_to_parallel(0.1, 10., 2.) + +def parallel_to_noninteractive(Kp, KI, KD): + r'''Converts `Kp`, `KI`, and `KD` in the parallel controller form to + `Kc`, `tauI`, and `tauD` in the ideal or standard controller form. + + .. math:: + K_c = K_p + + \tau_I = K_c/K_I + + \tau_D = K_D/K_c + + Parameters + ---------- + Kp : float + Proportional gain, [] + KI : float + Integral time constant, [1/time] + KD : float + Derivative time constant, [time] + + Returns + ------- + Kc : float + Controller gain, [] + tauI : float + Integral time constant, [time] + tauD : float + Derivative time constant, [time] + + Examples + -------- + >>> parallel_to_noninteractive(0.1, 0.01, 0.2) + (0.1, 10.0, 2.0) + ''' + Kc = Kp + tauI = Kc/KI + tauD = KD/Kc + return Kc, tauI, tauD + +#print parallel_to_noninteractive(0.1, 0.01, 0.2) + +def parallel_to_interactive(Kp, KI, KD): + r'''Converts Kp`, `KI`, and `KD` in the parallel controller form to `Kc`, + `tauI`, and `tauD` from the interactive (older) controller form. + + Parameters + ---------- + Kp : float + Proportional gain, [] + KI : float + Integral time constant, [1/time] + KD : float + Derivative time constant, [time] + + Returns + ------- + Kc : float + Controller gain, [] + tauI : float + Integral time constant, [time] + tauD : float + Derivative time constant, [time] + + Notes + ----- + Utilizes the functions parallel_to_noninteractive and + noninteractive_to_interactive to perform the conversions. + If tauD/tauI < 0.25, in the intermediary step, the conversion is not + possible - the square root of a negative number is being taken. + + Examples + -------- + >>> parallel_to_interactive(0.1, 0.01, 0.2) + (0.0723606797749979, 7.23606797749979, 2.7639320225002106) + ''' + Kc, tauI, tauD = parallel_to_noninteractive(Kp, KI, KD) + Kc, tauI, tauD = noninteractive_to_interactive(Kc, tauI, tauD) + return Kc, tauI, tauD + +#print parallel_to_interactive(0.1, 0.01, 0.2) + + +def interactive_to_parallel(Kc, tauI, tauD): + r'''Converts `Kc`, `tauI`, and `tauD` from the interactive (older) + controller form to `Kp`, `KI`, and `KD` in the parallel form. + + Parameters + ---------- + Kc : float + Controller gain, [] + tauI : float + Integral time constant, [time] + tauD : float + Derivative time constant, [time] + + Returns + ------- + Kp : float + Proportional gain, [] + KI : float + Integral time constant, [1/time] + KD : float + Derivative time constant, [time] + + Notes + ----- + Utilizes the functions interactive_to_noninteractive and + noninteractive_to_parallel to perform the conversions. + If tauD/tauI < 0.25, in the intermediary step, the conversion is not + possible - the square root of a negative number is being taken. + + Examples + -------- + >>> interactive_to_parallel(0.07236, 7.236, 2.764) + (0.09999999999999999, 0.009999999999999998, 0.20000303999999994) + ''' + Kc, tauI, tauD = interactive_to_noninteractive(Kc, tauI, tauD) + Kp, KI, KD = noninteractive_to_parallel(Kc, tauI, tauD) + return Kp, KI, KD + +#print interactive_to_parallel(0.07236, 7.236, 2.764) + + +### Controller form time unit conversions \ No newline at end of file diff --git a/controls/rtd.py b/controls/rtd.py new file mode 100644 index 0000000..46ebac3 --- /dev/null +++ b/controls/rtd.py @@ -0,0 +1,269 @@ +# -*- coding: utf-8 -*- +'''Chemical Engineering Design Library (ChEDL). Utilities for process modeling. +Copyright (C) 2016, Caleb Bell + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program. If not, see .''' + +from __future__ import division +from math import log, exp, sqrt +from numpy import roots, imag, real +from controls.sensor import Sensor + +# Coefs checked with original, from a NIST document +coefs_low = [-2.13534729, 3.18324720, -1.80143597, 0.71727204, 0.50344027, -0.61893395, -0.05332322, 0.28021362, 0.10715224, -0.29302865, 0.04459872, 0.11868632, -0.05248134] +coefs_high = [2.78157254, 1.64650916, -0.13714390, -0.00649767, -0.00234444, 0.00511868, 0.00187982, -0.00204472, -0.00046122, 0.00045724] +invcoefs_low = [0.183324722, 0.240975303, 0.209108771, 0.190439972, 0.142648498, 0.077993465, 0.012475611, -0.032267127, -0.075291522, -0.056470670, 0.076201285, 0.123893204, -0.029201193, -0.091173542, 0.001317696, 0.026025526] +invcoefs_high = [439.932854, 472.418020, 37.684494, 7.472018, 2.920828, 0.005184, -0.963864, -0.188732, 0.191203, 0.049025] + +def RTD_resistance_ITS90(T): + r'''Calculates the resisistance of a Resistance Temperature Detector, + relative to its value at 273.15 K. Two ranges are used, one below 273.16 K + until 13.8033 K and one above it, up to 961.78 degrees Celcius. + + .. math:: + \ln W_r(T_{90}) = A_0 + \sum_{i=1}^{12} A_i + \left[\frac{\ln(T_{90}/273.16K) + 1.5}{1.5}\right]^i + + W_r(T_{90}) = C_0 + \sum_{i=1}^9 C_i\left[\frac{T_{90}/K + - 754.15}{481}\right]^i + + Parameters + ---------- + T : float + Temperature, [K] + + Returns + ------- + w : float + Resistance of RTD relative to its value at 273.15 K, [-] + + Notes + ----- + No warnings if outside the suggested range. + + Examples + -------- + >>> RTD_resistance_ITS90(150.), RTD_resistance_ITS90(350.) + (0.4984000571593147, 1.302905074678503) + + References + ---------- + .. [1] Preston-Thomas, H. "The International Temperature Scale of 1990 + (ITS-90)." Metrologia 27, no. 1 (1990): 3. + doi:10.1088/0026-1394/27/1/002. + ''' + if T <= 273.16: + tot = [coefs_low[i]*((log(T/273.16) + 1.5)/1.5)**i for i in range(1,len(coefs_low))] + w = exp(coefs_low[0] + sum(tot)) + else: + tot = [coefs_high[i]*((T-754.15)/481.)**i for i in range(1,len(coefs_high))] + w = coefs_high[0] + sum(tot) + return w + + +def RTD_temperature_ITS90(w): + r'''Calculates the temperature of a Resistance Temperature Detector, + given its resistance relative to its value at 273.15 K. Two ranges are + used, one below 273.16 K until 13.8033 K and one above it, up to 961.78 + degrees Celcius. + + .. math:: + T_{90}/273.16K = B_0 + \sum_{i=1}^{15} B_i + \left[\frac{W_r(T_{90})^{1/6} - 0.65}{0.35}\right]^i + + T_{90}/K - 273.15 = D_0 + \sum_{i=1}^9 D_i\left[\frac{W_r(T_{90})- + 2.64}{1.64}\right]^i + + Parameters + ---------- + w : float + Resistance of RTD relative to its value at 273.15 K, [-] + + Returns + ------- + T : float + Temperature, [K] + + Notes + ----- + No warnings if outside the suggested range. + The first funcntino is equivalent to the normal function to within 0.1 mK, + and the second to 0.013 mK. + + Examples + -------- + >>> RTD_temperature_ITS90(0.4984000571593147), RTD_temperature_ITS90(1.302905074678503) + (150.00061624978588, 350.0000249661761) + + References + ---------- + .. [1] Preston-Thomas, H. "The International Temperature Scale of 1990 + (ITS-90)." Metrologia 27, no. 1 (1990): 3. + doi:10.1088/0026-1394/27/1/002. + ''' + if w < 1.0000599918: # low + tot = [invcoefs_low[i]*((w**(1/6.) -0.65)/0.35)**i for i in range(1,len(invcoefs_low))] + T = (invcoefs_low[0] + sum(tot))*273.16 + else: + tot = [invcoefs_high[i]*((w - 2.64)/1.64)**i for i in range(1,len(invcoefs_high))] + T = (invcoefs_high[0] + sum(tot))+273.15 + return T + + + +def Callendar_Van_Dusen_resistance(T, A=3.9083E-3, B=-5.775E-7, C=-4.183E-12): + r'''Calculates the resisistance of a Resistance Temperature Detector, + relative to its value at 273.15 K. Two ranges are used, one below 273.16 K + until 73.15 K and one above it, up to 850 degrees Celcius. + + .. math:: + W_r = 1 + At + Bt^2 + Ct^3(t-100) + + W_r = 1 + At + Bt^2 + + t = T - 273.15 + + Parameters + ---------- + T : float + Temperature, [K] + A, B, C : floats, optional + Coefficients for the Callendar Van Dusen equation, [various] + + Returns + ------- + w : float + Resistance of RTD relative to its value at 273.15 K, [-] + + Notes + ----- + No warnings if outside the suggested range. + + A = 3.9083E-3 + B = -5.775E-7, + C = -4.183E-12 + + Examples + -------- + >>> Callendar_Van_Dusen_resistance(150.), Callendar_Van_Dusen_resistance(350.) + (0.508191171034818, 1.2969421847562501) + + References + ---------- + .. [1] IEC 60751 Standard, as given in Omega RTD Specifications, Technical + Reference at http://www.omega.com/Temperature/pdf/RTDSpecs_Ref.pdf + ''' + t = T - 273.15 + if t < 0: + w = 1 + A*t + B*t**2 + C*t**3*(t-100) + else: + w = (1 + A*t + B*t**2) + return w + +#print RTD_resistance_ITS90(273.16), RTD_resistance_ITS90(200.) +#print [Callendar_Van_Dusen_temperature(150.), Callendar_Van_Dusen_temperature(350.)], 'hi' + +def Callendar_Van_Dusen_temperature(w, A=3.9083E-3, B=-5.775E-7, C=-4.183E-12): + r'''Calculates the temperature of a Resistance Temperature Detector, + from its resistance relative to its value at 273.15 K. Two ranges are used, + one below 273.16 K until 73.15 K and one above it, up to 850 degrees + Celcius. The following equations are solved by a root-finding algorithm, + or the quadratic equation is necessary. + + .. math:: + W_r = 1 + At + Bt^2 + Ct^3(t-100) + + W_r = 1 + At + Bt^2 + + t = T - 273.15 + + Parameters + ---------- + w : float + Resistance of RTD relative to its value at 273.15 K, [-] + A, B, C : floats, optional + Coefficients for the Callendar Van Dusen equation, [various] + + Returns + ------- + T : float + Temperature, [K] + + Notes + ----- + No warnings if outside the suggested range. + + A = 3.9083E-3 + B = -5.775E-7, + C = -4.183E-12 + + Examples + -------- + >>> Callendar_Van_Dusen_temperature(0.508191171034818), Callendar_Van_Dusen_temperature(1.2969421847562501) + (150.00000000000006, 349.9999999999998) + + References + ---------- + .. [1] IEC 60751 Standard, as given in Omega RTD Specifications, Technical + Reference at http://www.omega.com/Temperature/pdf/RTDSpecs_Ref.pdf + ''' + # Try the case of > 273.15 K + # IEC 60751 standard says neglect A for positive temperatures + t = (-A + (A**2 - 4*B*(1-w))**0.5)/(2*B) + if t > 0: + T = t +273.15 + else: + ts = roots([C, -100*C, B, A, 1-w]) + T = [real(i)+273.15 for i in ts if not imag(i)][-1] + return T +# + + +class RTD(Sensor): + '''RTD Class''' + def __init__(self, + variable_units='K', signal_units='ohm', + deadtime=0, offset=0, R0=100., method='IEC 60751'): + self.R0 = R0 + self.method = method + + if self.method == 'ITS-90': + variable_min = 13.81 + variable_max = 630.75 + 273.15 + signal_min = RTD_resistance_ITS90(variable_min) + signal_max = RTD_resistance_ITS90(variable_max) + self.signal_to_variable_function = RTD_temperature_ITS90 + self.variable_to_signal_function = RTD_resistance_ITS90 + elif 'IEC 60751' in self.method: + variable_min = 73.15 + variable_max = 850 + 273.15 + signal_min = Callendar_Van_Dusen_resistance(variable_min) + signal_max = Callendar_Van_Dusen_resistance(variable_max) + self.signal_to_variable_function = Callendar_Van_Dusen_temperature + self.variable_to_signal_function = Callendar_Van_Dusen_resistance + else: + raise Exception(ValueError) + + Sensor.__init__(self, variable_min=variable_min, variable_max=variable_max, + variable_units=variable_units, signal_min=signal_min*R0, + signal_max=signal_max*R0, signal_units=signal_units, + deadtime=deadtime, offset=offset) + + + def signal_to_variable(self, signal): + # Signal as sensed by the R0 resistor + + T = self.signal_to_variable_function(signal/self.R0) + return T + diff --git a/controls/sensor.py b/controls/sensor.py new file mode 100644 index 0000000..100e48c --- /dev/null +++ b/controls/sensor.py @@ -0,0 +1,154 @@ +# -*- coding: utf-8 -*- +'''Chemical Engineering Design Library (ChEDL). Utilities for process modeling. +Copyright (C) 2016, Caleb Bell + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program. If not, see .''' + +from __future__ import division +import ctypes +import scipy.stats as stats + + + + + + +__random_distributions = ['anglit', 'arcsine', 'cauchy', 'cosine', 'expon', + 'gilbrat', 'halfcauchy', 'halflogistic', 'halfnorm', + 'hypsecant', 'kstwobign', 'laplace', 'levy', + 'logistic', 'maxwell', 'norm', 'rayleigh', + 'semicircular', 'uniform', 'wald'] + + +class Sensor(object): + ''' + Sensor package includes transmitter, variable reader, and settings. + Can `output` the variable directly, or 4-20 mA or a-b V signal. + `sensed variable` is the desired information, i.e. T for a thermocouple. + `sensed signal + + The `output` max and min are different than the `variable` min and max; + functional relationships normally allow the `variable` to be calculated + in a much larger range than the `output` can be set. If the `output` + range is not set, its range is set to those of the variable. + ''' + t = 0 # Initialize at time = 0 + linear = False # Assume non-linear signal-variable relation + + def __init__(self, output_max=None, output_min=None, output_type='variable', + output_units=None, output_action='direct', + output_rounding_algorithm=None, + variable_min=None, variable_max=None, variable_units=None, + signal_min=None, signal_max=None, signal_units=None, + signal_rounding_algorithm=None, signal_error_distribution='norm', + deadtime=None, inaccuracy=None, offset=None): + + + self.variable_min = self.variable_zero = variable_min + self.variable_max = variable_max + self.variable_span = self.variable_max - self.variable_min + self.variable_units = variable_units + + self.signal_min = self.signal_zero = signal_min + self.signal_max = signal_max + self.signal_span = self.signal_max - self.signal_min + self.signal_units = signal_units + self.signal_rounding_algorithm = signal_rounding_algorithm + self.signal_error_distribution = signal_error_distribution + + self.output_max = output_max + self.output_min = output_min + self.output_type = output_type + self.output_units = output_units + self.output_action = output_action + self.output_rounding_algorithm = output_rounding_algorithm + self.__initialize_outputs__() + # Must add output units here too, based on output type + + self.deadtime = deadtime + self.inaccuracy = inaccuracy + self.offset = offset + + def __initialize_outputs__(self): + if self.output_type == 'variable': + self.output_min = self.output_zero = self.variable_min + self.output_max = self.variable_max + self.output_units = self.variable_units + elif self.output_type == '4-20 mA': + self.output_min = self.output_zero = 4. + self.output_max = 20. + self.output_span = self.output_max - self.output_min + self.output_units = 'mA' + + elif self.output_type == '1-5 V': + self.output_min = self.output_zero = 1. + self.output_max = 5. + self.output_span = self.output_max - self.output_min + self.output_units = 'V' + elif self.output_type == '0-5 V': + self.output_min = self.output_zero = 0. + self.output_max = 5. + self.output_span = self.output_max - self.output_min + self.output_units = 'V' + else: + raise ValueError(self.output_type) + + def output(self, signal): + self.variable = self.signal_to_variable(signal) + self.output = self.variable_to_output(self.variable) + return self.output + + # Dummy function, to be overwridden by those in subclasses. + return signal + + def variable_to_output(self, variable): + if self.output_type == 'variable': + output = variable + elif self.output_action == 'direct': + output = ((variable-self.variable_min) + /(self.variable_max - self.variable_min) + *(self.output_max - self.output_min)) + self.output_min + elif self.output_action == 'reverse': + output = -((variable-self.variable_min) + /(self.variable_max - self.variable_min) + *(self.output_max - self.output_min)) + self.output_max + else: + raise Exception(ValueError) + if self.output_rounding_algorithm: + output = self.rounding(output, output_rounding_algorithm, self.output_max, self.output_min) + + return output + + def rounding(self, unrounded, rounding_algorithm=None, rounded_max=None, rounded_min=None): + if rounding_algorithm == 'int': + output = int(unrounded) + elif rounding_algorithm == 'int8': + output = ctypes.c_int8(unrounded).value + elif rounding_algorithm == 'int16': + output = ctypes.c_int16(unrounded).value + elif rounding_algorithm == 'int32': + output = ctypes.c_int32(unrounded).value + elif rounding_algorithm == 'float': + output = ctypes.c_float(unrounded).value + elif 'bit' in rounding_algorithm: + bits = int(rounding_algorithm.replace('bit', '')) + output = unrounded - unrounded % (rounded_max - rounded_min)/2.**bits + else: + raise Exception(ValueError) + return output + + def signal_error(self, signal): + if self.inaccuracy and self.signal_error_distribution: + signal += signal*getattr(stats, dist).rvs(loc=inaccuracy) + return signal diff --git a/controls/thermistor.py b/controls/thermistor.py new file mode 100644 index 0000000..6aad9e2 --- /dev/null +++ b/controls/thermistor.py @@ -0,0 +1,198 @@ +# -*- coding: utf-8 -*- +'''Chemical Engineering Design Library (ChEDL). Utilities for process modeling. +Copyright (C) 2016, Caleb Bell + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program. If not, see .''' + +from __future__ import division +from math import log, exp, sqrt +# Thermistors are calibrated for specific resistances ONLY! + +def Steinhart_Hart_temperature(R, A, B, C): + r'''Calculates the temperature of a Thermistor, using its absolute + resistance and Steinhard-Hart coefficients according to [1]_. Coefficients + are calibrated for every single thermistor separately, due to large + variations in production methods. + + .. math:: + \frac{1}{T} = A + B(\ln R) + C(\ln R)^3 + + Parameters + ---------- + R : float + Resistance of thermistor at T [ohm] + A, B, C : floats + Coefficients for Steinhart-Hart function, [-] + + Returns + ------- + T : float + Temperature, [K] + + Notes + ----- + No warnings produced if outside the suggested range. + Thermistors are normally valid only for short ranges under 200 K. + + Examples + -------- + Coefficients for Omega 44004 Thermistor, 2252 ohm at STP: + + >>> Steinhart_Hart_temperature(2252, A=1.468E-3, B=2.383E-4, C=1.007E-7) + 298.1604531883543 + + References + ---------- + .. [1] Steinhart, John S., and Stanley R. Hart. "Calibration Curves for + Thermistors." Deep Sea Research and Oceanographic Abstracts 15, no. 4 + (August 1968): 497-503. doi:10.1016/0011-7471(68)90057-0. + ''' + T = (A + B*log(R) + C*(log(R))**3)**-1 + return T + + +def Steinhart_Hart_resistance(T, A, B, C): + r'''Calculates the resistance of a Thermistor, using its absolute + temperature and Steinhard-Hart coefficients according to [1]_. Coefficients + are calibrated for every single thermistor separately, due to large + variations in production methods. + + .. math:: + \frac{1}{T} = A + B(\ln R) + C(\ln R)^3 + + Parameters + ---------- + T : float + Temperature, [K] + A, B, C : floats + Coefficients for Steinhart-Hart function, [-] + + Returns + ------- + R : float + Resistance of thermistor at T [ohm] + + Notes + ----- + No warnings produced if outside the suggested range. + Thermistors are normally valid only for short ranges under 200 K. + This equation was produced by Sympy, and matches reverse calculations. + Examples + -------- + Coefficients for Omega 44004 Thermistor, 2252 ohm at STP: + + >>> Steinhart_Hart_resistance(2252, A=1.468E-3, B=2.383E-4, C=1.007E-7) + 2253.033420231679 + + References + ---------- + .. [1] Steinhart, John S., and Stanley R. Hart. "Calibration Curves for + Thermistors." Deep Sea Research and Oceanographic Abstracts 15, no. 4 + (August 1968): 497-503. doi:10.1016/0011-7471(68)90057-0. + ''' + R = exp(B/(3*C*(sqrt(B**3/(27*C**3) + (A*T - 1)**2/(4*C**2*T**2)) + + (A*T - 1)/(2*C*T))**(1/3)) - (sqrt(B**3/(27*C**3) + + (A*T - 1)**2/(4*C**2*T**2)) + (A*T - 1)/(2*C*T))**(1/3)) + return R + + +def beta_temperature(R, R0, B=3750., T0=298.15): + r'''Calculates the temperature of a thermistor, using its absolute + resistance and beta coefficient according to [1]_. Coefficients + are calibrated for every single thermistor separately, due to large + variations in production methods. Default value is provided as an example + only. Also necessary is a reference temperature, and resistivity at the + reference temperature. + + .. math:: + T = \frac{B}{\log[\frac{R_T}{R_{T0}}\exp(B/T_0)]} + + Parameters + ---------- + R : float + Resistance of thermistor at T [ohm] + R0 : float + Resistance of thermistor at T0 [ohm] + B : float, optional + Beta coefficient, [K] + T0 : float, optional + Reference temperature, [K] + + Returns + ------- + T : float + Temperature, [K] + + Notes + ----- + No warnings produced if outside the suggested range. + Thermistors are normally valid only for short ranges under 200 K. + + Examples + -------- + >>> beta_temperature(2.5E4, R0=1E4, B=3799.42) + 278.1500055180414 + + References + ---------- + .. [1] Liptak, Bela G. Instrument Engineers' Handbook, 4E, Volume One: + Process Measurement and Analysis. CRC Press, 2003. + ''' + T = B*(log(R/R0*exp(B/T0)))**-1 + return T + +def beta_resistance(T, R0, B=3750., T0=298.15): + r'''Calculates the resistance of a thermistor, using its absolute + temperature and beta coefficient according to [1]_. Coefficients + are calibrated for every single thermistor separately, due to large + variations in production methods. Default value is provided as an example + only. Also necessary is a reference temperature, and resistivity at the + reference temperature. + + .. math:: + R_T =R_{TO} \exp\left(B\left(\frac{1}{T} - \frac{1}{T_0}\right)\right) + + Parameters + ---------- + T : float + Temperature, [K] + R0 : float + Resistance of thermistor at T0 [ohm] + B : float, optional + Beta coefficient, [K] + T0 : float, optional + Reference temperature, [K] + + Returns + ------- + R : float + Resistance of thermistor at T [ohm] + + Notes + ----- + No warnings produced if outside the suggested range. + Thermistors are normally valid only for short ranges under 200 K. + + Examples + -------- + >>> beta_resistance(5+273.15, R0=1E4, B=3799.42, T0=298.15) + 25000.00677460831 + + References + ---------- + .. [1] Liptak, Bela G. Instrument Engineers' Handbook, 4E, Volume One: + Process Measurement and Analysis. CRC Press, 2003. + ''' + R = R0*exp(B*(1./T - 1./T0)) + return R diff --git a/controls/thermocouple.py b/controls/thermocouple.py new file mode 100644 index 0000000..4a9b1f6 --- /dev/null +++ b/controls/thermocouple.py @@ -0,0 +1,255 @@ +# -*- coding: utf-8 -*- +'''Chemical Engineering Design Library (ChEDL). Utilities for process modeling. +Copyright (C) 2016, Caleb Bell + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program. If not, see .''' + +from __future__ import division +from math import exp +from controls.sensor import Sensor +from numpy import roots, real, imag +from scipy.constants import * + + +def calc_EMF(T, coefs_splits, coefs_piecewise, coefs_exponential): + t = T - 273.15 + for i in range(len(coefs_splits)-1): + if T < coefs_splits[i+1]: + break + + E = sum([c*t**j for j, c in enumerate(coefs_piecewise[i])]) + if coefs_exponential and T > 273.15: + # Defined for K thermocouple only, otherwise a generic definition + # would be needed + E += coefs_exponential[0]*exp(coefs_exponential[1]*((t-coefs_exponential[2])**2)) + return E + +def calc_T(E, Tref, coefs_splits, coefs_piecewise, coefs_exponential): + Eref = calc_EMF(Tref, coefs_splits, coefs_piecewise, coefs_exponential) + # Brute force the polynomials to determine the solution + for i in range(len(coefs_splits)-1): + coefs = coefs_piecewise[i] + Tmin_coefs, Tmax_coefs = coefs_splits[i], coefs_splits[i+1] + + coefs.reverse() + coefs[-1] -= E + Eref + ans = roots(coefs) + for j in ans: + if not imag(j) and real(j) + 273.15 > Tmin_coefs and real(j) + 273.15 < Tmax_coefs: + return real(j) + 273.15 + # If the above didn't work, T is outside of safe limits. + # TODO: Try the lower limit, then the higher limit, with fsolve. + # No guarantees though. + return None + +#def sum_emfs(t, coefs): +# E = sum([c*t**i for i, c in enumerate(coefs)]) +# return E + +#def thermocouple(V=None, T=None, Tref=298.15, kind='K'): +## T = T - 273.15 +# if kind == 'K': +# ais = [0.118597600000E+00, -0.118343200000E-03, 0.126968600000E+03] +# if T < 0.000: +# if T < -270: +# coefs = [0.000000000000E+00, 0.394501280250E-01, +# 0.236223735980E-04, -0.328589067840E-06, +# -0.499048287770E-08, -0.675090591730E-10, +# -0.574103274280E-12, -0.310888728940E-14, +# -0.104516093650E-16, -0.198892668780E-19, +# -0.163226974860E-22] +# E = sum_emfs(T, Tref, coefs) +# else: +# coefs = [-0.176004136860E-01, 0.389212049750E-01, +# 0.185587700320E-04, -0.994575928740E-07, +# 0.318409457190E-09, -0.560728448890E-12, +# 0.560750590590E-15, -0.320207200030E-18, +# 0.971511471520E-22, -0.121047212750E-25] +# E = sum_emfs(T, Tref, coefs) + ais[0]*exp(ais[1]*((T-ais[2])**2)) - ais[0]*exp(ais[1]*((Tref-ais[2])**2)) +# return E + +#print thermocouple(T=70.6, Tref=21.3) + + +Thermocouple_list = ['Au-Pt', 'E', 'G', 'IrRh 0.4 vs Ir', 'J', 'K', 'K vs AuFe 0.07', 'M', 'N', 'P', 'PtMo 0.05 vs PtMo 0.001', 'Pt-Pd', 'PtRh 0.4 vs PtRh 0.2', 'R', 'S', 'T', 'B'] + +class Thermocouple(Sensor): + '''RTD Class''' + def __init__(self, + variable_units='K', signal_units='mV', + deadtime=0, offset=0, output_rounding=None, + method='K', Tref=298.15): + self.method = method + self.Tref = Tref + self.signal_to_variable_function = calc_T + self.variable_to_signal_function = calc_EMF + self.coefs_exponential = None + + if self.method == 'B': + self.coefs_piecewise = [[0.000000000000E+00, -0.246508183460E-03, 0.590404211710E-05, -0.132579316360E-08, 0.156682919010E-11, -0.169445292400E-14, 0.629903470940E-18], + [-0.389381686210E+01, 0.285717474700E-01, -0.848851047850E-04, 0.157852801640E-06, -0.168353448640E-09, 0.111097940130E-12, -0.445154310330E-16, 0.989756408210E-20, -0.937913302890E-24]] + self.coefs_splits = [273.15, 630.615 + 273.15, 1820.000 + 273.15] + self.description = 'Pt-30% Rh versus Pt-6% Rh' + elif self.method == 'E': + self.coefs_piecewise = [[0.000000000000E+00, 0.586655087080E-01, 0.454109771240E-04, -0.779980486860E-06, -0.258001608430E-07, -0.594525830570E-09, -0.932140586670E-11, -0.102876055340E-12, -0.803701236210E-15, -0.439794973910E-17, -0.164147763550E-19, -0.396736195160E-22, -0.558273287210E-25, -0.346578420130E-28], + [0.000000000000E+00, 0.586655087100E-01, 0.450322755820E-04, 0.289084072120E-07, -0.330568966520E-09, 0.650244032700E-12, -0.191974955040E-15, -0.125366004970E-17, 0.214892175690E-20, -0.143880417820E-23, 0.359608994810E-27]] + self.coefs_splits = [3.15, 273.15, 1000 + 273.15] + self.description = 'Ni-Cr alloy versus a Cu-Ni alloy' + elif self.method == 'J': + self.coefs_piecewise = [[0.000000000000E+00, 0.503811878150E-01, 0.304758369300E-04, -0.856810657200E-07, 0.132281952950E-09, -0.170529583370E-12, 0.209480906970E-15, -0.125383953360E-18, 0.156317256970E-22], + [0.296456256810E+03, -0.149761277860E+01, 0.317871039240E-02, -0.318476867010E-05, 0.157208190040E-08, -0.306913690560E-12]] + self.coefs_splits = [63.15, 760 + 273.15, 1200.000 + 273.15] + self.description = 'Fe versus a Cu-Ni alloy' + elif self.method == 'K': + self.coefs_piecewise = [[0.000000000000E+00, 0.394501280250E-01, 0.236223735980E-04, -0.328589067840E-06, -0.499048287770E-08, -0.675090591730E-10, -0.574103274280E-12, -0.310888728940E-14, -0.104516093650E-16, -0.198892668780E-19, -0.163226974860E-22], + [-0.176004136860E-01, 0.389212049750E-01, 0.185587700320E-04, -0.994575928740E-07, 0.318409457190E-09, -0.560728448890E-12, 0.560750590590E-15, -0.320207200030E-18, 0.971511471520E-22, -0.121047212750E-25]] + self.coefs_exponential = [0.118597600000E+00, -0.118343200000E-03, 0.126968600000E+03] + self.coefs_splits = [3.15, 273.15, 1372.000 + 273.15] + self.description = 'Ni-Cr alloy versus Ni-Al alloy' + elif self.method == 'N': + self.coefs_piecewise = [[0.000000000000E+00, 0.261591059620E-01, 0.109574842280E-04, -0.938411115540E-07, -0.464120397590E-10, -0.263033577160E-11, -0.226534380030E-13, -0.760893007910E-16, -0.934196678350E-19], + [0.000000000000E+00, 0.259293946010E-01, 0.157101418800E-04, 0.438256272370E-07, -0.252611697940E-09, 0.643118193390E-12, -0.100634715190E-14, 0.997453389920E-18, -0.608632456070E-21, 0.208492293390E-24, -0.306821961510E-28]] + self.coefs_splits = [3.15, 273.15, 1300.000 + 273.15] + self.description = 'Ni-Cr-Si alloy versus Ni-Si-Mg alloy' + elif self.method == 'R': + self.coefs_piecewise = [[0.000000000000E+00, 0.528961729765E-02, 0.139166589782E-04, -0.238855693017E-07, 0.356916001063E-10, -0.462347666298E-13, 0.500777441034E-16, -0.373105886191E-19, 0.157716482367E-22, -0.281038625251E-26], + [0.295157925316E+01, -0.252061251332E-02, 0.159564501865E-04, -0.764085947576E-08, 0.205305291024E-11, -0.293359668173E-15], + [0.152232118209E+03, -0.268819888545E+00, 0.171280280471E-03, -0.345895706453E-07, -0.934633971046E-14]] + self.coefs_splits = [223.15, 1064.180 + 273.15, 1664.500 + 273.15, 1768.1 + 273.15] + self.description = 'Pt-13% Rh versus Pt' + elif self.method == 'S': + self.coefs_piecewise = [[0.000000000000E+00, 0.540313308631E-02, 0.125934289740E-04, -0.232477968689E-07, 0.322028823036E-10, -0.331465196389E-13, 0.255744251786E-16, -0.125068871393E-19, 0.271443176145E-23], + [0.132900444085E+01, 0.334509311344E-02, 0.654805192818E-05, -0.164856259209E-08, 0.129989605174E-13], + [0.146628232636E+03, -0.258430516752E+00, 0.163693574641E-03, -0.330439046987E-07, -0.943223690612E-14]] + self.coefs_splits = [223.15, 1064.180+273.15, 1664.5 + 273.15, 1768.1 + 273.15] + self.description = 'Pt-10% Rh versus Pt' + elif self.method == 'T': + self.coefs_piecewise = [[0.000000000000E+00, 0.387481063640E-01, 0.441944343470E-04, 0.118443231050E-06, 0.200329735540E-07, 0.901380195590E-09, 0.226511565930E-10, 0.360711542050E-12, 0.384939398830E-14, 0.282135219250E-16, 0.142515947790E-18, 0.487686622860E-21, 0.107955392700E-23, 0.139450270620E-26, 0.797951539270E-30], + [0.000000000000E+00, 0.387481063640E-01, 0.332922278800E-04, 0.206182434040E-06, -0.218822568460E-08, 0.109968809280E-10, -0.308157587720E-13, 0.454791352900E-16, -0.275129016730E-19]] + self.coefs_splits = [3.15, 273.15, 400.000 + 273.15] + self.description = 'Cu versus a Cu-Ni alloy' + + # ASTM E1751 supplementary polynomials + elif self.method == 'G': + self.coefs_piecewise = [[0.0, 0.0012792201, 2.1634754e-05, -1.1393234e-08, 4.3850022e-12, -1.7089202e-15], + [-1.1064412, 0.0094962455, -3.6467516e-06, 3.114133e-08, -3.8615222e-11, 2.4455012e-14, -8.9888053e-18, 1.8120237e-21, -1.5534591e-25]] + self.coefs_splits = [273.15, 630.615 + 273.15, 2315 + 273.15] + self.description = 'Tungsten versus Tungsten - 26% Rhenium' + elif self.method == 'P': + self.coefs_piecewise = [[0.0, 0.029819716, 3.5175152e-05, -3.4878428e-08, 1.4851327e-11, -3.6375467e-15], + [-8.9621838, 0.0853772, -0.00010570233, 1.5424937e-07, -1.2855115e-10, 5.443876e-14, -9.3211269e-18]] + self.coefs_splits = [273.15, 746.6 + 273.15, 273.15 + 1395] + self.description = 'Platinel II' + elif self.method == 'K vs AuFe 0.07': + self.coefs_piecewise = [[0.0, 0.022272367466, 3.6406179664e-06, -1.5967928202e-07, -4.5260169888e-09, 4.0432555769e-11, 4.9063035769e-12, 1.2272348484e-13, 1.6829773697e-15, 1.4636450149e-17, 8.4287909747e-20, 3.2146639387e-22, 7.8225430483e-25, 1.1010930596e-27, 6.826366158e-31]] + self.coefs_splits = [0.15, 273.15+7] + self.description = 'K vs Gold - 0.07% Iron' + elif self.method == 'PtMo 0.05 vs PtMo 0.001': + self.coefs_piecewise = [[0.0, 0.010501456, 2.8410937e-05, -4.3368594e-08, 1.058577e-10, -2.384895e-13, 3.3574252e-16, -2.0186476e-19], + [6.8354086, -0.048776479, 0.00024913353, -4.9920472e-07, 6.4615219e-10, -5.3071212e-13, 2.6865173e-16, -7.6717268e-20, 9.4670862e-24]] + self.coefs_splits = [273.15, 273.15 + 491., 273.15 + 1600] + self.description = 'Platinum - 5% Molybdenum versus Platinum - 0.1% Molybdenum' + elif self.method == 'PtRh 0.4 vs PtRh 0.2': + self.coefs_piecewise = [[0.0, 0.00036246289, 3.936032e-07, 4.2594137e-10, 1.0382985e-12, -1.5406939e-15, 1.0033974e-18, -2.849716e-22], + [-0.91201877, 0.0035246931, -3.9077442e-06, 3.6728697e-09, -1.082471e-12, 1.151628e-16, -1.261964e-20]] + self.coefs_splits = [273.15, 951.7 + 273.15, 1888 + 273.15] + self.description = 'Platinum 40% Rhodium versus Platinum - 20% Rhodium' + + elif self.method == 'M': + self.coefs_piecewise = [[0.0, 0.03690092195, 4.408522682e-05, -3.142898226e-08, -1.02521613e-10, 1.846977453e-13, -9.738054601e-17, -3.3943879e-19], + [-11.45582129, 0.2059913943, -0.0008846963426, 2.650568429e-06, -4.958763813e-09, 6.145877457e-12, -5.041679909e-15, 2.627522669e-18, -7.864442961e-22, 1.027600874e-25]] + self.coefs_splits = [-50 + 273.15, 370.8 + 273.15, 1410 + 273.15] + self.description = 'Nickel-18% Molybdenum vs Nickel-0.008% Cobalt' + elif self.method == 'IrRh 0.4 vs Ir': + self.coefs_piecewise = [[0.0, 0.0030870016, 6.9649773e-06, -7.8890504e-09, 2.7700591e-12, 2.6762413e-14, -1.041804e-16, 1.5270867e-19, -7.9634082e-23], + [-0.096839082, 0.0036588615, 5.7455189e-06, -6.0547943e-09, 2.7235393e-12, -5.1797037e-16, 3.0821886e-20]] + self.coefs_splits = [273.15, 630.615 + 273.15, 2110 + 273.15] + self.description = 'Iridium - 40% Rhodium vs Iridium' + elif self.method == 'Au-Pt': + self.coefs_piecewise = [[0.0, 0.00603619861, 1.93672974e-05, -2.22998614e-08, 3.28711859e-11, -4.24206193e-14, 4.56927038e-17, -3.39430259e-20, 1.4298159e-23, -2.51672787e-27]] + self.coefs_splits = [273.15, 1273.15] + self.description = 'Gold vs. Platinum' + elif self.method == 'Pt-Pd': + # Coefficients originally in microvolts, converted to mV + self.coefs_piecewise = [[0.0, 0.005296958, 4.610494e-06, -9.602271e-09, 2.992243e-11, -2.012523e-14, -1.268514e-17, 2.257823e-20, -8.510068e-24], + [-0.4977137, 0.010182545, -1.5793515e-05, 3.63617e-08, -2.6901509e-11, 9.5627366e-15, -1.3570737e-18]] + self.coefs_splits = [273.15, 660.323 + 273.15, 1500 + 273.15] + self.description = 'Platinum vs. Palladium' + else: + raise Exception(ValueError) + + variable_min = self.coefs_splits[0] + variable_max = self.coefs_splits[-1] + + signal_min = self.variable_to_signal(variable_min, Tref=self.Tref) + signal_max = self.variable_to_signal(variable_max, Tref=self.Tref) + + Sensor.__init__(self, variable_min=variable_min, variable_max=variable_max, + variable_units=variable_units, signal_min=signal_min, + signal_max=signal_max, signal_units=signal_units, + deadtime=deadtime, offset=offset) + + + def signal_to_variable(self, signal, Tref=None): + if Tref: + self.Tref = Tref + T = self.signal_to_variable_function(signal, self.Tref, self.coefs_splits, self.coefs_piecewise, self.coefs_exponential) + return T + def variable_to_signal(self, variable, Tref=None): + if Tref: + self.Tref = Tref +# print variable, self.Tref +# print self.coefs_splits, self.coefs_piecewise, self.coefs_exponential + E = (self.variable_to_signal_function(variable, self.coefs_splits, self.coefs_piecewise, self.coefs_exponential) + - self.variable_to_signal_function(self.Tref, self.coefs_splits, self.coefs_piecewise, self.coefs_exponential)) + return E + + def variable_function(self, variable): + E = calc_EMF(variable, self.coefs_splits, self.coefs_piecewise, self.coefs_exponential) + return E + + + + +def plot_thermocouples(): + import numpy as np + import matplotlib.pyplot as plt + from matplotlib.font_manager import FontProperties + from matplotlib.backends.backend_pdf import PdfPages + fontP = FontProperties() + fontP.set_size('small') + fig = plt.figure(figsize=(11.69,8.27)) + ax = plt.subplot(111) + for ttype in Thermocouple_list: + t = Thermocouple(method=ttype) + Tmin, Tmax = t.variable_min, t.variable_max + Es = [] + Ts = np.linspace(Tmin, Tmax, 1000) + for T in Ts: + Es.append(t.variable_function(T)) + ax.plot(Ts, Es, label=t.method) + ax.text(Ts[-1], Es[-1], t.method) + box = ax.get_position() +# ax.set_position([box.x0, box.y0, box.width * 0.9, box.height]) +# ax.legend(loc='center left', bbox_to_anchor=(1, 0.5), prop=fontP) + plt.title('EMF vs temperature for most Thermocouple types') + plt.xlabel('Temperature, K') + plt.ylabel('EMF in mV') + plt.grid(True) + ax.xaxis.grid + plt.show() + pp = PdfPages('EMF vs T.pdf') + plt.savefig(pp, format='pdf') + pp.close() +