Temperature Test 1



In [1]:

    
%matplotlib inline

import disktemperature as dt
import numpy             as np
import matplotlib.pyplot as plt

from matplotlib.pyplot import style
from disktemperature.constants import R_sun,T_sun,M_sun,AU,pi,Grav,k_b,mu,m_p

style.use(['seaborn-dark',{'axes.grid': True,'font.size':18,'figure.figsize':(10,8),'lines.linewidth':2}]);

Define a disk model.



In [2]:

    
r        = np.logspace(np.log10(3e-2),np.log10(500),500)*AU
M_star  = M_sun
R_star  = R_sun
T_star  = T_sun
M_disk  = 0.1*M_star
alpha   = 1e-2
r_c     = 60*AU
eps     = 0.01
sig_g   = M_disk/(2*pi*r_c**2)*(r_c/r)*np.exp(-(r/r_c))
sig_g   = M_disk*sig_g/np.trapz(2*pi*r*sig_g,x=r)
sig_d_t = sig_g*eps

Initialize the mid-plane temperature model.



In [3]:

    
Temp = dt.tmid.tmid(_na=50,_nl=60,T_star=T_star,R_star=R_star,M_star=M_star,Tmin=10.)









    



Mie ... Done!



In [4]:

    
f,ax = plt.subplots()
#cols = brewer2mpl.get_map('Set2','Qualitative',8).mpl_colors
cols = plt.rcParams['axes.prop_cycle'].by_key()['color']

for iamax,amax in enumerate([1e-3,1e-1]):

    # Create dust size distribution
    
    a    = Temp.a*1e-4
    dist = (a*(a<amax))**0.5
    dist = dist/np.sum(dist)+1e-100
    sig_d = np.outer(dist,sig_d_t)
    sig_d = sig_d/sig_d.sum(0)*sig_d_t

    # calculate & plot temperature
    
    for ialpha,alpha in enumerate([0,1e-3,1e-2,1e-1]):

        T = Temp.get_t_mid(r,sig_g,sig_d,alpha*np.ones(r.shape))
        ax.loglog(r/AU,T,color=cols[ialpha],ls='-'+'-'*(iamax==1),\
                  label=r'$\alpha = $ {:3.2g}'.format(alpha)*(iamax==0))

    ax.plot(np.nan,np.nan,color='k',ls='-'+'-'*(iamax==1),label='$a_\mathrm{{max}}$ = {:3.2g} cm'.format(amax))
    l=ax.legend()
    ax.set_xlabel('$r$ [AU]')
    ax.set_ylabel('$T$ [K]')
    ax.set_ylim(ymin=1);

Try to iterate to set the irradiation angle with a radial dependence.



In [5]:

    
alpha    = 1e-3
phi      = 0.05
om       = np.sqrt(Grav*M_star/r**3)
T_of_phi = lambda phi: Temp.get_t_mid(r,sig_g,sig_d,alpha*np.ones(r.shape),phi=phi)
H_of_T   = lambda T:   np.sqrt(T*(k_b/mu/m_p))/om
T        = dt.tmid.temperature_iterator(r,T_of_phi,H_of_T,do_plot=True,phi_min=0.001)