##########################################################################
#                                                                        #
#    This program plots h(p(t)) as a function of t for a simple          #
#    2-out-of-3 system with Weibull-distributed component lifetimes.     #
#                                                                        #
##########################################################################

from math import *
import matplotlib.pyplot as plt
import numpy as np

# Time parameters
max_t = 120
steps = 120

# Weibull-parameters for the components
alpha = [2.0, 2.5, 3.0]
beta = [50.0, 60.0, 70.0]

# Component reliability vector
num_comps = 3
p = np.zeros(num_comps)

# System reliabilities
h = np.zeros(steps)


def reliability(pp):
    return pp[0] * pp[1] + pp[0] * pp[2] + pp[1] * pp[2] - 2 * pp[0] * pp[1] * pp[2]

# The survival function of a Weibull-distribution
def weibull(aa, bb, t):
    if t >= 0:
        return exp(- (t / bb) ** aa)
    else:
        return 1.0


time = np.linspace(0.0, max_t, steps)

for step in range(steps):
    for i in range(num_comps):
        p[i] = weibull(alpha[i], beta[i], time[step])
    h[step] = reliability(p)

fig = plt.figure(figsize = (7, 4))
plt.plot(time, h, label='h(p(t))')
plt.xlabel('Time')
plt.ylabel('Reliability')
plt.title("System reliability as a function of time")
plt.legend()
plt.show()