import numpy as np
import matplotlib.pyplot as plt
#366 for 0.1 (H) ,1,50,1200
#158 for 2 (H) ,1,50,1200
#values
D0= 20 #[micron]
K0= 10**10
n= 0.5
Q= 221000
R = 8.314
HR= 2 #input. vary 0.1 - 2.0
CR= 1 #input
t_iso = 50 #input
T_isoC = 1200 #input. vary 900-1400 [Celsius]
T_iso = T_isoC + 273 #[Kelvin]
delta_T = T_iso - 298
delta_t = delta_T/(100*HR) #100 steps with equal length for heating
dirta_t = delta_T/(100*CR) #100 steps with equal length for cooling
sum_before_iso = -np.exp(-Q/(R*298))*delta_t #Eventually the integral from 0 to t_1
sum_after_iso = -np.exp(-Q/(R*T_iso))*dirta_t #Eventually the integral from t_2 to t_3
for i in range(0,101):
    T_before = 298 + HR*delta_t*i
    sum_i = np.exp(-Q/(R*T_before))*delta_t
    sum_before_iso += sum_i
    T_after = T_iso - CR*dirta_t*i
    Sum_i = np.exp(-Q/(R*T_after))*dirta_t
    sum_after_iso += Sum_i

a = D0
b = K0*sum_before_iso
c = K0*np.exp(-Q/(R*T_iso))*t_iso
d = K0*sum_after_iso

D = (a**(1/n)+b+c+d)**n
print(a,b,c,d)
print(D)



#t_iso_list = []
#print(x1)
#y1 = T_before_iso
#y2 = [T_iso] * 10
#y2 = [T_iso]*int(t_iso*HR/delta_T)
#y3 = T_after_iso
#y = [*y1,*y2,*y3]
#print(y)
#plt.plot(x1,y1)
#plt.plot(x2,y2)
#plt.show()

#if sub == 1:
    #HR = input('Choose HR')
#elif sub == 2:
    #T_iso = input('Choose T_iso')
#elif sub == 3:
    #das
#elif sub == 4:
    #sdasd