init_spiral_arms() - Code Metrics - benbaror/ne2001 - Measure and Improve Code Quality continuously with Scrutinizer

init_spiral_arms() B
last analyzed 2017-04-25 12:30 UTC

↳ Parent: Project

Complexity

Conditions

Size

Total Lines

Duplication

Lines	0
Ratio	0 %

Importance

Changes	3
Bugs	0	Features	0

Metric	Value
cc	6
c	3
b	0
f	0
dl	0
loc	81
rs	7.3726

How to fix Long Method

""" Module for I/O
"""
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import json
import os
from builtins import super

import numpy as np
from astropy.table import Table
from scipy.interpolate import CubicSpline

from . import __path__

DATA_PATH = os.path.join(__path__[0], 'data')


def numpify_dict(d):
    """
    Recursively make lists in a dictionary into numpy array
    """
    def numpify(d):
        for k, v in d.items():
            if isinstance(v, list):
                d[k] = np.array(v)
            elif isinstance(v, dict):
                numpify(v)
    new_dict = d.copy()
    numpify(new_dict)
    return new_dict


class Params(dict):
    """
    Input parameters
    """

    def __init__(self, ifile='ne2001_params.json', path=None, **new_params):
        """
        """
        if path is None:
            path = DATA_PATH
        self.path = path
        self.ifile = ifile
        try:
            params = numpify_dict(parse_json(os.path.join(self.path,
                                                          self.ifile)))
            params['spiral_arms']['adict'] = init_spiral_arms()
        except IOError:
            params = {}
        params.update(new_params)
        super().__init__(params)


def parse_json(json_file):
    "Parse json file"
    with open(json_file, 'rt') as json_data:
        data = json.load(json_data)
    return data


def read_galparam(ifile='gal_param.json'):
    """
    Read Galaxy parameters

    Parameters
    ----------
    ifile : str, optional

    Returns
    -------
    gal_param : dict

    """
    old_param = parse_json(os.path.join(DATA_PATH, ifile))
    gal_param = {}

    gal_param['thick_disk'] = dict(e_density=(old_param['n1h1'] /
                                              old_param['h1']),
                                   height=old_param['h1'],
                                   radius=old_param['A1'],
                                   F=old_param['F1'])

    gal_param['thin_disk'] = dict(e_density=old_param['n2'],
                                  height=old_param['h2'],
                                  radius=old_param['A2'],
                                  F=old_param['F2'])

    return gal_param


def read_gc(ifile='ne_gc.json'):
    """ Read Galactic Center parameters
    Returns
    -------
    gc_param : dict
      dict of parameters

    """
    old_param = parse_json(os.path.join(DATA_PATH, ifile))
    gc_param = {}

    gc_param['galactic_center'] = dict(e_density=old_param['negc0'],
                                       center=tuple(old_param['centroid'].
                                                    values()),
                                       F=old_param['Fgc0'],
                                       height=old_param['hgc'],
                                       radius=old_param['rgc'])

    return gc_param


def read_lism(ifile='ne_lism.json'):
    """
    Parameters
    ----------
    ifile : str, optional

    Returns
    -------
    lism_dict : dict

    """
    # Read
    with open(os.path.join(DATA_PATH, ifile), 'rt') as fh:
        lism_dict = json.load(fh)
    # Return
    return lism_dict


