ne_spiral_arm() - Code Metrics - Inspection of "Spiral arms" - benbaror/ne2001 - Measure and Improve Code Quality continuously with Scrutinizer

Completed

Pull Request — master (#8)

unknown

created 2017-03-04 15:18 UTC

ne_spiral_arm() F

↳ Parent: Project

Complexity

Conditions

Size

Total Lines

205

Duplication

Lines	0
Ratio	0 %

Importance

Changes	1
Bugs	0	Features	0

Metric	Value
cc	17
c	1
b	0
f	0
dl	0
loc	205
rs	2

How to fix Long Method Complexity

""" Density for spiral arms
"""
from __future__ import print_function, absolute_import, division, unicode_literals
import numpy as np
import pdb

from .utils import rad3d2
from .utils import rad2d2


def ne_spiral_arm(xyz, gal_param, adict):
    """
    Parameters
    ----------
    xyz
    gal_param
    adict

    Returns
    -------
    nea : ndarray or float
    whicharm : ndarray or int

    c-----------------------------------------------------------------------
    c  Spiral arms are defined as logarithmic spirals using the
    c    parameterization in Wainscoat et al. 1992, ApJS, 83, 111-146.
    c  But arms are modified selectively at various places to distort them
    c    as needed (08 Aug 2000).
    c  Note that arm numbering follows that of TC93 for the four large arms
    c (after remapping).
    c  The local spiral arm is number 5.
    c  06 Apr 02:   removed TC type modifications of arms 2,3 (fac calculations)
    c  		and replaced with new versions.  Data for these are hard wired.
    """
    from scipy.interpolate import CubicSpline
    x,y,z = xyz[0], xyz[1], xyz[-1]
    if isinstance(x,float):
        x = np.array([x])
        y = np.array([y])
        z = np.array([z])
        flg_float = True
    else:
        flg_float = False

    # see get_parameters for definitions of narm, warm, harm.
    narmsmax=5
    #common/armfactors/
    # .     harm(narmsmax),narm(narmsmax),warm(narmsmax),farm(narmsmax)

    #   parameter(rad=57.29577 95130 823)
    rad = 180/np.pi
    ks=3
    NN=7
    nfine = 1000

    # rr
    rr = rad2d2(xyz)
    if flg_float:
        rr = np.array([rr])

    #c
    #c Get spiral arm component:  30 do loop finds a coarse minimum distance
    #c from line of sight to arm; 40 do loop finds a fine minimum distance
    #c from line of sight to arm; line 35 ensures that arm limits are not
    #c exceeded; linear interpolation beginning at line 41 finds the
    #c minimum distance from line of sight to arm on a finer scale than gridding
    #c of arms allows (TJL)

    # Init
    whicharm = np.zeros_like(x).astype(int)
    nea = np.zeros_like(x)

    # thxy
    thxy = np.arctan2(-x, y) * rad		#! measured ccw from +y axis
                                        #! (different from tc93 theta)
    neg_th = thxy < 0.
    thxy[neg_th] += 360.
    # Cut on values near the disk
    icutz = np.where(np.abs(z/gal_param['ha']) < 10.)[0]
    if len(icutz) > 0:
        cutx = x[icutz]
        cuty = y[icutz]
        cutz = z[icutz]
        sminmin = 1.e10 * np.ones_like(cutx)
    else:
        # Time to return
        return nea
    # Find closest distance to each arm and then assign arm
    #   We will brute force with a spline
    #   Could get memory intensive
    #   Could refine
    for j in range(adict['narms']):
        #do 50 j=1,narms
        jj = adict['armmap'][j]
        # Crude
        xarm = adict['arm'][j,:adict['kmax'][j],0]
        yarm = adict['arm'][j,:adict['kmax'][j],1]

        xtmp = np.outer(cutx, np.ones_like(xarm))
        ytmp = np.outer(cuty, np.ones_like(yarm))
        xatmp = np.outer(np.ones_like(cutx), xarm)
        yatmp = np.outer(np.ones_like(cuty), yarm)
        dist = (xtmp-xatmp)**2 + (ytmp-yatmp)**2
        kmin = np.argmin(dist,axis=1)

        # Refine with a Spline? -- Requires a loop
        csp_xarm = CubicSpline(np.arange(adict['kmax'][j]),
                               adict['arm'][j,:adict['kmax'][j],0])
        csp_yarm = CubicSpline(np.arange(adict['kmax'][j]),
                               adict['arm'][j,:adict['kmax'][j],1])
        jtmp = np.zeros((len(cutx),nfine))
        for ii,ix in enumerate(cutx):
            j0 = max(0, kmin[ii]-1)
            j1 = min(adict['kmax'][j], kmin[ii]+1)
            jtmp[ii,:] = np.linspace(j0,j1,num=nfine)
        xjarm = csp_xarm(jtmp)
        yjarm = csp_yarm(jtmp)
        xtmp = np.outer(cutx, np.ones(nfine))
        ytmp = np.outer(cuty, np.ones(nfine))
        dist = (xtmp-xjarm)**2 + (ytmp-yjarm)**2
        kmin2 = np.argmin(dist,axis=1)
        min_dist2 = np.amin(dist,axis=1)

        #
        smin=np.sqrt(min_dist2)		 # ! Distance of (x,y,z) from this arm's axis
        # Close enough?
        gd_wa = np.where(smin < (3*gal_param['wa']))[0]
        if len(gd_wa) > 0:
            # ga
            ga = np.exp(-(smin[gd_wa]/(
                gal_param['warm{:d}'.format(jj)]*gal_param['wa']))**2)	#! arm, get the arm weighting factor
            # Galactocentric radial dependence of arms
            tmp_rr = rr[icutz[gd_wa]]
            lg_rr = tmp_rr > gal_param['Aa']
            if np.sum(lg_rr) > 0:
                ga[lg_rr] *= 1. / (np.cosh((tmp_rr[lg_rr]-gal_param['Aa'])/2.0))**2

            # arm3 reweighting:
            if adict['armmap'][j] == 3:
                th3a=320.
                th3b=390.
                th3b=370.
                th3a=290.
                th3b=363.
                th3b=363.
                fac3min=0.0
                test3 = thxy[icutz[gd_wa]]-th3a
                neg_t3 = test3 < 0
                test3[neg_t3] += 360.
                gd_t3 = (0. <= test3) & (test3 < (th3b-th3a))
                if np.sum(gd_t3) > 0:
                    arg = 6.2831853*(thxy[icutz[gd_wa]][gd_t3]-th3a)/(th3b-th3a)
                    fac = (1.+fac3min + (1.-fac3min)*np.cos(arg))/2.
                    fac = fac**4.0
                    # Update ga
                    ga[gd_t3] *= fac

            # arm2 reweighting:
            if adict['armmap'][j] == 2:
                th2a=35.
                th2b=55.
                test2 = thxy[icutz[gd_wa]]-th2a
                fac = 1.
                if False:
                    #    first: as in tc93 (note different definition of theta)
                    gd_t2_first =  (0 <= test2) & (test2 < (th2b-th2a))
                    if np.sum(gd_t2_first) > 0:
                        fac=1.+ test2[gd_t2_first]/(th2b-th2a)
                        fac = 1.		# !**** note turned off
                        #ga=ga*fac
                    gd_t2_first2 =  test2 > (th2b-th2a)
                    if np.sum(gd_t2_first2) > 0:
                        fac = 2.
                        fac = 1.		#!**** note turned off
                        #ga=ga*fac
                # c    second:  weaken the arm in a short range:
                th2a=340.
                th2b=370.
                # c note fix does nothing if fac2min = 1.0
                fac2min=0.1
                neg_t2 = test2 < 0.
                test2[neg_t2] += 360.
                gd_t2_snd = (0. <= test2) & (test2 < (th2b-th2a))
                if np.sum(gd_t2_snd) > 0:
                    arg=6.2831853*(thxy[icutz[gd_wa]][gd_t2_snd]-th2a)/(th2b-th2a)
                    fac = (1.+fac2min + (1.-fac2min)*np.cos(arg))/2.
                    # Update ga
                    ga[gd_t2_snd] *= fac

            # Find the ones which are closer than any previous arm
            iclosest = smin[gd_wa] < sminmin[gd_wa]
            if np.sum(iclosest) > 0:
                whicharm[icutz[gd_wa[iclosest]]] = adict['armmap'][j]

            # Update nea
            nea[icutz[gd_wa]] += gal_param['narm{:d}'.format(jj)] * gal_param['na'] * (
                ga / (np.cosh(cutz[gd_wa]/(gal_param['harm{:d}'.format(jj)]*gal_param['ha']))**2))
    if flg_float:
        nea = float(nea[0])
        whicharm = int(whicharm[0])

    # Return
    return nea # , whicharm

    '''
    Farms = 0
    if(whicharm_spiralmodel .eq. 0) then
    whicharm = 0
    else
      whicharm = armmap(whicharm_spiralmodel)	! remap arm number
      Farms = Fa * farm(whicharm)
    endif
    return
    end
    '''


1			""" Density for spiral arms
2			"""
3			from __future__ import print_function, absolute_import, division, unicode_literals
4			import numpy as np
5			import pdb
6
7			from .utils import rad3d2
8			from .utils import rad2d2
9
10
11			def ne_spiral_arm(xyz, gal_param, adict):
12			"""
13			Parameters
14			----------
15			xyz
16			gal_param
17			adict
18
19			Returns
20			-------
21			nea : ndarray or float
22			whicharm : ndarray or int
23
24			c-----------------------------------------------------------------------
25			c Spiral arms are defined as logarithmic spirals using the
26			c parameterization in Wainscoat et al. 1992, ApJS, 83, 111-146.
27			c But arms are modified selectively at various places to distort them
28			c as needed (08 Aug 2000).
29			c Note that arm numbering follows that of TC93 for the four large arms
30			c (after remapping).
31			c The local spiral arm is number 5.
32			c 06 Apr 02: removed TC type modifications of arms 2,3 (fac calculations)
33			c and replaced with new versions. Data for these are hard wired.
34			"""
35			from scipy.interpolate import CubicSpline
36			x,y,z = xyz[0], xyz[1], xyz[-1]
37			if isinstance(x,float):
38			x = np.array([x])
39			y = np.array([y])
40			z = np.array([z])
41			flg_float = True
42			else:
43			flg_float = False
44
45			# see get_parameters for definitions of narm, warm, harm.
46			narmsmax=5
47			#common/armfactors/
48			# . harm(narmsmax),narm(narmsmax),warm(narmsmax),farm(narmsmax)
49
50			# parameter(rad=57.29577 95130 823)
51			rad = 180/np.pi
52			ks=3
53			NN=7
54			nfine = 1000
55
56			# rr
57			rr = rad2d2(xyz)
58			if flg_float:
59			rr = np.array([rr])
60
61			#c
62			#c Get spiral arm component: 30 do loop finds a coarse minimum distance
63			#c from line of sight to arm; 40 do loop finds a fine minimum distance
64			#c from line of sight to arm; line 35 ensures that arm limits are not
65			#c exceeded; linear interpolation beginning at line 41 finds the
66			#c minimum distance from line of sight to arm on a finer scale than gridding
67			#c of arms allows (TJL)
68
69			# Init
70			whicharm = np.zeros_like(x).astype(int)
71			nea = np.zeros_like(x)
72
73			# thxy
74			thxy = np.arctan2(-x, y) * rad #! measured ccw from +y axis
75			#! (different from tc93 theta)
76			neg_th = thxy < 0.
77			thxy[neg_th] += 360.
78			# Cut on values near the disk
79			icutz = np.where(np.abs(z/gal_param['ha']) < 10.)[0]
80			if len(icutz) > 0:
81			cutx = x[icutz]
82			cuty = y[icutz]
83			cutz = z[icutz]
84			sminmin = 1.e10 * np.ones_like(cutx)
85			else:
86			# Time to return
87			return nea
88			# Find closest distance to each arm and then assign arm
89			# We will brute force with a spline
90			# Could get memory intensive
91			# Could refine
92			for j in range(adict['narms']):
93			#do 50 j=1,narms
94			jj = adict['armmap'][j]
95			# Crude
96			xarm = adict['arm'][j,:adict['kmax'][j],0]
97			yarm = adict['arm'][j,:adict['kmax'][j],1]
98
99			xtmp = np.outer(cutx, np.ones_like(xarm))
100			ytmp = np.outer(cuty, np.ones_like(yarm))
101			xatmp = np.outer(np.ones_like(cutx), xarm)
102			yatmp = np.outer(np.ones_like(cuty), yarm)
103			dist = (xtmp-xatmp)2 + (ytmp-yatmp)2
104			kmin = np.argmin(dist,axis=1)
105
106			# Refine with a Spline? -- Requires a loop
107			csp_xarm = CubicSpline(np.arange(adict['kmax'][j]),
108			adict['arm'][j,:adict['kmax'][j],0])
109			csp_yarm = CubicSpline(np.arange(adict['kmax'][j]),
110			adict['arm'][j,:adict['kmax'][j],1])
111			jtmp = np.zeros((len(cutx),nfine))
112			for ii,ix in enumerate(cutx):
113			j0 = max(0, kmin[ii]-1)
114			j1 = min(adict['kmax'][j], kmin[ii]+1)
115			jtmp[ii,:] = np.linspace(j0,j1,num=nfine)
116			xjarm = csp_xarm(jtmp)
117			yjarm = csp_yarm(jtmp)
118			xtmp = np.outer(cutx, np.ones(nfine))
119			ytmp = np.outer(cuty, np.ones(nfine))
120			dist = (xtmp-xjarm)2 + (ytmp-yjarm)2
121			kmin2 = np.argmin(dist,axis=1)
122			min_dist2 = np.amin(dist,axis=1)
123
124			#
125			smin=np.sqrt(min_dist2) # ! Distance of (x,y,z) from this arm's axis
126			# Close enough?
127			gd_wa = np.where(smin < (3*gal_param['wa']))[0]
128			if len(gd_wa) > 0:
129			# ga
130			ga = np.exp(-(smin[gd_wa]/(
131			gal_param['warm{:d}'.format(jj)]gal_param['wa']))*2) #! arm, get the arm weighting factor
132			# Galactocentric radial dependence of arms
133			tmp_rr = rr[icutz[gd_wa]]
134			lg_rr = tmp_rr > gal_param['Aa']
135			if np.sum(lg_rr) > 0:
136			ga[lg_rr] = 1. / (np.cosh((tmp_rr[lg_rr]-gal_param['Aa'])/2.0))*2
137
138			# arm3 reweighting:
139			if adict['armmap'][j] == 3:
140			th3a=320.
141			th3b=390.
142			th3b=370.
143			th3a=290.
144			th3b=363.
145			th3b=363.
146			fac3min=0.0
147			test3 = thxy[icutz[gd_wa]]-th3a
148			neg_t3 = test3 < 0
149			test3[neg_t3] += 360.
150			gd_t3 = (0. <= test3) & (test3 < (th3b-th3a))
151			if np.sum(gd_t3) > 0:
152			arg = 6.2831853*(thxy[icutz[gd_wa]][gd_t3]-th3a)/(th3b-th3a)
153			fac = (1.+fac3min + (1.-fac3min)*np.cos(arg))/2.
154			fac = fac**4.0
155			# Update ga
156			ga[gd_t3] *= fac
157
158			# arm2 reweighting:
159			if adict['armmap'][j] == 2:
160			th2a=35.
161			th2b=55.
162			test2 = thxy[icutz[gd_wa]]-th2a
163			fac = 1.
164			if False:
165			# first: as in tc93 (note different definition of theta)
166			gd_t2_first = (0 <= test2) & (test2 < (th2b-th2a))
167			if np.sum(gd_t2_first) > 0:
168			fac=1.+ test2[gd_t2_first]/(th2b-th2a)
169			fac = 1. # !**** note turned off
170			#ga=ga*fac
171			gd_t2_first2 = test2 > (th2b-th2a)
172			if np.sum(gd_t2_first2) > 0:
173			fac = 2.
174			fac = 1. #!**** note turned off
175			#ga=ga*fac
176			# c second: weaken the arm in a short range:
177			th2a=340.
178			th2b=370.
179			# c note fix does nothing if fac2min = 1.0
180			fac2min=0.1
181			neg_t2 = test2 < 0.
182			test2[neg_t2] += 360.
183			gd_t2_snd = (0. <= test2) & (test2 < (th2b-th2a))
184			if np.sum(gd_t2_snd) > 0:
185			arg=6.2831853*(thxy[icutz[gd_wa]][gd_t2_snd]-th2a)/(th2b-th2a)
186			fac = (1.+fac2min + (1.-fac2min)*np.cos(arg))/2.
187			# Update ga
188			ga[gd_t2_snd] *= fac
189
190			# Find the ones which are closer than any previous arm
191			iclosest = smin[gd_wa] < sminmin[gd_wa]
192			if np.sum(iclosest) > 0:
193			whicharm[icutz[gd_wa[iclosest]]] = adict['armmap'][j]
194
195			# Update nea
196			nea[icutz[gd_wa]] += gal_param['narm{:d}'.format(jj)] * gal_param['na'] * (
197			ga / (np.cosh(cutz[gd_wa]/(gal_param['harm{:d}'.format(jj)]gal_param['ha']))*2))
198			if flg_float:
199			nea = float(nea[0])
200			whicharm = int(whicharm[0])
201
202			# Return
203			return nea # , whicharm
204
205			'''
206			Farms = 0
207			if(whicharm_spiralmodel .eq. 0) then
208			whicharm = 0
209			else
210			whicharm = armmap(whicharm_spiralmodel) ! remap arm number
211			Farms = Fa * farm(whicharm)
212			endif
213			return
214			end
215			'''
216

benbaror / ne2001

Pull Request — master (#8)

ne_spiral_arm() F

Complexity

Size

Duplication

Importance

How to fix Long Method Complexity

Long Method

Complexity

Duplication Side-by-Side

Filter issues like