MRN  = graindist.sizedist.Powerlaw(amax=AMAX)
ECUT = graindist.sizedist.ExpCutoff(acut=AMAX)
SIL  = graindist.composition.CmSilicate(rho=RHO)

gd_mrn  = graindist.GrainDist(MRN, SIL, md=MD, custom=True)
gd_ecut = graindist.GrainDist(ECUT, SIL, md=MD, custom=True)

ax = plt.subplot(111)
gd_mrn.plot(ax, color='k', lw=2, label='Powerlaw')
gd_ecut.plot(ax, color='b', lw=2, alpha=0.8, label='Exponential cut-off')
plt.legend(loc='upper left', frameon=False)

"""
Explanation: Making a GrainDist object from scratch
End of explanation
"""


gd_mrn2 = graindist.GrainDist('Powerlaw', 'Silicate', amax=AMAX, rho=RHO, md=MD)
gd_ecut2 = graindist.GrainDist('ExpCutoff', 'Silicate', amax=AMAX, rho=RHO, md=MD)

ax = plt.subplot(111)
gd_mrn.plot(ax, color='k', lw=2, label='')
gd_mrn2.plot(ax, color='r', lw=2, ls='--', label='Powerlaw')

gd_ecut.plot(ax, color='k', lw=2, label='')
gd_ecut2.plot(ax, color='b', lw=2, ls='--', label='Exponential Cutoff')
plt.legend(loc='upper left', frameon=False)

"""
Explanation: Making a GrainDist object with helper function
This is a shortcut for making some common grain distributions
You can change AMAX and RHO from here, too
For the different grain size distributions, amax acts in a different way:


graindist.sizedist.Grain: amax sets the singular grain size


graindist.sizedist.Powerlaw: amax sets the maximum grain size in the distribution


graindist.sizedist.ExpCutoff: amax sets the "acut" value
End of explanation
"""


gd_mrn3 = graindist.GrainDist('Powerlaw', 'Silicate', md=MD)
gd_ecut3 = graindist.GrainDist('ExpCutoff', 'Silicate', md=MD)

ax = plt.subplot(111)
gd_mrn.plot(ax, color='k', lw=2, label='')
gd_mrn3.plot(ax, color='r', lw=2, ls='--', label='Powerlaw')

gd_ecut.plot(ax, color='k', lw=2, label='')
gd_ecut3.plot(ax, color='b', lw=2, ls='--', label='Exponential Cutoff')
plt.legend(loc='upper left', frameon=False)

"""
Explanation: Note that the values I've used in this example are different from the defaults
Silicate has a default grain material density of 3.8 g cm^-3
End of explanation
"""


plt.plot(gd_mrn.a, gd_mrn.mdens, label='Powerlaw')
plt.plot(gd_ecut.a, gd_ecut.mdens, label='ExpCutoff')
plt.loglog()
plt.xlabel('Radius (um)')
plt.ylabel('Mass density (g cm$^{-2}$ um$^{-1}$)')
plt.legend(loc='upper right', frameon=False)

"""
Explanation: You can also get the mass density of dust grains
But the GrainDist.plot function only does number density
End of explanation
"""


gd_single = graindist.GrainDist('Grain', 'Silicate', amax=AMAX)

gd_single.plot()

"""
Explanation: You can't plot single grain sizes
It prints the number density instead
End of explanation
"""


