init_spiral_arms() - Code Metrics - Inspection of "Galactic sightline script (and more)" - benbaror/ne2001 - Measure and Improve Code Quality continuously with Scrutinizer

Completed

Pull Request — master (#4)

unknown

created 2017-02-26 14:03 UTC

init_spiral_arms() B

↳ Parent: Project

Complexity

Conditions

Size

Total Lines

Duplication

Lines	0
Ratio	0 %

Importance

Changes	1
Bugs	0	Features	0

Metric	Value
cc	6
c	1
b	0
f	0
dl	0
loc	82
rs	7.3557

How to fix Long Method

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

import numpy as np
import json
import pdb

import ne2001
data_path = ne2001.__path__[0]+'/data/'


def read_params(ifile='ne2001_params.json'):
    """ Read parameter file"""
    # Read
    with open(data_path+ifile, 'rt') as fh:
        PARAMS = json.load(fh)
    # Recast?
    return PARAMS


def read_galparam(ifile='gal_param.json'):
    """ Read Galaxy parameters
    Parameters
    ----------
    ifile : str, optional

    Returns
    -------
    gal_param : dict

    """
    # Read
    with open(data_path+ifile, 'rt') as fh:
        galparam_dict = json.load(fh)
    # Thick disk
    thick_dict = dict(e_density=galparam_dict['n1h1']/galparam_dict['h1'],
                      height=galparam_dict['h1'],
                      radius=galparam_dict['A1'],
                      F=galparam_dict['F1'],
                      )
    # Thin disk
    thin_dict = dict(e_density=galparam_dict['n2'],
                      height=galparam_dict['h2'],
                      radius=galparam_dict['A2'],
                      F=galparam_dict['F2'],
                      )
    # Return
    return galparam_dict


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

    """
    # Read
    with open(data_path+ifile, 'rt') as fh:
        gc_dict = json.load(fh)
    # Return
    return gc_dict


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

    Returns
    -------
    lism_dict : dict

    """
    # Read
    with open(data_path+ifile, 'rt') as fh:
        lism_dict = json.load(fh)
    # Return
    return lism_dict


def init_spiral_arms():

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

benbaror / ne2001

Pull Request — master (#4)

init_spiral_arms() B

Complexity

Size

Duplication

Importance

How to fix Long Method

Long Method

Duplication Side-by-Side

Filter issues like