agndisks/testmod.py at master · ajdittmann/agndisks · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
import numpy as np
import models


from scipy import interpolate as scint
opacTable = "combined-opacity/opacitysolar09dustq3p5amax0p1new.txt"
with open(opacTable,'r') as fp:
    line=fp.readline()
    numbers_str = line.split()
    nrho=int(numbers_str[1])
    nt=int(numbers_str[3])
    rhoread=np.zeros(nrho)
    treadnew=np.zeros(nt)
    rosscombineread=np.zeros([nt,nrho])
    planckcombineread=np.zeros([nt,nrho])
    for irho in range(nrho):
        for it in range(nt):
            line=fp.readline()
            numbers_str = line.split()
            rhoread[irho]=float(numbers_str[0])
            treadnew[it]=float(numbers_str[1])
            rosscombineread[it,irho]=float(numbers_str[2])
    opacityinter = scint.RectBivariateSpline(np.log10(rhoread), np.log10(treadnew), np.log10(rosscombineread.T), kx=3, ky=3, s=0)
    maxT = np.max(treadnew)

def opacTabulated(rho, T):
    if T > maxT: return kes
    else: return 10.0**opacityinter.ev(np.log10(rho), np.log10(T))


#thisMod = SG03() # opacityFunction = opacTabulated)
thisMod = models.SG03( opacityFunction = opacTabulated)
#thisMod = SS73( opacityFunction = opacTabulated)
test = thisMod.calcModel(7.0, 0.1, 0.1, Nr = 10**4)

import matplotlib.pyplot as plt
#plt.semilogx(test['radius']/test['Rg'], test['badcells'])
#plt.show()
plt.loglog(test['radius']/test['Rg'], test['temp'])
#plt.show()
plt.loglog(test['radius']/test['Rg'], test['kappa'])
#plt.show()
plt.loglog(test['radius']/test['Rg'], test['tau'])
#plt.show()
plt.loglog(test['radius']/test['Rg'], test['toomreQ'])
plt.show()