-
Notifications
You must be signed in to change notification settings - Fork 1
/
tinti.py
executable file
·119 lines (87 loc) · 4.44 KB
/
tinti.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
import numpy as np
def i1(t,Tr):
p1 = np.pi*(4*t - Tr)*Tr + 8*t**1.5*np.sqrt(-t + Tr)
p2 = 4*Tr*np.sqrt(t*(-t + Tr))
with np.errstate(divide='ignore'):
p3 = - 2*(4*t - Tr)*Tr*np.arctan((-2*t + Tr)/(2.*np.sqrt(t*(-t + Tr))))
return (p1+p2+p3)/16.0
def i2(t,Tr,Ts):
p1 = (Tr*np.sqrt(t*(-t + Tr)) + 2*np.sqrt(t**3*(-t + Tr)) - t*Tr*np.sqrt((-t + Tr + Ts)/(t - Ts)) +
Tr*Ts*np.sqrt((-t + Tr + Ts)/(t - Ts)) - 2*t*np.sqrt((t - Ts)*(-t + Tr + Ts)) -
2*Ts*np.sqrt((t - Ts)*(-t + Tr + Ts)) + (4*t - Tr)*Tr*np.arcsin(np.sqrt(t/Tr)) +
Tr*(-4*t + Tr)*np.arcsin(np.sqrt((t - Ts)/Tr)))/4.
p2 = (-((2*t + Tr - 6*Ts)*np.sqrt((t - Ts)*(-t + Tr + Ts))) +
Tr*(-4*t + Tr + 8*Ts)*np.arcsin(np.sqrt((t - Ts)/Tr)))/4.
return p1+p2
def i3(t,Tr,Ts):
p1 = (Tr*np.sqrt(t*(-t + Tr)) + 2*np.sqrt(t**3*(-t + Tr)) - t*Tr*np.sqrt((-t + Tr + Ts)/(t - Ts)) +
Tr*Ts*np.sqrt((-t + Tr + Ts)/(t - Ts)) - 2*t*np.sqrt((t - Ts)*(-t + Tr + Ts)) -
2*Ts*np.sqrt((t - Ts)*(-t + Tr + Ts)) + (4*t - Tr)*Tr*np.arcsin(np.sqrt(t/Tr)) +
Tr*(-4*t + Tr)*np.arcsin(np.sqrt((t - Ts)/Tr)))/4.
p2 = (2*(-2*t*np.sqrt((t - Ts)*(-t + Tr + Ts)) - Tr*np.sqrt((t - Ts)*(-t + Tr + Ts)) +
6*Ts*np.sqrt((t - Ts)*(-t + Tr + Ts)) + 2*t*np.sqrt((t - 2*Ts)*(-t + Tr + 2*Ts)) +
Tr*np.sqrt((t - 2*Ts)*(-t + Tr + 2*Ts)) - 4*Ts*np.sqrt((t - 2*Ts)*(-t + Tr + 2*Ts))) -
Tr*(-4*t + Tr + 8*Ts)*np.arctan((-2*t + Tr + 2*Ts)/(2.*np.sqrt((t - Ts)*(-t + Tr + Ts)))) +
Tr*(-4*t + Tr + 8*Ts)*np.arctan((-2*t + Tr + 4*Ts)/(2.*np.sqrt((t - 2*Ts)*(-t + Tr + 2*Ts)))))/8.
return p1+p2
def i4(t,Tr,Ts):
p1 = t*Tr*np.arccos(np.sqrt((t - Ts)/Tr))
p2 = (-(np.sqrt((t - Ts)*(-t + Tr + Ts))*(2*t + Tr + 2*Ts)) -
Tr**2*np.arccos(1.0/(np.sqrt(Tr/(t - Ts)))))/4.
p3 = (2*(-2*t*np.sqrt((t - Ts)*(-t + Tr + Ts)) - Tr*np.sqrt((t - Ts)*(-t + Tr + Ts)) +
6*Ts*np.sqrt((t - Ts)*(-t + Tr + Ts)) + 2*t*np.sqrt((t - 2*Ts)*(-t + Tr + 2*Ts)) +
Tr*np.sqrt((t - 2*Ts)*(-t + Tr + 2*Ts)) - 4*Ts*np.sqrt((t - 2*Ts)*(-t + Tr + 2*Ts))) -
Tr*(-4*t + Tr + 8*Ts)*np.arctan((-2*t + Tr + 2*Ts)/(2.*np.sqrt((t - Ts)*(-t + Tr + Ts)))) +
Tr*(-4*t + Tr + 8*Ts)*np.arctan((-2*t + Tr + 4*Ts)/(2.*np.sqrt((t - 2*Ts)*(-t + Tr + 2*Ts)))))/8.
return p1+p2+p3
def i5(t,Tr,Ts):
p1 = (4*(2*t + Tr - 4*Ts)*np.sqrt(-((t - 2*Ts)*(t - Tr - 2*Ts))) + np.pi*Tr*(-4*t + Tr + 8*Ts) +
2*Tr*(-4*t + Tr + 8*Ts)*np.arctan((-2*t + Tr + 4*Ts)/(2.*np.sqrt((t - 2*Ts)*(-t + Tr + 2*Ts)))))/16.
return p1
def si3(t, Tr, Ts):
p1 = t*Tr*np.arccos(np.sqrt((t - Ts)/Tr)) + (-(np.sqrt((t - Ts)*(-t + Tr + Ts))*(2*t + Tr + 2*Ts)) -
Tr**2*np.arccos(1.0/np.sqrt(Tr/(t - Ts))))/4.
p2 = (-4*(2*t + Tr - 6*Ts)*np.sqrt((t - Ts)*(-t + Tr + Ts)) + np.pi*Tr*(-4*t + Tr + 8*Ts) -
2*Tr*(-4*t + Tr + 8*Ts)*np.arctan((-2*t + Tr + 2*Ts)/(2.*np.sqrt((t - Ts)*(-t + Tr + Ts)))))/16.
return p1+p2
def tinti(t, Ts, Tr, t0):
n = len(t)
dt = t[1]-t[0]
# normalizing constant
k = 2 / (np.pi * Tr * Ts**2)
stf = np.zeros(n)
if Tr > 2*Ts:
stf[(t >= 0) & (t <= Ts)] = i1(t[(t >= 0) & (t <= Ts)], Tr)
stf[(t > Ts) & (t < 2*Ts)] = i2(t[(t > Ts) & (t < 2*Ts)], Tr, Ts)
stf[(t >= 2*Ts) & (t < Tr)] = i3(t[(t >= 2*Ts) & (t < Tr)], Tr, Ts)
stf[(t >= Tr) & (t < Tr+Ts)] = i4(t[(t >= Tr) & (t < Tr+Ts)], Tr, Ts)
stf[(t >= Tr+Ts) & (t < Tr + 2*Ts)] = i5(t[(t >= Tr+Ts) & (t < Tr + 2*Ts)], Tr, Ts)
elif Tr > Ts and Tr <= 2*Ts:
stf[(t >= 0) & (t <= Ts)] = i1(t[(t >= 0) & (t <= Ts)], Tr)
stf[(t > Ts) & (t < Tr)] = i2(t[(t > Ts) & (t < Tr)], Tr, Ts)
stf[(t >= Tr) & (t < 2*Ts)] = si3(t[(t >= Tr) & (t < 2*Ts)], Tr, Ts)
stf[(t >= 2*Ts) & (t < Ts+Tr)] = i4(t[(t >= 2*Ts) & (t < Ts+Tr)], Tr, Ts)
stf[(t >= Ts+Tr) & (t<2*Ts+Tr)] = i5(t[(t >= Ts+Tr) & (t<2*Ts+Tr)], Tr, Ts)
else:
print(f"Invalid parameters assigned. (tr={Tr:.3f}, ts={Ts:.3f})")
d = int(np.ceil((t0)/dt))
stf = np.roll(stf, d)
return k*stf
if __name__ == "__main__":
# testing goes here.
from pylab import *
dt = 0.002
tt = 5.0
t = arange(0, tt, dt)
ts = 0.1
tr = 1.0
t0 = 0.5
stf = tinti(t, ts, tr, t0)
figure()
plot(t, stf)
ts = 0.6
tr = 1.0
t0 = 0.5
stf = tinti(t, ts, tr, t0)
plot(t, stf)
show()