-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathbound.R
123 lines (82 loc) · 3.37 KB
/
bound.R
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
120
121
122
source("routinesHO.R")
###############################################################################
# Parameter values
###############################################################################
eta <- 0.5
w0 <- 1.2
hb <- 0.5
h0 <- 0.1
r0 <- 0.5
w1 <- w0*sqrt(abs(eta^2-1))
A <- rep2ap(eta, w0, hb, h0, r0)$A
phi <- rep2ap(eta, w0, hb, h0, r0)$phi
mu <- w1/(w0*eta)
###############################################################################
# Under-damped case (eta < 1)
###############################################################################
K1 <- (abs(A)/(w0*eta))*(1+mu)/(1+mu^2)
K2 <- w0*eta
K = exp( -(A/(w0*eta))*( sin(phi) + mu*cos(phi) )/(mu^2+1) )
bound <- Vectorize( function(t) K*exp(-hb*t + K1*exp(-K2*t)) )
survf <- Vectorize(function(t) exp(-chHO(t, eta, w0, hb, h0, r0)) )
curve(bound,0,5, lwd = 2, n = 1000)
curve(survf, 0, 5, lwd = 2, lty = 2, col = "red", add = TRUE, n = 1000)
curve(bound,20,25, lwd = 2)
curve(survf, 20, 25, lwd = 1, lty = 2, col = "red", add = TRUE)
curve(bound,20,21, lwd = 2)
curve(survf, 20, 21, lwd = 1, lty = 2, col = "red", add = TRUE)
###############################################################################
# Parameter values
###############################################################################
eta <- 1.5
w0 <- 1.2
hb <- 0.5
h0 <- 0.1
r0 <- 0.5
w1 <- w0*sqrt(abs(eta^2-1))
A <- rep2ap(eta, w0, hb, h0, r0)$A
phi <- rep2ap(eta, w0, hb, h0, r0)$phi
mu <- w1/(w0*eta)
a <- (h0-hb)/mu + r0/w1
###############################################################################
# Over-damped case (eta > 1)
###############################################################################
K1 <- max( c((hb-h0-a)/(w1-w0*eta), (h0-hb-a)/(w1+w0*eta) ))
K2 <- w0*eta - w1
K = exp( -0.5*(hb-h0-a)/(w1-w0*eta) - 0.5*(h0-hb-a)/(w1+w0*eta) )
bound <- Vectorize( function(t) K*exp(-hb*t + K1*exp(-K2*t)) )
survf <- Vectorize(function(t) exp(-chHO(t, eta, w0, hb, h0, r0)) )
curve(bound,0,5, lwd = 2, n = 1000, ylim=c(0,2))
curve(survf, 0, 5, lwd = 2, lty = 2, col = "red", add = TRUE, n = 1000)
curve(bound,20,25, lwd = 2)
curve(survf, 20, 25, lwd = 1, lty = 2, col = "red", add = TRUE)
curve(bound,20,21, lwd = 2)
curve(survf, 20, 21, lwd = 1, lty = 2, col = "red", add = TRUE)
###############################################################################
# Parameter values
###############################################################################
eta <- 1
w0 <- 1.2
hb <- 0.5
h0 <- 0.1
r0 <- 0.5
w1 <- w0*sqrt(abs(eta^2-1))
A <- rep2ap(eta, w0, hb, h0, r0)$A
phi <- rep2ap(eta, w0, hb, h0, r0)$phi
mu <- w1/(w0*eta)
a <- (h0-hb)/mu + r0/w1
###############################################################################
# Critically-damped case (eta = 1)
###############################################################################
K1 <- max( c((h0-hb)/w0, (r0 + w0*(h0-hb))/(w0^2) ))
K2 <- w0
K3 <- (r0 + w0*(h0-hb))/w0
K = exp( -(h0-hb)/w0 -(r0 + w0*(h0-hb))/(w0^2) )
bound <- Vectorize( function(t) K*exp(-hb*t + K1*exp(-K2*t) + K3*t*exp(-K2*t) ) )
survf <- Vectorize(function(t) exp(-chHO(t, eta, w0, hb, h0, r0)) )
curve(bound,0,5, lwd = 2, n = 1000, ylim=c(0,2))
curve(survf, 0, 5, lwd = 2, lty = 2, col = "red", add = TRUE, n = 1000)
curve(bound,20,25, lwd = 2)
curve(survf, 20, 25, lwd = 1, lty = 2, col = "red", add = TRUE)
curve(bound,20,21, lwd = 2)
curve(survf, 20, 21, lwd = 1, lty = 2, col = "red", add = TRUE)