# HW 34
import numpy as np
from scipy import interpolate 
from scipy import integrate
import matplotlib.pyplot as plt

def qfun(x):
    A = (5+15.5*x-19.5*x**2)/1000
    T = 100*(1-np.sinh(x))
    return T*A
    
def invinterp(x,y,ystar):
    n = len(x)
    xstar = []
    numxstar= 0
    for i in range(n-1):
        if (y[i]-ystar)*(y[i+1]-ystar) < 0:
          numxstar = numxstar+1;
          xs = [x[i-1],x[i],x[i+1],x[i+2]]
          ys = [y[i-1],y[i],y[i+1],y[i+2]]          
          S = interpolate.splrep(ys,xs);
          temp = interpolate.splev(ystar,S,der=0)
          print('temp = ',temp)
          xstar.append(temp)    
    return xstar   
    
 
rho = 8830
c = 385
n = 21
x = np.linspace(0,1,n)
Q = []
for i in range(n):
    temp = integrate.quad(qfun,0,x[i])
    Q.append(temp[0]*rho*c)

Qtotal = Q[n-1]
ystar1 = 0.17*Qtotal
xstar1 = invinterp(x,Q,ystar1)
print('x* for 17% = ',xstar1)

ystar2 = 0.46*Qtotal
xstar2 = invinterp(x,Q,ystar2)
print('x* for 47% = ',xstar2) 

ystar3 = 0.83*Qtotal
xstar3 = invinterp(x,Q,ystar3)
print('x* for 84% = ',xstar3)

xs = [xstar1, xstar2, xstar3]
ys = [ystar1, ystar2, ystar3]

plt.figure(1)
plt.plot(x,Q,'r-')
plt.plot(x,Q,'bo')
plt.plot(xs,ys,'k*')
plt.plot