https://github.com/QTB-HHU/ModelHeatShock
Tip revision: 7ead643f56ee9bfa3077d5f0c05f0912cbf9a4ee authored by StefanoMagni on 05 April 2018, 15:38:57 UTC
Delete .DS_Store
Delete .DS_Store
Tip revision: 7ead643
HSM_Temperature.py
from matplotlib.pylab import *
############ TEMPERATURE ############
# IMPORTANT: note difference with the original T(t) using steps, due to .ode crashing for t>220s.
# Once the new shape for the temperature is used .ode works fine, and this is more physical.
# But be aware of the difference!!!
def T(t, TparamsSet):
#print(t/60.)
Ttype = TparamsSet.CurrentParams["Ttype"] # 1 --> Tin, Tup 2 --> Tin, Tup, Tin 3 --> Tin, Tup, Tin, Tup
Tin = TparamsSet.CurrentParams["Tin"]
# REMARK!!! tau should be ways smaller than (ti-tj) approx few 1000s.
# Around 1s to 10s. Smaller -> .ode crashes, bigger -> to big difference with step T(t).
# Best value would be 7s for all the code (but then canavnine gives problems)
if Ttype == 0:
TT = Tin
elif Ttype == 1:
Tup = TparamsSet.CurrentParams["Tup"]
DT = Tup - Tin
ta = TparamsSet.CurrentParams["ta"]
tau = TparamsSet.CurrentParams["tau"]
if t <= ta:
TT = Tin
elif t > ta:
TT = Tin + DT * (1. - exp(-(t - ta) / tau))
else:
print("Error in Temperature function!")
elif Ttype == 2:
Tup = TparamsSet.CurrentParams["Tup"]
DT = Tup - Tin
ta = TparamsSet.CurrentParams["ta"]
tb = TparamsSet.CurrentParams["tb"]
tau = TparamsSet.CurrentParams["tau"]
if t <= ta:
TT = Tin
elif t > ta and t <= tb:
TT = Tin + DT * (1. - exp(-(t - ta) / tau))
elif t > tb:
TT = Tin + DT * (exp(-(t - tb) / tau))
else:
print("Error in Temperature function!")
elif Ttype == 3:
Tup = TparamsSet.CurrentParams["Tup"]
DT = Tup - Tin
ta = TparamsSet.CurrentParams["ta"]
tb = TparamsSet.CurrentParams["tb"]
tc = TparamsSet.CurrentParams["tc"]
tau = TparamsSet.CurrentParams["tau"]
if t <= ta:
TT = Tin
elif t > ta and t <= tb:
TT = Tin + DT * (1. - exp(-(t - ta) / tau))
elif t > tb and t <= tc:
TT = Tin + DT * (exp(-(t - tb) / tau))
elif t > tc:
TT = Tin + DT * (1. - exp(-(t - tc) / tau))
else:
print("Error in Temperature function!")
elif Ttype == 4:
Tup = TparamsSet.CurrentParams["Tup"]
DT = Tup - Tin
ta = TparamsSet.CurrentParams["ta"]
tb = TparamsSet.CurrentParams["tb"]
tc = TparamsSet.CurrentParams["tc"]
td = TparamsSet.CurrentParams["td"]
tau = TparamsSet.CurrentParams["tau"]
if t <= ta:
TT = Tin
elif t > ta and t <= tb:
TT = Tin + DT * (1. - exp(-(t - ta) / tau))
elif t > tb and t <= tc:
TT = Tin + DT * (exp(-(t - tb) / tau))
elif t > tc and t <= td:
TT = Tin + DT * (1. - exp(-(t - tc) / tau))
elif t > td:
TT = Tin + DT * (exp(-(t - td) / tau))
else:
print("Error in Temperature function!")
elif Ttype == 5:
Tup1 = TparamsSet.CurrentParams["Tup1"]
Tup2 = TparamsSet.CurrentParams["Tup2"]
DT1 = Tup1 - Tin
DT2 = Tup2 - Tup1
t1 = TparamsSet.CurrentParams["t1"]
t2 = TparamsSet.CurrentParams["t2"]
tau = TparamsSet.CurrentParams["tau"]
if t <= t1:
TT = Tin
elif t > t1 and t <= t2:
TT = Tin + DT1 * (1. - exp(-(t - t1) / tau))
elif t > t2:
TT = Tup1 + DT2 * (1. - exp(-(t - t2) / tau))
else:
print("Error in Temperature function!")
elif Ttype == 6:
Tup = TparamsSet.CurrentParams["Tup"]
DT = Tup - Tin
tUP = TparamsSet.CurrentParams["Dt"]
NumHS = TparamsSet.CurrentParams["NumHS"]
tau = TparamsSet.CurrentParams["tau"]
TT = Tin
number = int((t-vorl)/tUP)
if t >= vorl:
if number == 0:
pass
elif number % 2 == 0:
TT = Tin + DT * (exp(-(t - (vorl + (number) * tUP)) / tau))
else:
TT = Tin + DT * (1. - exp(-(t - (vorl + (number) * tUP)) / tau))
#for i in range(0, NumHS):
# if t > vorl + (1 + 2 * i) * tUP and t <= vorl + (2 + 2 * i) * tUP:
# TT = Tin + DT * (1. - exp(-(t - (vorl + (1 + 2 * i) * tUP)) / tau))
# elif t > vorl + (2 + 2 * i) * tUP and t <= vorl + (3 + 2 * i) * tUP:
# TT = Tin + DT * (exp(-(t - (vorl + (2 + 2 * i) * tUP)) / tau))
elif Ttype == 7:
Tup = TparamsSet.CurrentParams["Tup"]
Period = TparamsSet.CurrentParams["Period"]
# Phase = TparamsSet.CurrentParams["Phase"]
DT = Tup - Tin
if t < vorl:
TT = Tin # Tin+DT/2.+DT/2.*sin( (-Phase)/Period*2*3.14152 - 3.14152/2.)
elif t >= vorl:
TT = Tin + DT / 2. + DT / 2. * sin((t - vorl) / Period * 2 * math.pi - math.pi / 2.)
else:
print("Error in Temperature function!")
elif Ttype == 8:
Tup = TparamsSet.CurrentParams["Tup"]
DT = Tup - Tin
ta = TparamsSet.CurrentParams["ta"]
tau = TparamsSet.CurrentParams["tau"]
#print("!!! "+str(t)+" "+str(tau)+"\n")
if t <= ta:
TT = Tin
elif t > ta and t <= (ta + tau):
TT = Tin + 0.5*DT + 0.5*DT*math.cos((t-ta)/tau*math.pi-math.pi)
elif t > (ta + tau):
TT = Tup
else:
print("Error in Temperature function!")
else:
print("Error, wrong Temperature index!")
#print( str(Ttype) + " " + str((t-vorl)/60.) + " " + str(TT) )
return TT
vorl = (2500.) * 60. # (changed from 1 hour=60*60 in Alex's code!!!)
# time after which we start to perturb/look at the system (i.e. we assume steady state has been reached)
def MakeGrayBackgroundTemperature(ax, timeset, TparamsSet):
""" Put gray background to the plots when the temperature is high (Heat Shock) """
tend = (timeset.CurrentParams["t_stop"] - vorl) / 60.
Ttype = TparamsSet.CurrentParams["Ttype"]
ColorOfFace = (255. / 255, 220. / 255., 203. / 255.) # (0.8, 0.8, 0.8)
ColorOfEdge = (255. / 255, 210. / 255., 200. / 255.) # (0.8, 0.8, 0.8)
if Ttype == 3:
ta = (TparamsSet.CurrentParams["ta"] - vorl) / 60.
tb = (TparamsSet.CurrentParams["tb"] - vorl) / 60.
tc = (TparamsSet.CurrentParams["tc"] - vorl) / 60.
ax.axvspan(ta, tb, facecolor=ColorOfFace, edgecolor=ColorOfEdge, linewidth=0.8, alpha=1)
ax.axvspan(tc, tend, facecolor=ColorOfFace, edgecolor=ColorOfEdge, linewidth=0.8, alpha=1)
elif Ttype == 2:
ta = (TparamsSet.CurrentParams["ta"] - vorl) / 60.
tb = (TparamsSet.CurrentParams["tb"] - vorl) / 60.
ax.axvspan(ta, tb, facecolor=ColorOfFace, edgecolor=ColorOfEdge, linewidth=0.8, alpha=1)
elif Ttype == 1:
ta = (TparamsSet.CurrentParams["ta"] - vorl) / 60.
ax.axvspan(ta, tend, facecolor=ColorOfFace, edgecolor=ColorOfEdge, linewidth=0.8, alpha=1)
elif Ttype == 4:
ta = (TparamsSet.CurrentParams["ta"] - vorl) / 60.
tb = (TparamsSet.CurrentParams["tb"] - vorl) / 60.
tc = (TparamsSet.CurrentParams["tc"] - vorl) / 60.
td = (TparamsSet.CurrentParams["td"] - vorl) / 60.
ax.axvspan(ta, tb, facecolor=ColorOfFace, edgecolor=ColorOfEdge, linewidth=0.8, alpha=1)
ax.axvspan(tc, td, facecolor=ColorOfFace, edgecolor=ColorOfEdge, linewidth=0.8, alpha=1)
elif Ttype == 0 or Ttype == 5 or Ttype == 6 or Ttype == 7 or Ttype == 8:
pass
else:
print("Error, wrong Temperature index in function PlotTrajectories!")