def smooth_fft_vsini(dv,spec,sigma):
# The Fourier coordinate
ss = rfftfreq(len(spec), d=dv)
# Make the fourier space taper
ss[0] = 0.01 #junk so we don't get a divide by zero error
ub = 2. * np.pi * sigma * ss
sb = j1(ub) / ub - 3 * np.cos(ub) / (2 * ub ** 2) + 3. * np.sin(ub) / (2 * ub ** 3)
#set zeroth frequency to 1 separately (DC term)
sb[0] = 1.
# Fourier transform the spectrum
FF = np.fft.rfft(spec)
# Multiply in fourier space
FF_tap = FF * sb
# Fourier transform back
spec_conv = np.fft.irfft(FF_tap)
return spec_conv
python类j1()的实例源码
def jinc(r):
''' The Jinc function.
Args:
r (`number`): radial distance.
Returns:
`float`: the value of j1(x)/x for x != 0, 0.5 at 0.
'''
if r == 0:
return 0.5
else:
return j1(r) / r
def LambDipole(model, U=.01,R = 1.):
""" Generate Lamb's dipole vorticity field.
Parameters
-----------
U: float
Translation speed (dipole's strength).
R: float
Diple's radius.
Return
-------
q: array of floats
Vorticity (physical space).
"""
N = model.nx
x, y = model.x, model.y
x0,y0 = x[N//2,N//2],y[N//2,N//2]
r = np.sqrt( (x-x0)**2 + (y-y0)**2 )
s = np.zeros_like(r)
for i in range(N):
for j in range(N):
if r[i,j] == 0.:
s[i,j] = 0.
else:
s[i,j] = (y[i,j]-y0)/r[i,j]
lam = (3.8317)/R
C = -(2.*U*lam)/(special.j0(lam*R))
q = np.zeros_like(r)
q[r<=R] = C*special.j1(lam*r[r<=R])*s[r<=R]
return q
def airy(th, B, lam=None):
"""
Return the visibility value (unsquared), given
- th the angular diameter in radian
- B the baseline in m (alternatively in m/lambda)
- lam the wavelength (alternatively None if B is already given as m/lambda)
"""
if lam is None:
x = np.pi*th*B
else:
x = np.pi*th*B/lam
return 2*airyJ1(x)/x
def lanczos(n, fc=0.02):
"""
Compute the coefficients of a Lanczos window
Parameters
----------
n : int
Number of points in the output window, must be an odd integer.
fc : float
Cutoff frequency
Returns
-------
w : ndarray
The weights associated to the boxcar window
"""
if not isinstance(n, int):
try:
n = int(n)
except:
TypeError, "n must be an integer"
if not n % 2 == 1:
raise ValueError, "n must an odd integer"
dim = 1
if dim == 1:
k = np.arange(- n / 2 + 1, n / 2 + 1)
# k = np.arange(0, n + 1)
# w = (np.sin(2 * pi * fc * k) / (pi * k) *
# np.sin(pi * k / n) / (pi * k / n))
w = (np.sin(2. * np.pi * fc * k) / (np.pi * k) *
np.sin(np.pi * k / (n / 2.)) / (np.pi * k / (n / 2.)))
# Particular case where k=0
w[n / 2] = 2. * fc
elif dim == 2:
#TODO: Test this bidimensional window
fcx, fcy = fc
nx, ny = n
# Grid definition according to the number of weights
kx, ky = np.meshgrid(np.arange(-nx, nx + 1), np.arange(-ny, ny + 1), indexing='ij')
# Computation of the response weight on the grid
z = np.sqrt((fcx * kx) ** 2 + (fcy * ky) ** 2)
w_rect = 1 / z * fcx * fcy * spec.j1(2 * np.pi * z)
w = (w_rect * np.sin(np.pi * kx / nx) / (np.pi * kx / nx) *
np.sin(np.pi * ky / ny) / (np.pi * ky / ny))
# Particular case where z=0
w[nx, :] = (w_rect[nx, :] * 1. / (np.pi * ky[nx, :] / ny) *
np.sin(np.pi * ky[nx, :] / ny))
w[:, ny] = (w_rect[:, ny] * 1. / (np.pi * kx[:, ny] / nx) *
np.sin(np.pi * kx[:, ny] / nx))
w[nx, ny] = np.pi * fcx * fcy
return w
# Define a list of available windows
def quad(sPJ0r, sPJ0i, sPJ1r, sPJ1i, sPJ0br, sPJ0bi, ab, off, factAng, iinp):
"""Quadrature for Hankel transform.
This is the kernel of the QUAD method, used for the Hankel transforms
`hquad` and `hqwe` (where the integral is not suited for QWE).
"""
# Define the quadrature kernels
def quad0(klambd, sPJ, sPJb, koff, kang):
"""Quadrature for J0."""
tP0 = sPJ(np.log10(klambd)) + kang*sPJb(np.log10(klambd))
return tP0*special.j0(koff*klambd)
def quad1(klambd, sPJ, ab, koff, kang):
"""Quadrature for J1."""
tP1 = kang*sPJ(np.log10(klambd))
if ab in [11, 12, 21, 22, 14, 24, 15, 25]: # Because of J2
# J2(kr) = 2/(kr)*J1(kr) - J0(kr)
tP1 /= koff
return tP1*special.j1(koff*klambd)
# Carry out quadrature of J0
iargs = (sPJ0r, sPJ0br, off, factAng)
fr0 = integrate.quad(quad0, args=iargs, full_output=1, **iinp)
iargs = (sPJ0i, sPJ0bi, off, factAng)
fi0 = integrate.quad(quad0, args=iargs, full_output=1, **iinp)
# Carry out quadrature of J1
iargs = (sPJ1r, ab, off, factAng)
fr1 = integrate.quad(quad1, args=iargs, full_output=1, **iinp)
iargs = (sPJ1i, ab, off, factAng)
fi1 = integrate.quad(quad1, args=iargs, full_output=1, **iinp)
# If there is a fourth output from QUAD, it means it did not converge
if np.any(np.array([len(fr0), len(fi0), len(fr1), len(fi1)]) > 3):
conv = False
else:
conv = True
# Collect the results
return fr0[0] + fr1[0] + 1j*(fi0[0] + fi1[0]), conv
