tugas 3a

import numpy as np
import matplotlib.pyplot as plt

#Data soal
Nb=4e4
rhos=2
Nr=11
Rspan=np.linspace(0.1,0.6,Nr)
W=1e5
lamda=100
Cp=1
kx=0.01
theta=200
n=40
tol=1e-6
T0guess=330
eps=0.1

#Sub-Routine 1
def fun(To,ri,mo):
    m=4/3*np.pi*ri**3*rhos*Nb
    T=To+lamda*(mo-m)/(W+mo)*Cp
    x=(mo-m)/W
    xs=np.exp(8.8053-3333/T)
    return 1/(xs-x)

#Sub-routine 2
def funfun(To,Ro):
    mo=4/3*np.pi*Ro**3*rhos*Nb
    Rispan=np.linspace(0,Ro,n+1)
    dr=Rispan[1]-Rispan[0]
    f=np.zeros(n+1)
    fsim=np.zeros(n+1)
    for i in range(n+1):
        f[i]=fun(To,Rispan[i],mo)
        if i==1 or i==n:
            fsim[i]=f[i]
        elif (-1)**i<0:
            fsim[i]=4*f[i]
        else:
            fsim[i]=2*f[i]
    thetacalc=rhos*dr/3/kx*sum(fsim)
    return (theta-thetacalc)


#Sub-routine 3
def funfunfun(Ro):
    Told=T0guess
    eror=1
    while eror>tol:
        fold=funfun(Told,Ro)
        dfold=(funfun(Told+eps,Ro)-funfun(Told-eps,Ro))/2/eps
        Tnew=Told-fold/dfold
        eror=np.abs((Told-Tnew)/Told)
        Told=Tnew
    return Tnew

#Main Program
Tstore=np.zeros(Nr)
Fx=np.zeros(Nr)
for i in range(Nr):
    Tstore[i]=funfunfun(Rspan[i])
    Fx[i]=funfun(Tstore[i],Rspan[i])

#Tabel
tabel=np.zeros([Nr,3])
tabel[:,0],tabel[:,1],tabel[:,2]=Rspan,Tstore,Fx
garis='='*40
print(garis)
header=(['R0 (cm)','T0 opt (K)','Nilai F(x)'])
print('{:^12s}{:^15s}{:^15s}'.format(*header))
print(garis)
for baris in tabel:
    print('{:^12.2f}{:^15.2f}{:^15.2e}'.format(*baris))
print(garis)

#plotting
plt.plot(Rspan,Tstore,'-or')
plt.grid()
plt.ylabel('Nilai T0 optimum (K)')
plt.xlabel('jejari awal partikel (cm)')
plt.show()