1 #include <complex>
2 #include <cmath>
3 #include <new>
4
5 #include "../utils.hpp"
6 #include "../numpypp/numpy.hpp"
7
8
9 namespace{
10
11 const char TypeErrorMsg[] =
12 "Type not understood. "
13 "This is caused by either a direct call to _zernike (which is dangerous: types are not checked!) or a bug in zernike.py.\n";
14
15 double _factorialtable[] = {
16 1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600
17 };
18
19 inline
fact(int n)20 double fact(int n) {
21 if (unsigned(n) < sizeof(_factorialtable)/sizeof(double)) return _factorialtable[n];
22 return double(n) * fact(n - 1);
23 }
24
py_znl(PyObject * self,PyObject * args)25 PyObject* py_znl(PyObject* self, PyObject* args) {
26 using numpy::ndarray_cast;
27 using std::pow;
28 using std::atan;
29 using std::atan2;
30 using std::polar;
31 using std::conj;
32 using std::complex;
33
34 const double pi = atan(1.0)*4;
35
36 PyArrayObject* Da;
37 PyArrayObject* Aa;
38 PyArrayObject* Pa;
39 int n;
40 int l;
41 if (!PyArg_ParseTuple(args,"OOOii", &Da, &Aa, &Pa, &n, &l)) return NULL;
42 if (!PyArray_Check(Da) || !PyArray_Check(Aa) || !PyArray_Check(Pa) ||
43 PyArray_TYPE(Da) != NPY_DOUBLE || PyArray_TYPE(Aa) != NPY_CDOUBLE || PyArray_TYPE(Pa) != NPY_DOUBLE) {
44 PyErr_SetString(PyExc_RuntimeError, TypeErrorMsg);
45 return NULL;
46 }
47 holdref Da_hr(Da);
48 holdref Aa_hr(Aa);
49 holdref Pa_hr(Pa);
50
51 double* D = ndarray_cast<double*>(Da);
52 complex<double>* A = ndarray_cast<complex<double>*>(Aa);
53 double* P = ndarray_cast<double*>(Pa);
54 const int Nelems = PyArray_SIZE(Da);
55 complex<double> v = 0.;
56 try {
57 gil_release nogil;
58 complex<double> Vnl = 0.0;
59 double * g_m = new double[ int( (n-l)/2 ) + 1];
60 for(int m = 0; m <= (n-l)/2; m++) {
61 double f = (m & 1) ? -1 : 1;
62 g_m[m] = f * fact(n-m) /
63 ( fact(m) * fact((n - 2*m + l) / 2) * fact((n - 2*m - l) / 2) );
64 }
65
66 for (int i = 0; i != Nelems; ++i) {
67 double d=D[i];
68 complex<double> a=A[i];
69 double p=P[i];
70 Vnl = 0.;
71 for(int m = 0; m <= (n-l)/2; m++) {
72 Vnl += g_m[m] * pow(d, double(n - 2*m)) * a;
73 }
74 v += p * conj(Vnl);
75 }
76 v *= (n+1)/pi;
77 delete [] g_m;
78 } catch (const std::bad_alloc&) {
79 PyErr_NoMemory();
80 return NULL;
81 }
82 return PyComplex_FromDoubles(v.real(), v.imag());
83 }
84 PyMethodDef methods[] = {
85 {"znl",(PyCFunction)py_znl, METH_VARARGS, NULL},
86 {NULL, NULL,0,NULL},
87 };
88
89 } // namespace
90 DECLARE_MODULE(_zernike)
91
92