agndisks/grid.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
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
import numpy as np
import models

from scipy import interpolate as scint
#opacTable = "combined-opacity/opacitysolar09dustq3p5amax0p1new.txt"
opacTable = "combined-opacity/opacitysolar09M1dustq3p5amax10new.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)
  R, T = np.meshgrid(np.log10(rhoread), np.log10(treadnew))
  R = R.flatten()
  T = T.flatten()
  #opacityinter = scint.SmoothBivariateSpline(R, T, np.log10(rosscombineread.T).flatten(), kx=3, ky=3)
  maxT = np.max(treadnew)

def opacTabulated(rho, T):
    if T > maxT: T = maxT
    if T < 1: T = 1.0
    if rho < 0: rho = 10**-24

    #else: return 10.0**opacityinter.ev(np.log10(rho), np.log10(T))
    return 10.0**opacityinter.ev(np.log10(rho), np.log10(T))

#sg = models.SG03()
sg = models.SG03( opacityFunction = opacTabulated)

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

#Nr = 3*10**3
#Nbh = 30

Nr = 1*10**4
Nbh = 200


mbh = np.geomspace(6.0, 9.0, Nbh)
data = np.empty((Nbh, Nr, 10))

rhos = np.empty((Nr, Nbh))
taus = np.empty((Nr, Nbh))
for i in range(Nbh):
  model = sg.calcModel(mbh[i], 0.1, 0.1, Nr = Nr, rout = 10**7.0, rin=7)
  rhos[:,i] = model['dens']
  taus[:,i] = np.log10(model['tau']) #prad']/model['pgas'])
  data[i,:,0] = model['radius']
  data[i,:,1] = model['omega']
  data[i,:,2] = model['dens']
  data[i,:,3] = model['temp']
  data[i,:,4] = model['height']
  data[i,:,5] = model['tau']
  data[i,:,6] = model['kappa']
  data[i,:,7] = model['prad']
  data[i,:,8] = model['pgas']
  data[i,:,9] = model['cs']

np.savez('grid_alpha0.1_eta0.1.npz', data=data, mbh=mbh)

import matplotlib.pyplot as plt
plt.pcolor(np.log10(model['radius']/model['Rg']), mbh, np.log10(rhos.T))
plt.show()
plt.pcolor(np.log10(model['radius']/model['Rg']), mbh, taus.T, vmin = -np.max(np.abs(taus)), vmax = np.max(np.abs(taus)), cmap='RdBu' )
plt.show()