def init_spiral_arms(ifile='ne_arms_log_mod.inp'):
    armsinp = os.path.join(DATA_PATH, ifile)
    # logarms = DATA_PATH + 'log_arms.out'

    narms = 5
    # integer armmap(5)		! for remapping from Wainscoat
    # data armmap/1, 3, 4, 2, 5/	! order to TC93 order, which is
    #                ! from GC outwards toward Sun.
    armmap = [1, 3, 4, 2, 5]
    NNj = [20, 20, 20, 20, 20]
    narmpoints = 500
    ncoord = 2
    NNmax = 20
    rad = 180/np.pi
    # Arms
    arms_tbl = Table.read(armsinp, format='ascii')  # a, rmin, thmin, extent
    assert len(arms_tbl) == narms

    r1 = np.zeros((NNmax, narms))
    th1 = np.zeros((NNmax, narms))
    kmax = np.zeros(narms).astype(int)
    arm = np.zeros((narms, narmpoints, ncoord))

    for j, row in enumerate(arms_tbl):
        th1[0:NNj[j], j] = (row['thmin'] +
                            np.arange(NNj[j])*row['extent']/(NNj[j]-1.))  # rad
        r1[:, j] = row['rmin'] * np.exp((th1[:, j]-row['thmin'])/row['a'])
        th1[:, j] *= rad  # ! deg
        # c *** begin sculpting spiral arm 2 == TC arm 3***
        if armmap[j] == 3:
            cut1 = (th1[:, j] > 370.) & (th1[:, j] <= 410.)
            r1[cut1, j] *= (1. + 0.04 * np.cos((th1[cut1, j]-390.)*180 /
                                               (40.*rad)))
            # c . (1. + 0.01*cos((th1(n,j)-390.)*180./(40.*rad)))
            cut2 = (th1[:, j] > 315.) & (th1[:, j] <= 370.)
            r1[cut2, j] *= (1. - 0.07 * np.cos((th1[cut2, j]-345.)*180 /
                                               (55.*rad)))
            # c   .   (1.0 - 0.08*cos((th1(n,j)-345.)*180./(55.*rad)))
            cut3 = (th1[:, j] > 180.) & (th1[:, j] <= 315.)
            r1[cut3, j] *= (1 + 0.16 * np.cos((th1[cut3, j]-260.)*180 /
                                              (135.*rad)))
            # (1 + 0.13* np.cos((th1[cut3,j]-260.)*180./(135.*rad)))
        # c *** begin sculpting spiral arm 4 == TC arm 2***
        if armmap[j] == 2:
            cut1 = (th1[:, j] > 290.) & (th1[:, j] <= 395.)
            r1[cut1, j] *= (1. - 0.11 * np.cos((th1[cut1, j]-350.)*180 /
                                               (105.*rad)))
        # c *** end arm sculpting ***

    """
    open(11,file=logarms, status='unknown')
        write(11,*) 'arm  n   xa     ya'
    """
    # do 21 j=1,narms
    for j in range(narms):
        dth = 5.0/r1[0, j]  # Python indexing
        th = th1[0, j]-0.999*dth
        # Generate spline
        cspline = CubicSpline(th1[:NNj[j], j], r1[:NNj[j], j])
        # call cspline(th1(1,j),r1(1,j),-NNj(j),th,r)
        # for k in range(narmpoints):
        # do 10 k=1,narmpoints-1
        th = th + dth * np.arange(narmpoints)
        gd_th = np.where(th <= th1[NNj[j]-1, j])[0]
        kmax[j] = np.max(gd_th) + 1  # Python indexing (we will use arange)
        r = cspline(th[gd_th])
        # x,y of each arm
        arm[j, gd_th, 0] = -r*np.sin(th[gd_th]/rad)  # Python indexing
        arm[j, gd_th, 1] = r*np.cos(th[gd_th]/rad)

    # Wrap into a dict
    arms_dict = {}
    arms_dict['table'] = arms_tbl
    arms_dict['r1'] = r1
    arms_dict['th1'] = r1
    arms_dict['kmax'] = kmax
    arms_dict['narms'] = narms
    arms_dict['narmpoints'] = narmpoints
    arms_dict['armmap'] = armmap
    arms_dict['arm'] = arm
    return arms_dict


