random_code_snippets

def windacclaw(ur,a,uh,x):
    Gamma = 0.1 
    m  = 1     #point mass
    y  = a*a   #fixed solid angle
    dur_da = (Gamma * y * np.exp(-x) * (1. - ur/uh)**2 - m) /a**2/2./ur
    return dur_da
    
    
    
 #initial condition
ur = 1.e-5
a = np.linspace(1.,10.,200)
uh_arr = [2., 5., 10.]

plt.figure(figsize=(12,6))
ls = ['-', '--', '-.']
uh = 5.
for x in [-3,-4,-5]:
    i=0
    for uh in uh_arr:
        Uax = odeint(windacclaw, ur, a, args=(uh,x))
        plt.plot(a, Uax, label=r'$x=$%.f'%(x), ls=ls[i])
        i+=1

plt.ylim(0.0,5.0)
plt.xlim(1.0, 10.)
# plt.legend()
plt.xlabel('a')
plt.ylabel('Ua(x)')



###Get time from acceleration law
uh = 5
v0 = 500 * kmps
r0 = kpc/2.
x  = -5
a  = np.linspace(1,10.,200)
Uax = odeint(windacclaw, ur, a, args=(uh,x))
da  = a[1] - a[0]
ur_arr = np.zeros(a.shape[0]-1)
t = np.zeros(a.shape[0]-1)
dt = np.zeros(a.shape[0]-1)
i=0
ur_prev = ur
t_prev  = 0.0
for acc in Uax[1:]:
    ur_arr[i] = ur_prev + da * acc
    ur_prev = ur_arr[i] 
    t[i] = t_prev + da/ur_prev
    dt[i] = da/ur_prev
    t_prev = t[i]
    i+=1
    
###plot time

t_unit = r0/v0
t_phy = t*t_unit/Myr

Sigma_gas = 10. #Msun pc^-2
t_diss = 7. * (np.exp(x)*Sigma_gas/10.) #in Myr

plt.figure(figsize=(12,6))
plt.plot(a[1:]*r0/kpc,t_phy)
plt.axhline(t_diss, label=r'$t_\rm{diss}$', ls='--')
plt.xlim(r0/kpc, 10.*r0/kpc)
# plt.ylim(0.0,20.)
plt.xlabel('r [kpc]')
plt.ylabel('t [Myr]')