Skip to content

Commit

Permalink
Merge pull request #280 from SeaifanAladdin/SCF_C
Browse files Browse the repository at this point in the history 
Merge: SCF Potential updated
  • Loading branch information
@jobovy
jobovy committed Jul 18, 2016
2 parents dd70069 + 56d3ec5 commit 39fe3f5
Show file tree
Hide file tree
Showing 5 changed files with 529 additions and 489 deletions.
15 changes: 7 additions & 8 deletions galpy/orbit_src/integratePlanarOrbit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -182,15 +182,14 @@ def _parse_pot(pot):
*p._Pot._glx[ii]**2.))
for ii in range(p._Pot._glorder)])
pot_args.extend([p._Pot._glw[ii] for ii in range(p._Pot._glorder)])
pot_args.extend([0.,0.,0.,0.,0.,0.])
pot_args.extend([p._RZPot._amp,p._RZPot._a])
elif isinstance(p,potential_src.planarPotential.planarPotentialFromRZPotential) \
and isinstance(p._RZPot,potential.SCFPotential):
pot_args.extend([0.,0.,0.,0.,0.,0.])
elif (isinstance(p,potential_src.planarPotential.planarPotentialFromFullPotential) or isinstance(p,potential_src.planarPotential.planarPotentialFromRZPotential)) \
and isinstance(p._Pot,potential.SCFPotential):
pot_type.append(24)
pot_args.extend([p._RZPot._a])
pot_args.extend(p._RZPot._Acos.shape)
pot_args.extend(p._RZPot._amp*p._RZPot._Acos.flatten(order='C'))
pot_args.extend(p._RZPot._amp*p._RZPot._Asin.flatten(order='C'))
pot_args.extend([p._Pot._a])
pot_args.extend(p._Pot._Acos.shape)
pot_args.extend(p._Pot._amp*p._Pot._Acos.flatten(order='C'))
pot_args.extend(p._Pot._amp*p._Pot._Asin.flatten(order='C'))
pot_args.extend([-1., 0, 0, 0, 0, 0, 0])
pot_type= nu.array(pot_type,dtype=nu.int32,order='C')
pot_args= nu.array(pot_args,dtype=nu.float64,order='C')
Expand Down
88 changes: 32 additions & 56 deletions galpy/potential_src/SCFPotential.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,6 @@
from scipy.special import gammaln

import hashlib
import inspect


