import numpy as np
import matplotlib.pyplot as plt
from grains2 import PlaneParallelIsotropicLTE, FractalPorosity, amcarbon

D = [2.5, 2.6, 2.7, 2.8, 2.9, 3.0]
a = np.logspace(-1, 2)
ac = amcarbon()
rh = 1.5
T = np.zeros((6, len(a)))

fig, ax = plt.subplots()

for i in range(len(D)):
    porosity = FractalPorosity(0.1, D[i])
    lte = PlaneParallelIsotropicLTE(a, ac, rh, porosity=porosity)
    ax.plot(a, lte.T, label=D[i])

ax.legend()
plt.setp(
    ax,
    xlabel="Radius (μm)",
    xscale="log",
    xlim=(0.1, 100),
    ylabel="$T$",
)
plt.tight_layout()