1			""" Module for I/O
2			"""
3			from __future__ import absolute_import
4			from __future__ import division
5			from __future__ import print_function
6
7			import json
8			import os
9			from builtins import super
10
11			import numpy as np
12			from astropy.table import Table
13			from scipy.interpolate import CubicSpline
14
15			from . import __path__
16
17			DATA_PATH = os.path.join(__path__[0], 'data')
18
19
20			def numpify_dict(d):
21			"""
22			Recursively make lists in a dictionary into numpy array
23			"""
24			def numpify(d):
25			for k, v in d.items():
26			if isinstance(v, list):
27			d[k] = np.array(v)
28			elif isinstance(v, dict):
29			numpify(v)
30			new_dict = d.copy()
31			numpify(new_dict)
32			return new_dict
33
34
35			class Params(dict):
36			"""
37			Input parameters
38			"""
39
40			def __init__(self, ifile='ne2001_params.json', path=None, **new_params):
41			"""
42			"""
43			if path is None:
44			path = DATA_PATH
45			self.path = path
46			self.ifile = ifile
47			try:
48			params = numpify_dict(parse_json(os.path.join(self.path,
49			self.ifile)))
50			params['spiral_arms']['adict'] = init_spiral_arms()
51			except IOError:
52			params = {}
53			params.update(new_params)
54			super().__init__(params)
55
56
57			def parse_json(json_file):
58			"Parse json file"
59			with open(json_file, 'rt') as json_data:
60			data = json.load(json_data)
61			return data
62
63
64			def read_galparam(ifile='gal_param.json'):
65			"""
66			Read Galaxy parameters
67
68			Parameters
69			----------
70			ifile : str, optional
71
72			Returns
73			-------
74			gal_param : dict
75
76			"""
77			old_param = parse_json(os.path.join(DATA_PATH, ifile))
78			gal_param = {}
79
80			gal_param['thick_disk'] = dict(e_density=(old_param['n1h1'] /
81			old_param['h1']),
82			height=old_param['h1'],
83			radius=old_param['A1'],
84			F=old_param['F1'])
85
86			gal_param['thin_disk'] = dict(e_density=old_param['n2'],
87			height=old_param['h2'],
88			radius=old_param['A2'],
89			F=old_param['F2'])
90
91			return gal_param
92
93
94			def read_gc(ifile='ne_gc.json'):
95			""" Read Galactic Center parameters
96			Returns
97			-------
98			gc_param : dict
99			dict of parameters
100
101			"""
102			old_param = parse_json(os.path.join(DATA_PATH, ifile))
103			gc_param = {}
104
105			gc_param['galactic_center'] = dict(e_density=old_param['negc0'],
106			center=tuple(old_param['centroid'].
107			values()),
108			F=old_param['Fgc0'],
109			height=old_param['hgc'],
110			radius=old_param['rgc'])
111
112			return gc_param
113
114
115			def read_lism(ifile='ne_lism.json'):
116			"""
117			Parameters
118			----------
119			ifile : str, optional
120
121			Returns
122			-------
123			lism_dict : dict
124
125			"""
126			# Read
127			with open(os.path.join(DATA_PATH, ifile), 'rt') as fh:
128			lism_dict = json.load(fh)
129			# Return
130			return lism_dict
131
132
133			def init_spiral_arms(ifile='ne_arms_log_mod.inp'):
134			armsinp = os.path.join(DATA_PATH, ifile)
135			# logarms = DATA_PATH + 'log_arms.out'
136
137			narms = 5
138			# integer armmap(5) ! for remapping from Wainscoat
139			# data armmap/1, 3, 4, 2, 5/ ! order to TC93 order, which is
140			# ! from GC outwards toward Sun.
141			armmap = [1, 3, 4, 2, 5]
142			NNj = [20, 20, 20, 20, 20]
143			narmpoints = 500
144			ncoord = 2
145			NNmax = 20
146			rad = 180/np.pi
147			# Arms
148			arms_tbl = Table.read(armsinp, format='ascii') # a, rmin, thmin, extent
149			assert len(arms_tbl) == narms
150
151			r1 = np.zeros((NNmax, narms))
152			th1 = np.zeros((NNmax, narms))
153			kmax = np.zeros(narms).astype(int)
154			arm = np.zeros((narms, narmpoints, ncoord))
155
156			for j, row in enumerate(arms_tbl):
157			th1[0:NNj[j], j] = (row['thmin'] +
158			np.arange(NNj[j])*row['extent']/(NNj[j]-1.)) # rad
159			r1[:, j] = row['rmin'] * np.exp((th1[:, j]-row['thmin'])/row['a'])
160			th1[:, j] *= rad # ! deg
161			# c * begin sculpting spiral arm 2 == TC arm 3*
162			if armmap[j] == 3:
163			cut1 = (th1[:, j] > 370.) & (th1[:, j] <= 410.)
164			r1[cut1, j] = (1. + 0.04 np.cos((th1[cut1, j]-390.)*180 /
165			(40.*rad)))
166			# c . (1. + 0.01cos((th1(n,j)-390.)180./(40.*rad)))
167			cut2 = (th1[:, j] > 315.) & (th1[:, j] <= 370.)
168			r1[cut2, j] = (1. - 0.07 np.cos((th1[cut2, j]-345.)*180 /
169			(55.*rad)))
170			# c . (1.0 - 0.08cos((th1(n,j)-345.)180./(55.*rad)))
171			cut3 = (th1[:, j] > 180.) & (th1[:, j] <= 315.)
172			r1[cut3, j] = (1 + 0.16 np.cos((th1[cut3, j]-260.)*180 /
173			(135.*rad)))
174			# (1 + 0.13* np.cos((th1[cut3,j]-260.)180./(135.rad)))
175			# c * begin sculpting spiral arm 4 == TC arm 2*
176			if armmap[j] == 2:
177			cut1 = (th1[:, j] > 290.) & (th1[:, j] <= 395.)
178			r1[cut1, j] = (1. - 0.11 np.cos((th1[cut1, j]-350.)*180 /
179			(105.*rad)))
180			# c * end arm sculpting *
181
182			"""
183			open(11,file=logarms, status='unknown')
184			write(11,*) 'arm n xa ya'
185			"""
186			# do 21 j=1,narms
187			for j in range(narms):
188			dth = 5.0/r1[0, j] # Python indexing
189			th = th1[0, j]-0.999*dth
190			# Generate spline
191			cspline = CubicSpline(th1[:NNj[j], j], r1[:NNj[j], j])
192			# call cspline(th1(1,j),r1(1,j),-NNj(j),th,r)
193			# for k in range(narmpoints):
194			# do 10 k=1,narmpoints-1
195			th = th + dth * np.arange(narmpoints)
196			gd_th = np.where(th <= th1[NNj[j]-1, j])[0]
197			kmax[j] = np.max(gd_th) + 1 # Python indexing (we will use arange)
198			r = cspline(th[gd_th])
199			# x,y of each arm
200			arm[j, gd_th, 0] = -r*np.sin(th[gd_th]/rad) # Python indexing
201			arm[j, gd_th, 1] = r*np.cos(th[gd_th]/rad)
202
203			# Wrap into a dict
204			arms_dict = {}
205			arms_dict['table'] = arms_tbl
206			arms_dict['r1'] = r1
207			arms_dict['th1'] = r1
208			arms_dict['kmax'] = kmax
209			arms_dict['narms'] = narms
210			arms_dict['narmpoints'] = narmpoints
211			arms_dict['armmap'] = armmap
212			arms_dict['arm'] = arm
213			return arms_dict
214

benbaror / ne2001

init_spiral_arms() B last analyzed 2017-04-25 12:30 UTC

Complexity

Size

Duplication

Importance

How to fix Long Method

Long Method

Duplication Side-by-Side

Filter issues like

init_spiral_arms() B
last analyzed 2017-04-25 12:30 UTC