class SCFPotential(Potential):
Expand All @@ -29,15 +28,15 @@ def __init__(self, amp=1., Acos=nu.array([[[1]]]),Asin=nu.array([[[0]]]), a = 1.
INPUT:
amp - amplitude to be applied to the potential (default: 1); can be a Quantity with units of mass density or Gxmass density
amp - amplitude to be applied to the potential (default: 1); can be a Quantity with units of mass density or Gxmass density
Acos - The real part of the expansion coefficent (NxLxL matrix)
Acos - The real part of the expansion coefficent (NxLxL matrix)
Asin - The imaginary part of the expansion coefficent (NxLxL matrix)
Asin - The imaginary part of the expansion coefficent (NxLxL matrix)
a - scale length (can be Quantity)
a - scale length (can be Quantity)
normalize - if True, normalize such that vc(1.,0.)=1., or, if given as a number, such that the force is this fraction of the force necessary to make vc(1.,0.)=1.
normalize - if True, normalize such that vc(1.,0.)=1., or, if given as a number, such that the force is this fraction of the force necessary to make vc(1.,0.)=1.
ro=, vo= distance and velocity scales for translation into internal units (default from configuration file)
Expand All @@ -48,7 +47,6 @@ def __init__(self, amp=1., Acos=nu.array([[[1]]]),Asin=nu.array([[[0]]]), a = 1.
HISTORY:
2016-05-13 - Written - Aladdin
"""
Potential.__init__(self,amp=amp/2.,ro=ro,vo=vo,amp_units='unitless')
if _APY_LOADED and isinstance(a,units.Quantity):
Expand Down Expand Up @@ -192,33 +190,6 @@ def _compute(self, funcTilde, R, z, phi):
func = func_tilde[:,:,None]*(Acos[:,:,:]*mcos + Asin[:,:,:]*msin)*PP[None,:,:]
return func

def _getShape(self, R, z, phi):
"""
NAME:
_getShape
PURPOSE:
get the shape of R, z and phi
INPUT:
R - Cylindrical Galactocentric radius
z - vertical height
phi - azimuth
OUTPUT:
R, z and phi as numpy arrays
HISTORY:
2016-06-02 - Written - Aladdin
"""
shape=None
##Determine which of these are arrays
if type(R).__name__ == nu.ndarray.__name__:
shape = R.shape
elif type(z).__name__ == nu.ndarray.__name__:
shape = z.shape
elif type(phi).__name__ == nu.ndarray.__name__:
shape = phi.shape
else:
shape = None

return shape

def _computeArray(self, funcTilde, R, z, phi):
"""
Expand All @@ -236,13 +207,14 @@ def _computeArray(self, funcTilde, R, z, phi):
HISTORY:
2016-06-02 - Written - Aladdin
"""
shape = self._getShape(R,z,phi)
if shape == None: return nu.sum(self._compute(funcTilde, R,z,phi))
R = nu.array(R); z = nu.array(z); phi = nu.array(phi);
shape = (R*z*phi).shape
if shape == (): return nu.sum(self._compute(funcTilde, R,z,phi))
R *= nu.ones(shape); z *= nu.ones(shape); phi *= nu.ones(shape);
func = nu.zeros(shape, float)


li = cartesian(shape)
li = _cartesian(shape)
for i in range(li.shape[0]):
func[li[i]] = nu.sum(self._compute(funcTilde, R[li[i]][0],z[li[i]][0],phi[li[i]][0]))
return func
Expand Down Expand Up @@ -323,7 +295,7 @@ def _computeforce(self,R,z,phi=0,t=0):
phi - azimuth
t - time
OUTPUT:
the derivative of phiTilde with respect to r
The force in the x direction
HISTORY:
2016-06-07 - Written - Aladdin
"""
Expand Down Expand Up @@ -372,12 +344,13 @@ def _computeforceArray(self,dr_dx, dtheta_dx, dphi_dx, R, z, phi):
phi - azimuth
t - time
OUTPUT:
density or potential evaluated at (R,z, phi)
The forces in the x direction
HISTORY:
2016-06-02 - Written - Aladdin
"""
shape = self._getShape(R,z,phi)
if shape == None:
R = nu.array(R); z = nu.array(z); phi = nu.array(phi);
shape = (R*z*phi).shape
if shape == ():
dPhi_dr,dPhi_dtheta,dPhi_dphi = \
self._computeforce(R,z,phi)
return dr_dx*dPhi_dr + dtheta_dx*dPhi_dtheta +dPhi_dphi*dphi_dx
Expand All @@ -386,7 +359,7 @@ def _computeforceArray(self,dr_dx, dtheta_dx, dphi_dx, R, z, phi):
force = nu.zeros(shape, float)


li = cartesian(shape)
li = _cartesian(shape)

for i in range(li.shape[0]):
dPhi_dr,dPhi_dtheta,dPhi_dphi = \
Expand Down Expand Up @@ -456,8 +429,11 @@ def _phiforce(self, R,z,phi=0,t=0):
dr_dphi = 0; dtheta_dphi = 0; dphi_dphi = 1
return self._computeforceArray(dr_dphi, dtheta_dphi, dphi_dphi, R,z,phi)

def OmegaP(self):
return 0


def xiToR(xi, a =1):
def _xiToR(xi, a =1):
return a*nu.divide((1. + xi),(1. - xi))


Expand Down Expand Up @@ -525,7 +501,7 @@ def scf_compute_coeffs_spherical(dens, N, a=1.):


def integrand(xi):
r = xiToR(xi, a)
r = _xiToR(xi, a)
R = r
param[0] = R
return a**3. * dens(*param)*(1 + xi)**2. * (1 - xi)**-3. * _C(xi, N, 1)[:,0]
Expand All @@ -534,7 +510,7 @@ def integrand(xi):
Asin = nu.zeros((N,1,1), float)

Ksample = [max(N + 1, 20)]
integrated = gaussianQuadrature(integrand, [[-1., 1.]], Ksample=Ksample)
integrated = _gaussianQuadrature(integrand, [[-1., 1.]], Ksample=Ksample)
n = nu.arange(0,N)
K = 16*nu.pi*(n + 3./2)/((n + 2)*(n + 1)*(1 + n*(n + 3.)/2.))
Acos[n,0,0] = 2*K*integrated
Expand Down Expand Up @@ -566,7 +542,7 @@ def scf_compute_coeffs_axi(dens, N, L, a=1.,radial_order=None, costheta_order=No
param = [0]*numOfParam;
def integrand(xi, costheta):
l = nu.arange(0, L)[nu.newaxis, :]
r = xiToR(xi,a)
r = _xiToR(xi,a)
R = r*nu.sqrt(1 - costheta**2.)
z = r*costheta
Legandre = lpmn(0,L-1,costheta)[0].T[nu.newaxis,:,0]
Expand All @@ -581,14 +557,14 @@ def integrand(xi, costheta):
Asin = nu.zeros((N,L,L), float)

##This should save us some computation time since we're only taking the double integral once, rather then L times
Ksample = [max(N + 3*L/2 + 1, 20) , max(L + 1,20) ]
Ksample = [max(N + 3*L//2 + 1, 20) , max(L + 1,20) ]
if radial_order != None:
Ksample[0] = radial_order
if costheta_order != None:
Ksample[1] = costheta_order


integrated = gaussianQuadrature(integrand, [[-1., 1.], [-1, 1]], Ksample = Ksample)*(2*nu.pi)
integrated = _gaussianQuadrature(integrand, [[-1, 1], [-1, 1]], Ksample = Ksample)*(2*nu.pi)
n = nu.arange(0,N)[:,nu.newaxis]
l = nu.arange(0,L)[nu.newaxis,:]
K = .5*n*(n + 4*l + 3) + (l + 1)*(2*l + 1)
Expand Down Expand Up @@ -623,7 +599,7 @@ def scf_compute_coeffs(dens, N, L, a=1., radial_order=None, costheta_order=None,
def integrand(xi, costheta, phi):
l = nu.arange(0, L)[nu.newaxis, :, nu.newaxis]
m = nu.arange(0, L)[nu.newaxis,nu.newaxis,:]
r = xiToR(xi, a)
r = _xiToR(xi, a)
R = r*nu.sqrt(1 - costheta**2.)
z = r*costheta
Legandre = lpmn(L - 1,L-1,costheta)[0].T[nu.newaxis,:,:]
Expand All @@ -638,15 +614,15 @@ def integrand(xi, costheta, phi):
Acos = nu.zeros((N,L,L), float)
Asin = nu.zeros((N,L,L), float)

##This should save us some computation time since we're only taking the Triple integral once, rather then L times

Ksample = [max(N + 3*L//2 + 1,20), max(L + 1,20 ), max(L + 1,20)]
if radial_order != None:
Ksample[0] = radial_order
if costheta_order != None:
Ksample[1] = costheta_order
if phi_order != None:
Ksample[2] = phi_order
integrated = gaussianQuadrature(integrand, [[-1., 1.], [-1., 1.], [0, 2*nu.pi]], Ksample = Ksample)
integrated = _gaussianQuadrature(integrand, [[-1., 1.], [-1., 1.], [0, 2*nu.pi]], Ksample = Ksample)
n = nu.arange(0,N)[:,nu.newaxis, nu.newaxis]
l = nu.arange(0,L)[nu.newaxis,:, nu.newaxis]
m = nu.arange(0,L)[nu.newaxis,nu.newaxis,:]
Expand All @@ -666,7 +642,7 @@ def integrand(xi, costheta, phi):

return Acos, Asin

def cartesian(arraySizes, out=None):
def _cartesian(arraySizes, out=None):
"""
NAME:
cartesian
Expand Down Expand Up @@ -696,12 +672,12 @@ def cartesian(arraySizes, out=None):
m = n / arrays[0].size
out[:,0] = nu.repeat(arrays[0], m)
if arrays[1:]:
cartesian(arraySizes[1:], out=out[0:m,1:])
for j in xrange(1, arrays[0].size):
_cartesian(arraySizes[1:], out=out[0:m,1:])
for j in range(1, arrays[0].size):
out[j*m:(j+1)*m,1:] = out[0:m,1:]
return out

def gaussianQuadrature(integrand, bounds, Ksample=[20], roundoff=0):
def _gaussianQuadrature(integrand, bounds, Ksample=[20], roundoff=0):
"""
NAME:
_gaussianQuadrature
Expand Down Expand Up @@ -742,7 +718,7 @@ def gaussianQuadrature(integrand, bounds, Ksample=[20], roundoff=0):


#gets all combinations of indices from each integrand
li = cartesian(Ksample)
li = _cartesian(Ksample)

##Performs the actual integration
for i in range(li.shape[0]):
Expand Down

0 comments on commit 39fe3f5

Please sign in to comment.