import math
from pfdLib import PFD

c_I =	10000
c_H =	10000000
c_R =	50000
lmbda =	0.01752
f_D =	0.2
beta =	0.1


print("Analytical solution:")
print("\nC(tau) = c_I/tau + c_R(lmbda - lmbda**2*tau/2)+ c_H*lmbda*tau/2*f_D")
print("C'(tau) = -c_I/tau**2 + c_R( - lmbda**2/2)+ c_H*lmbda/2*f_D=0")
print("==>tau* = sqrt(2*c_I/(lmbda*(f_D*c_H-lmbda*c_R)))")
print("\ntau* =",math.sqrt(2*c_I/(lmbda*(f_D*c_H-lmbda*c_R))))

def C(tau,lmbda,f_D,c_I,c_R,c_H):
    return c_I/tau+c_R*(lmbda-lmbda**2*tau/2)+c_H*lmbda*tau/2*f_D

from scipy.optimize import minimize    
tau0 = 1.1 # Initial guess

result = minimize(C, tau0, args = (lmbda,f_D,c_I,c_R,c_H))
print("\nNumerical solution from scipy minimize:")
print("tau* =",result.x[0])
print("C(tau*) =",result.fun)

import matplotlib.pyplot as plt
from numpy import arange
# Create empty vectors to hold values to plot
tau_i  = []
c_I_i= []
c_R_i= []
c_F_i= []
c_T_i=[]
for tau in arange (0.25,2.5,0.1):
    tau_i.append(tau)
    c_I_i.append(c_I/tau)
    c_R_i.append(c_R*(lmbda-lmbda**2*tau/2))
    c_F_i.append(c_H*lmbda*tau/2*f_D)
    c_T_i.append(c_I/tau+c_R*(lmbda-lmbda**2*tau/2)+c_H*lmbda*tau/2*f_D)

plt.plot (tau_i , c_I_i, label= 'Inspection')
plt.plot (tau_i , c_R_i, label= 'Repair')
plt.plot (tau_i , c_F_i, label= 'Hazard')
plt.plot (tau_i , c_T_i, label= 'Total')
plt.xlabel (r'$\tau$ (years)')
plt.ylabel (r'$C (\tau )$')
plt.title ('Cost as function of maintenance interval ')
plt.legend()
#plt.savefig ('output.svg')
plt.show ()