ENH: scipy.special, Implemented Ellipsoidal harmonic function: the first kind, the second kind and normalization constant #3811
Conversation
… first kind, the second kind and normalization constant
|
@jennystone for next time: it would be better to have more smaller commits. Looks like you squashed everything together now in one large commit. |
| >>> from scipy.special import ellip_harm_2 | ||
| >>> w = ellip_normal(5,8,3,7) | ||
| >>> w | ||
| >>> 1723.38796997 |
There was a problem hiding this comment.
Output doesn't have >>>, should be only the number on this line.
| """ | ||
| return _ellip_harm(h2, k2, n, p, s, signm, signn) | ||
|
|
||
| def ellip_harm_2(h2, k2, n, p, s): |
There was a problem hiding this comment.
The definition should be changed so that it can evaluate several values at once:
# np.vectorize does not work on Cython functions on Numpy < 1.8, so a wrapper is needed
def _ellip_harm_2_vec(h2, k2, n, p, s):
return _ellipsoid(h2, k2, n, p, s)
_ellip_harm_2_vec = np.vectorize(_ellip_harm_2_vec, otypes='d')
def ellip_harm_2(h2, k2, n, p, s):
...
with _ellip_lock:
return _ellip_harm_2_vec(h2, k2, n, p, s)
The issue is that _ellipsoid is not an ufunc, and evaluates only one value at a time, so it needs to be wrapped.
| return 1/res | ||
| assert_almost_equal(recursion(100, 5, 2, 41, 5, 2, 15, 25), potential(100, 5, 2, 41, 5, 2, 15, 25)) | ||
| assert_almost_equal(recursion(120, sqrt(19), 2, 41, sqrt(17), 2, 15, 25), potential(120, sqrt(19), 2, 41, sqrt(17), 2, 15, 25)) | ||
| def test_ellip_norm(): |
There was a problem hiding this comment.
Style: empty line(s) before def ...:
|
I think we are getting quite close to finishing this up. The potential test proves that the implemented functions work as they should (at least within sensible parameter ranges). |
…es, Implementation of test with Electric Potential values, Improved implementation of ellipsoidal harmonics
| _F_integrand4_t = ctypes.CFUNCTYPE(ctypes.c_double, ctypes.c_double) | ||
| _F_integrand4_ctypes = ctypes.cast(<size_t>&_F_integrand4, _F_integrand4_t) | ||
|
|
||
| cdef double * calc_coefficients(double h2, double k2, int n, int p): |
There was a problem hiding this comment.
This routine is not necessary. You also don't need to store bufferp as a global variable as it can be a local variable in each of the routines.
| signn = signn/fabs(signn) | ||
|
|
||
| if n < 0: | ||
| sf_error.error("ellip_harm", sf_error.ARG, "invalid value for p") |
|
Merged in gh-4182. Thanks! |
This is an attempt to implement ellipsoidal harmonic functionality.
The functions implemented are
ellip_harm, ellip_harm_2, ellip_normal
The sources referred to are:
http://www.geology.osu.edu/~jekeli.1/OSUReports/reports/report_499.pdf
http://arxiv.org/abs/1204.0267