|
1
|
|
|
""" Module for I/O |
|
2
|
|
|
""" |
|
3
|
|
|
from __future__ import absolute_import, division, print_function |
|
4
|
|
|
|
|
5
|
|
|
import numpy as np |
|
6
|
|
|
import json |
|
7
|
|
|
import pdb |
|
8
|
|
|
|
|
9
|
|
|
import ne2001 |
|
10
|
|
|
data_path = ne2001.__path__[0]+'/data/' |
|
11
|
|
|
|
|
12
|
|
|
|
|
13
|
|
|
def read_params(ifile='ne2001_params.json'): |
|
14
|
|
|
""" Read parameter file""" |
|
15
|
|
|
# Read |
|
16
|
|
|
with open(data_path+ifile, 'rt') as fh: |
|
17
|
|
|
PARAMS = json.load(fh) |
|
18
|
|
|
# Convert lists to array (maybe should save these as YAML not JSON) |
|
19
|
|
|
for key in PARAMS: |
|
20
|
|
|
if isinstance(PARAMS[key], dict): |
|
21
|
|
|
for key2,item in PARAMS[key].items(): |
|
22
|
|
|
if isinstance(item,list): |
|
23
|
|
|
PARAMS[key][key2] = np.array(item) |
|
24
|
|
|
# Add |
|
25
|
|
|
PARAMS['spiral_arms']['adict'] = init_spiral_arms() |
|
26
|
|
|
PARAMS['spiral_arms']['gal_param'] = read_galparam() |
|
27
|
|
|
# Recast? |
|
28
|
|
|
return PARAMS |
|
29
|
|
|
|
|
30
|
|
|
|
|
31
|
|
|
def read_galparam(ifile='gal_param.json'): |
|
32
|
|
|
""" Read Galaxy parameters |
|
33
|
|
|
Parameters |
|
34
|
|
|
---------- |
|
35
|
|
|
ifile : str, optional |
|
36
|
|
|
|
|
37
|
|
|
Returns |
|
38
|
|
|
------- |
|
39
|
|
|
gal_param : dict |
|
40
|
|
|
|
|
41
|
|
|
""" |
|
42
|
|
|
# Read |
|
43
|
|
|
with open(data_path+ifile, 'rt') as fh: |
|
44
|
|
|
galparam_dict = json.load(fh) |
|
45
|
|
|
# Thick disk |
|
46
|
|
|
thick_dict = dict(e_density=galparam_dict['n1h1']/galparam_dict['h1'], |
|
47
|
|
|
height=galparam_dict['h1'], |
|
48
|
|
|
radius=galparam_dict['A1'], |
|
49
|
|
|
F=galparam_dict['F1'], |
|
50
|
|
|
) |
|
51
|
|
|
# Thin disk |
|
52
|
|
|
thin_dict = dict(e_density=galparam_dict['n2'], |
|
53
|
|
|
height=galparam_dict['h2'], |
|
54
|
|
|
radius=galparam_dict['A2'], |
|
55
|
|
|
F=galparam_dict['F2'], |
|
56
|
|
|
) |
|
57
|
|
|
# Return |
|
58
|
|
|
return galparam_dict |
|
59
|
|
|
|
|
60
|
|
|
|
|
61
|
|
|
def read_gc(ifile='ne_gc.json'): |
|
62
|
|
|
""" Read Galactic Center parameters |
|
63
|
|
|
Returns |
|
64
|
|
|
------- |
|
65
|
|
|
gc_dict : dict |
|
66
|
|
|
dict of parameters |
|
67
|
|
|
|
|
68
|
|
|
""" |
|
69
|
|
|
# Read |
|
70
|
|
|
with open(data_path+ifile, 'rt') as fh: |
|
71
|
|
|
gc_dict = json.load(fh) |
|
72
|
|
|
# Return |
|
73
|
|
|
return gc_dict |
|
74
|
|
|
|
|
75
|
|
|
|
|
76
|
|
|
def read_lism(ifile='ne_lism.json'): |
|
77
|
|
|
""" |
|
78
|
|
|
Parameters |
|
79
|
|
|
---------- |
|
80
|
|
|
ifile : str, optional |
|
81
|
|
|
|
|
82
|
|
|
Returns |
|
83
|
|
|
------- |
|
84
|
|
|
lism_dict : dict |
|
85
|
|
|
|
|
86
|
|
|
""" |
|
87
|
|
|
# Read |
|
88
|
|
|
with open(data_path+ifile, 'rt') as fh: |
|
89
|
|
|
lism_dict = json.load(fh) |
|
90
|
|
|
# Return |
|
91
|
|
|
return lism_dict |
|
92
|
|
|
|
|
93
|
|
|
|
|
94
|
|
|
def init_spiral_arms(): |
|
95
|
|
|
|
|
96
|
|
|
from astropy.table import Table |
|
97
|
|
|
from scipy.interpolate import CubicSpline |
|
98
|
|
|
armsinp= data_path + 'ne_arms_log_mod.inp' |
|
99
|
|
|
logarms= data_path + 'log_arms.out' |
|
100
|
|
|
|
|
101
|
|
|
narms=5 |
|
102
|
|
|
#integer armmap(5) ! for remapping from Wainscoat |
|
103
|
|
|
#data armmap/1, 3, 4, 2, 5/ ! order to TC93 order, which is |
|
104
|
|
|
# ! from GC outwards toward Sun. |
|
105
|
|
|
armmap = [1,3,4,2,5] |
|
106
|
|
|
NNj = [20, 20, 20, 20, 20] |
|
107
|
|
|
narmpoints=500 |
|
108
|
|
|
ncoord=2 |
|
109
|
|
|
NNmax=20 |
|
110
|
|
|
rad = 180/np.pi |
|
111
|
|
|
# Arms |
|
112
|
|
|
arms_tbl = Table.read(armsinp, format='ascii') # a, rmin, thmin, extent |
|
113
|
|
|
assert len(arms_tbl) == narms |
|
114
|
|
|
|
|
115
|
|
|
r1 = np.zeros((NNmax, narms)) |
|
116
|
|
|
th1 = np.zeros((NNmax, narms)) |
|
117
|
|
|
kmax = np.zeros(narms).astype(int) |
|
118
|
|
|
arm = np.zeros((narms, narmpoints, ncoord)) |
|
119
|
|
|
|
|
120
|
|
|
for j, row in enumerate(arms_tbl): |
|
121
|
|
|
th1[0:NNj[j],j] = row['thmin'] + np.arange(NNj[j])*row['extent']/(NNj[j]-1.) #! rad |
|
122
|
|
|
r1[:,j] = row['rmin'] * np.exp((th1[:,j]-row['thmin'])/row['a']) |
|
123
|
|
|
th1[:,j] *= rad # ! deg |
|
124
|
|
|
#c *** begin sculpting spiral arm 2 == TC arm 3*** |
|
125
|
|
|
if armmap[j] == 3: |
|
126
|
|
|
cut1 = (th1[:,j] > 370.) & (th1[:,j] <= 410.) |
|
127
|
|
|
r1[cut1,j] *= (1. + 0.04* np.cos((th1[cut1,j]-390.)*180./(40.*rad))) |
|
128
|
|
|
#c . (1. + 0.01*cos((th1(n,j)-390.)*180./(40.*rad))) |
|
129
|
|
|
cut2 = (th1[:,j] > 315.) & (th1[:,j] <= 370.) |
|
130
|
|
|
r1[cut2,j] *= (1. - 0.07* np.cos((th1[cut2,j]-345.)*180./(55.*rad))) |
|
131
|
|
|
#c . (1.0 - 0.08*cos((th1(n,j)-345.)*180./(55.*rad))) |
|
132
|
|
|
cut3 = (th1[:,j] > 180.) & (th1[:,j] <= 315.) |
|
133
|
|
|
r1[cut3,j] *= (1 + 0.16* np.cos((th1[cut3,j]-260.)*180./(135.*rad))) |
|
134
|
|
|
# (1 + 0.13* np.cos((th1[cut3,j]-260.)*180./(135.*rad))) |
|
135
|
|
|
#c *** begin sculpting spiral arm 4 == TC arm 2*** |
|
136
|
|
|
if armmap[j] == 2: |
|
137
|
|
|
cut1 = (th1[:,j] > 290.) & (th1[:,j] <= 395.) |
|
138
|
|
|
r1[cut1,j] *= (1. - 0.11* np.cos((th1[cut1,j]-350.)*180./(105.*rad))) |
|
139
|
|
|
#c *** end arm sculpting *** |
|
140
|
|
|
|
|
141
|
|
|
|
|
142
|
|
|
|
|
143
|
|
|
""" |
|
144
|
|
|
open(11,file=logarms, status='unknown') |
|
145
|
|
|
write(11,*) 'arm n xa ya' |
|
146
|
|
|
""" |
|
147
|
|
|
#do 21 j=1,narms |
|
148
|
|
|
for j in range(narms): |
|
149
|
|
|
dth = 5.0/r1[0,j] # Python indexing |
|
150
|
|
|
th = th1[0,j]-0.999*dth |
|
151
|
|
|
# Generate spline |
|
152
|
|
|
cspline = CubicSpline(th1[:NNj[j],j],r1[:NNj[j],j]) |
|
153
|
|
|
#call cspline(th1(1,j),r1(1,j),-NNj(j),th,r) |
|
154
|
|
|
#for k in range(narmpoints): |
|
155
|
|
|
#do 10 k=1,narmpoints-1 |
|
156
|
|
|
th = th + dth * np.arange(narmpoints) |
|
157
|
|
|
gd_th = np.where(th <= th1[NNj[j]-1, j])[0] |
|
158
|
|
|
kmax[j] = np.max(gd_th) + 1 # Python indexing (we will use arange) |
|
159
|
|
|
r = cspline(th[gd_th]) |
|
160
|
|
|
# x,y of each arm |
|
161
|
|
|
arm[j,gd_th,0] = -r*np.sin(th[gd_th]/rad) # Python indexing |
|
162
|
|
|
arm[j,gd_th,1] = r*np.cos(th[gd_th]/rad) |
|
163
|
|
|
|
|
164
|
|
|
# Wrap into a dict |
|
165
|
|
|
arms_dict = {} |
|
166
|
|
|
arms_dict['table'] = arms_tbl |
|
167
|
|
|
arms_dict['r1'] = r1 |
|
168
|
|
|
arms_dict['th1'] = r1 |
|
169
|
|
|
arms_dict['kmax'] = kmax |
|
170
|
|
|
arms_dict['narms'] = narms |
|
171
|
|
|
arms_dict['narmpoints'] = narmpoints |
|
172
|
|
|
arms_dict['armmap'] = armmap |
|
173
|
|
|
arms_dict['arm'] = arm |
|
174
|
|
|
return arms_dict |
|
175
|
|
|
|
benbaror /
ne2001
