Generated by Cython 3.0.3
Yellow lines hint at Python interaction.
Click on a line that starts with a "+" to see the C code that Cython generated for it.
Raw output: linkingnumber_cython.c
+001: cimport numpy as np
__pyx_t_7 = __Pyx_PyDict_NewPresized(0); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 1, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_7); if (PyDict_SetItem(__pyx_d, __pyx_n_s_test, __pyx_t_7) < 0) __PYX_ERR(0, 1, __pyx_L1_error) __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;
002: import cython
003: # from libc.math cimport M_PI, asin,
004: from libc.math cimport sqrt, atan, isnan
005:
006: ctypedef np.float32_t DTYPE_t
007: from cpython.array cimport array, clone
008: from cython cimport view
009: from cython.view cimport array as cvarray
010:
011:
012: @cython.boundscheck(False)
013: @cython.wraparound(False)
014: @cython.nonecheck(False)
015: @cython.profile(True)
016: @cython.cdivision(True)
+017: cpdef double _eval_lk_cython(double[:, ::1] curve1, double[:, ::1] curve2):
static PyObject *__pyx_pw_4pylk_20linkingnumber_cython_1_eval_lk_cython(PyObject *__pyx_self,
#if CYTHON_METH_FASTCALL
PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds
#else
PyObject *__pyx_args, PyObject *__pyx_kwds
#endif
); /*proto*/
static double __pyx_f_4pylk_20linkingnumber_cython__eval_lk_cython(__Pyx_memviewslice __pyx_v_curve1, __Pyx_memviewslice __pyx_v_curve2, CYTHON_UNUSED int __pyx_skip_dispatch) {
int __pyx_v_N1;
int __pyx_v_N2;
int __pyx_v_i;
int __pyx_v_j;
double __pyx_v_lk;
double __pyx_v_dlk;
double __pyx_v_nv;
double __pyx_v_nw;
double __pyx_v_cosbeta;
double __pyx_v_sinbetasq;
double __pyx_v_a0;
double __pyx_v_a1;
double __pyx_v_a2;
double __pyx_v_val;
__Pyx_memviewslice __pyx_v_p1 = { 0, 0, { 0 }, { 0 }, { 0 } };
__Pyx_memviewslice __pyx_v_p2 = { 0, 0, { 0 }, { 0 }, { 0 } };
__Pyx_memviewslice __pyx_v_q1 = { 0, 0, { 0 }, { 0 }, { 0 } };
__Pyx_memviewslice __pyx_v_q2 = { 0, 0, { 0 }, { 0 }, { 0 } };
__Pyx_memviewslice __pyx_v_v = { 0, 0, { 0 }, { 0 }, { 0 } };
__Pyx_memviewslice __pyx_v_w = { 0, 0, { 0 }, { 0 }, { 0 } };
__Pyx_memviewslice __pyx_v_e1 = { 0, 0, { 0 }, { 0 }, { 0 } };
__Pyx_memviewslice __pyx_v_e2 = { 0, 0, { 0 }, { 0 }, { 0 } };
__Pyx_memviewslice __pyx_v_r12 = { 0, 0, { 0 }, { 0 }, { 0 } };
__Pyx_memviewslice __pyx_v_e2cosbetamine1 = { 0, 0, { 0 }, { 0 }, { 0 } };
__Pyx_memviewslice __pyx_v_e2mine1cosbeta = { 0, 0, { 0 }, { 0 }, { 0 } };
__Pyx_memviewslice __pyx_v_crosse1e2 = { 0, 0, { 0 }, { 0 }, { 0 } };
double __pyx_v_minoneover4pi;
double __pyx_v_a0sq;
double __pyx_v_num;
double __pyx_v_denom;
double __pyx_v_t1;
double __pyx_v_t2;
CYTHON_UNUSED double __pyx_v_a0sqcosbeta;
arrayobject *__pyx_v_templatemv = 0;
double __pyx_v_a0sqsinbetasq;
__Pyx_LocalBuf_ND __pyx_pybuffernd_templatemv;
__Pyx_Buffer __pyx_pybuffer_templatemv;
double __pyx_r;
__Pyx_TraceDeclarations
__Pyx_TraceFrameInit(__pyx_codeobj__11)
__Pyx_TraceCall("_eval_lk_cython", __pyx_f[0], 17, 0, __PYX_ERR(0, 17, __pyx_L1_error));
__pyx_pybuffer_templatemv.pybuffer.buf = NULL;
__pyx_pybuffer_templatemv.refcount = 0;
__pyx_pybuffernd_templatemv.data = NULL;
__pyx_pybuffernd_templatemv.rcbuffer = &__pyx_pybuffer_templatemv;
/* … */
/* function exit code */
__pyx_L1_error:;
__Pyx_XDECREF(__pyx_t_1);
__PYX_XCLEAR_MEMVIEW(&__pyx_t_6, 1);
{ PyObject *__pyx_type, *__pyx_value, *__pyx_tb;
__Pyx_PyThreadState_declare
__Pyx_PyThreadState_assign
__Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb);
__Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_templatemv.rcbuffer->pybuffer);
__Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
__Pyx_AddTraceback("pylk.linkingnumber_cython._eval_lk_cython", __pyx_clineno, __pyx_lineno, __pyx_filename);
__pyx_r = -1;
goto __pyx_L2;
__pyx_L0:;
__Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_templatemv.rcbuffer->pybuffer);
__pyx_L2:;
__PYX_XCLEAR_MEMVIEW(&__pyx_v_p1, 1);
__PYX_XCLEAR_MEMVIEW(&__pyx_v_p2, 1);
__PYX_XCLEAR_MEMVIEW(&__pyx_v_q1, 1);
__PYX_XCLEAR_MEMVIEW(&__pyx_v_q2, 1);
__PYX_XCLEAR_MEMVIEW(&__pyx_v_v, 1);
__PYX_XCLEAR_MEMVIEW(&__pyx_v_w, 1);
__PYX_XCLEAR_MEMVIEW(&__pyx_v_e1, 1);
__PYX_XCLEAR_MEMVIEW(&__pyx_v_e2, 1);
__PYX_XCLEAR_MEMVIEW(&__pyx_v_r12, 1);
__PYX_XCLEAR_MEMVIEW(&__pyx_v_e2cosbetamine1, 1);
__PYX_XCLEAR_MEMVIEW(&__pyx_v_e2mine1cosbeta, 1);
__PYX_XCLEAR_MEMVIEW(&__pyx_v_crosse1e2, 1);
__Pyx_XDECREF((PyObject *)__pyx_v_templatemv);
__Pyx_TraceReturn(Py_None, 0);
__Pyx_RefNannyFinishContext();
return __pyx_r;
}
/* Python wrapper */
static PyObject *__pyx_pw_4pylk_20linkingnumber_cython_1_eval_lk_cython(PyObject *__pyx_self,
#if CYTHON_METH_FASTCALL
PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds
#else
PyObject *__pyx_args, PyObject *__pyx_kwds
#endif
); /*proto*/
static PyMethodDef __pyx_mdef_4pylk_20linkingnumber_cython_1_eval_lk_cython = {"_eval_lk_cython", (PyCFunction)(void*)(__Pyx_PyCFunction_FastCallWithKeywords)__pyx_pw_4pylk_20linkingnumber_cython_1_eval_lk_cython, __Pyx_METH_FASTCALL|METH_KEYWORDS, 0};
static PyObject *__pyx_pw_4pylk_20linkingnumber_cython_1_eval_lk_cython(PyObject *__pyx_self,
#if CYTHON_METH_FASTCALL
PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds
#else
PyObject *__pyx_args, PyObject *__pyx_kwds
#endif
) {
__Pyx_memviewslice __pyx_v_curve1 = { 0, 0, { 0 }, { 0 }, { 0 } };
__Pyx_memviewslice __pyx_v_curve2 = { 0, 0, { 0 }, { 0 }, { 0 } };
#if !CYTHON_METH_FASTCALL
CYTHON_UNUSED Py_ssize_t __pyx_nargs;
#endif
CYTHON_UNUSED PyObject *const *__pyx_kwvalues;
PyObject *__pyx_r = 0;
__Pyx_RefNannyDeclarations
__Pyx_RefNannySetupContext("_eval_lk_cython (wrapper)", 0);
#if !CYTHON_METH_FASTCALL
#if CYTHON_ASSUME_SAFE_MACROS
__pyx_nargs = PyTuple_GET_SIZE(__pyx_args);
#else
__pyx_nargs = PyTuple_Size(__pyx_args); if (unlikely(__pyx_nargs < 0)) return NULL;
#endif
#endif
__pyx_kwvalues = __Pyx_KwValues_FASTCALL(__pyx_args, __pyx_nargs);
{
PyObject **__pyx_pyargnames[] = {&__pyx_n_s_curve1,&__pyx_n_s_curve2,0};
PyObject* values[2] = {0,0};
if (__pyx_kwds) {
Py_ssize_t kw_args;
switch (__pyx_nargs) {
case 2: values[1] = __Pyx_Arg_FASTCALL(__pyx_args, 1);
CYTHON_FALLTHROUGH;
case 1: values[0] = __Pyx_Arg_FASTCALL(__pyx_args, 0);
CYTHON_FALLTHROUGH;
case 0: break;
default: goto __pyx_L5_argtuple_error;
}
kw_args = __Pyx_NumKwargs_FASTCALL(__pyx_kwds);
switch (__pyx_nargs) {
case 0:
if (likely((values[0] = __Pyx_GetKwValue_FASTCALL(__pyx_kwds, __pyx_kwvalues, __pyx_n_s_curve1)) != 0)) {
(void)__Pyx_Arg_NewRef_FASTCALL(values[0]);
kw_args--;
}
else if (unlikely(PyErr_Occurred())) __PYX_ERR(0, 17, __pyx_L3_error)
else goto __pyx_L5_argtuple_error;
CYTHON_FALLTHROUGH;
case 1:
if (likely((values[1] = __Pyx_GetKwValue_FASTCALL(__pyx_kwds, __pyx_kwvalues, __pyx_n_s_curve2)) != 0)) {
(void)__Pyx_Arg_NewRef_FASTCALL(values[1]);
kw_args--;
}
else if (unlikely(PyErr_Occurred())) __PYX_ERR(0, 17, __pyx_L3_error)
else {
__Pyx_RaiseArgtupleInvalid("_eval_lk_cython", 1, 2, 2, 1); __PYX_ERR(0, 17, __pyx_L3_error)
}
}
if (unlikely(kw_args > 0)) {
const Py_ssize_t kwd_pos_args = __pyx_nargs;
if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_kwvalues, __pyx_pyargnames, 0, values + 0, kwd_pos_args, "_eval_lk_cython") < 0)) __PYX_ERR(0, 17, __pyx_L3_error)
}
} else if (unlikely(__pyx_nargs != 2)) {
goto __pyx_L5_argtuple_error;
} else {
values[0] = __Pyx_Arg_FASTCALL(__pyx_args, 0);
values[1] = __Pyx_Arg_FASTCALL(__pyx_args, 1);
}
__pyx_v_curve1 = __Pyx_PyObject_to_MemoryviewSlice_d_dc_double(values[0], PyBUF_WRITABLE); if (unlikely(!__pyx_v_curve1.memview)) __PYX_ERR(0, 17, __pyx_L3_error)
__pyx_v_curve2 = __Pyx_PyObject_to_MemoryviewSlice_d_dc_double(values[1], PyBUF_WRITABLE); if (unlikely(!__pyx_v_curve2.memview)) __PYX_ERR(0, 17, __pyx_L3_error)
}
goto __pyx_L6_skip;
__pyx_L5_argtuple_error:;
__Pyx_RaiseArgtupleInvalid("_eval_lk_cython", 1, 2, 2, __pyx_nargs); __PYX_ERR(0, 17, __pyx_L3_error)
__pyx_L6_skip:;
goto __pyx_L4_argument_unpacking_done;
__pyx_L3_error:;
{
Py_ssize_t __pyx_temp;
for (__pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) {
__Pyx_Arg_XDECREF_FASTCALL(values[__pyx_temp]);
}
}
__PYX_XCLEAR_MEMVIEW(&__pyx_v_curve1, 1);
__PYX_XCLEAR_MEMVIEW(&__pyx_v_curve2, 1);
__Pyx_AddTraceback("pylk.linkingnumber_cython._eval_lk_cython", __pyx_clineno, __pyx_lineno, __pyx_filename);
__Pyx_RefNannyFinishContext();
return NULL;
__pyx_L4_argument_unpacking_done:;
__pyx_r = __pyx_pf_4pylk_20linkingnumber_cython__eval_lk_cython(__pyx_self, __pyx_v_curve1, __pyx_v_curve2);
int __pyx_lineno = 0;
const char *__pyx_filename = NULL;
int __pyx_clineno = 0;
/* function exit code */
__PYX_XCLEAR_MEMVIEW(&__pyx_v_curve1, 1);
__PYX_XCLEAR_MEMVIEW(&__pyx_v_curve2, 1);
{
Py_ssize_t __pyx_temp;
for (__pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) {
__Pyx_Arg_XDECREF_FASTCALL(values[__pyx_temp]);
}
}
__Pyx_RefNannyFinishContext();
return __pyx_r;
}
static PyObject *__pyx_pf_4pylk_20linkingnumber_cython__eval_lk_cython(CYTHON_UNUSED PyObject *__pyx_self, __Pyx_memviewslice __pyx_v_curve1, __Pyx_memviewslice __pyx_v_curve2) {
PyObject *__pyx_r = NULL;
__Pyx_TraceDeclarations
__Pyx_TraceFrameInit(__pyx_codeobj__11)
__Pyx_TraceCall("_eval_lk_cython (wrapper)", __pyx_f[0], 17, 0, __PYX_ERR(0, 17, __pyx_L1_error));
__Pyx_XDECREF(__pyx_r);
if (unlikely(!__pyx_v_curve1.memview)) { __Pyx_RaiseUnboundLocalError("curve1"); __PYX_ERR(0, 17, __pyx_L1_error) }
if (unlikely(!__pyx_v_curve2.memview)) { __Pyx_RaiseUnboundLocalError("curve2"); __PYX_ERR(0, 17, __pyx_L1_error) }
__pyx_t_1 = __pyx_f_4pylk_20linkingnumber_cython__eval_lk_cython(__pyx_v_curve1, __pyx_v_curve2, 0); if (unlikely(__pyx_t_1 == ((double)-1) && PyErr_Occurred())) __PYX_ERR(0, 17, __pyx_L1_error)
__pyx_t_2 = PyFloat_FromDouble(__pyx_t_1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 17, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__pyx_r = __pyx_t_2;
__pyx_t_2 = 0;
goto __pyx_L0;
/* function exit code */
__pyx_L1_error:;
__Pyx_XDECREF(__pyx_t_2);
__Pyx_AddTraceback("pylk.linkingnumber_cython._eval_lk_cython", __pyx_clineno, __pyx_lineno, __pyx_filename);
__pyx_r = NULL;
__pyx_L0:;
__Pyx_XGIVEREF(__pyx_r);
__Pyx_TraceReturn(__pyx_r, 0);
__Pyx_RefNannyFinishContext();
return __pyx_r;
}
/* … */
__pyx_tuple__24 = PyTuple_Pack(2, __pyx_n_s_curve1, __pyx_n_s_curve2); if (unlikely(!__pyx_tuple__24)) __PYX_ERR(0, 17, __pyx_L1_error)
__Pyx_GOTREF(__pyx_tuple__24);
__Pyx_GIVEREF(__pyx_tuple__24);
/* … */
__pyx_t_7 = __Pyx_CyFunction_New(&__pyx_mdef_4pylk_20linkingnumber_cython_1_eval_lk_cython, 0, __pyx_n_s_eval_lk_cython, NULL, __pyx_n_s_pylk_linkingnumber_cython, __pyx_d, ((PyObject *)__pyx_codeobj__11)); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 17, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_7);
if (PyDict_SetItem(__pyx_d, __pyx_n_s_eval_lk_cython, __pyx_t_7) < 0) __PYX_ERR(0, 17, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;
018: cdef int N1,N2,i,j
019: cdef double lk,dlk,nv,nw,cosbeta,sinbetasq,a0,a1,a2,val
020: cdef double[::1] p1,p2,q1,q2,v,w,e1,e2,r12
021: cdef double[::1] e2cosbetamine1,e2mine1cosbeta,crosse1e2
022:
023: cdef double minoneover4pi, a0sq, num, denom
024: cdef double t1,t2,a0sqcosbeta
025:
026: cdef array[double] templatemv
+027: templatemv = array('d')
__pyx_t_1 = __Pyx_PyObject_Call(((PyObject *)__pyx_ptype_7cpython_5array_array), __pyx_tuple__12, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 27, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); { __Pyx_BufFmt_StackElem __pyx_stack[1]; __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_templatemv.rcbuffer->pybuffer); __pyx_t_2 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_templatemv.rcbuffer->pybuffer, (PyObject*)((arrayobject *)__pyx_t_1), &__Pyx_TypeInfo_double, PyBUF_FORMAT| PyBUF_INDIRECT, 1, 0, __pyx_stack); if (unlikely(__pyx_t_2 < 0)) { PyErr_Fetch(&__pyx_t_3, &__pyx_t_4, &__pyx_t_5); if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_templatemv.rcbuffer->pybuffer, (PyObject*)__pyx_v_templatemv, &__Pyx_TypeInfo_double, PyBUF_FORMAT| PyBUF_INDIRECT, 1, 0, __pyx_stack) == -1)) { Py_XDECREF(__pyx_t_3); Py_XDECREF(__pyx_t_4); Py_XDECREF(__pyx_t_5); __Pyx_RaiseBufferFallbackError(); } else { PyErr_Restore(__pyx_t_3, __pyx_t_4, __pyx_t_5); } __pyx_t_3 = __pyx_t_4 = __pyx_t_5 = 0; } __pyx_pybuffernd_templatemv.diminfo[0].strides = __pyx_pybuffernd_templatemv.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_templatemv.diminfo[0].shape = __pyx_pybuffernd_templatemv.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_templatemv.diminfo[0].suboffsets = __pyx_pybuffernd_templatemv.rcbuffer->pybuffer.suboffsets[0]; if (unlikely((__pyx_t_2 < 0))) __PYX_ERR(0, 27, __pyx_L1_error) } __pyx_v_templatemv = ((arrayobject *)__pyx_t_1); __pyx_t_1 = 0; /* … */ __pyx_tuple__12 = PyTuple_Pack(1, __pyx_n_s_d); if (unlikely(!__pyx_tuple__12)) __PYX_ERR(0, 27, __pyx_L1_error) __Pyx_GOTREF(__pyx_tuple__12); __Pyx_GIVEREF(__pyx_tuple__12);
028:
+029: p1 = clone(templatemv,3,False)
__pyx_t_1 = ((PyObject *)__pyx_f_7cpython_5array_clone(((arrayobject *)__pyx_v_templatemv), 3, 0)); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 29, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); __pyx_t_6 = __Pyx_PyObject_to_MemoryviewSlice_dc_double(__pyx_t_1, PyBUF_WRITABLE); if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 29, __pyx_L1_error) __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; __pyx_v_p1 = __pyx_t_6; __pyx_t_6.memview = NULL; __pyx_t_6.data = NULL;
+030: p2 = clone(templatemv,3,False)
__pyx_t_1 = ((PyObject *)__pyx_f_7cpython_5array_clone(((arrayobject *)__pyx_v_templatemv), 3, 0)); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 30, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); __pyx_t_6 = __Pyx_PyObject_to_MemoryviewSlice_dc_double(__pyx_t_1, PyBUF_WRITABLE); if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 30, __pyx_L1_error) __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; __pyx_v_p2 = __pyx_t_6; __pyx_t_6.memview = NULL; __pyx_t_6.data = NULL;
+031: q1 = clone(templatemv,3,False)
__pyx_t_1 = ((PyObject *)__pyx_f_7cpython_5array_clone(((arrayobject *)__pyx_v_templatemv), 3, 0)); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 31, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); __pyx_t_6 = __Pyx_PyObject_to_MemoryviewSlice_dc_double(__pyx_t_1, PyBUF_WRITABLE); if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 31, __pyx_L1_error) __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; __pyx_v_q1 = __pyx_t_6; __pyx_t_6.memview = NULL; __pyx_t_6.data = NULL;
+032: q2 = clone(templatemv,3,False)
__pyx_t_1 = ((PyObject *)__pyx_f_7cpython_5array_clone(((arrayobject *)__pyx_v_templatemv), 3, 0)); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 32, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); __pyx_t_6 = __Pyx_PyObject_to_MemoryviewSlice_dc_double(__pyx_t_1, PyBUF_WRITABLE); if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 32, __pyx_L1_error) __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; __pyx_v_q2 = __pyx_t_6; __pyx_t_6.memview = NULL; __pyx_t_6.data = NULL;
+033: v = clone(templatemv,3,False)
__pyx_t_1 = ((PyObject *)__pyx_f_7cpython_5array_clone(((arrayobject *)__pyx_v_templatemv), 3, 0)); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 33, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); __pyx_t_6 = __Pyx_PyObject_to_MemoryviewSlice_dc_double(__pyx_t_1, PyBUF_WRITABLE); if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 33, __pyx_L1_error) __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; __pyx_v_v = __pyx_t_6; __pyx_t_6.memview = NULL; __pyx_t_6.data = NULL;
+034: w = clone(templatemv,3,False)
__pyx_t_1 = ((PyObject *)__pyx_f_7cpython_5array_clone(((arrayobject *)__pyx_v_templatemv), 3, 0)); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 34, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); __pyx_t_6 = __Pyx_PyObject_to_MemoryviewSlice_dc_double(__pyx_t_1, PyBUF_WRITABLE); if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 34, __pyx_L1_error) __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; __pyx_v_w = __pyx_t_6; __pyx_t_6.memview = NULL; __pyx_t_6.data = NULL;
+035: e1 = clone(templatemv,3,False)
__pyx_t_1 = ((PyObject *)__pyx_f_7cpython_5array_clone(((arrayobject *)__pyx_v_templatemv), 3, 0)); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 35, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); __pyx_t_6 = __Pyx_PyObject_to_MemoryviewSlice_dc_double(__pyx_t_1, PyBUF_WRITABLE); if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 35, __pyx_L1_error) __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; __pyx_v_e1 = __pyx_t_6; __pyx_t_6.memview = NULL; __pyx_t_6.data = NULL;
+036: e2 = clone(templatemv,3,False)
__pyx_t_1 = ((PyObject *)__pyx_f_7cpython_5array_clone(((arrayobject *)__pyx_v_templatemv), 3, 0)); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 36, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); __pyx_t_6 = __Pyx_PyObject_to_MemoryviewSlice_dc_double(__pyx_t_1, PyBUF_WRITABLE); if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 36, __pyx_L1_error) __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; __pyx_v_e2 = __pyx_t_6; __pyx_t_6.memview = NULL; __pyx_t_6.data = NULL;
+037: r12 = clone(templatemv,3,False)
__pyx_t_1 = ((PyObject *)__pyx_f_7cpython_5array_clone(((arrayobject *)__pyx_v_templatemv), 3, 0)); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 37, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); __pyx_t_6 = __Pyx_PyObject_to_MemoryviewSlice_dc_double(__pyx_t_1, PyBUF_WRITABLE); if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 37, __pyx_L1_error) __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; __pyx_v_r12 = __pyx_t_6; __pyx_t_6.memview = NULL; __pyx_t_6.data = NULL;
+038: e2cosbetamine1 = clone(templatemv,3,False)
__pyx_t_1 = ((PyObject *)__pyx_f_7cpython_5array_clone(((arrayobject *)__pyx_v_templatemv), 3, 0)); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 38, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); __pyx_t_6 = __Pyx_PyObject_to_MemoryviewSlice_dc_double(__pyx_t_1, PyBUF_WRITABLE); if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 38, __pyx_L1_error) __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; __pyx_v_e2cosbetamine1 = __pyx_t_6; __pyx_t_6.memview = NULL; __pyx_t_6.data = NULL;
+039: e2mine1cosbeta = clone(templatemv,3,False)
__pyx_t_1 = ((PyObject *)__pyx_f_7cpython_5array_clone(((arrayobject *)__pyx_v_templatemv), 3, 0)); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 39, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); __pyx_t_6 = __Pyx_PyObject_to_MemoryviewSlice_dc_double(__pyx_t_1, PyBUF_WRITABLE); if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 39, __pyx_L1_error) __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; __pyx_v_e2mine1cosbeta = __pyx_t_6; __pyx_t_6.memview = NULL; __pyx_t_6.data = NULL;
+040: crosse1e2 = clone(templatemv,3,False)
__pyx_t_1 = ((PyObject *)__pyx_f_7cpython_5array_clone(((arrayobject *)__pyx_v_templatemv), 3, 0)); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 40, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); __pyx_t_6 = __Pyx_PyObject_to_MemoryviewSlice_dc_double(__pyx_t_1, PyBUF_WRITABLE); if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 40, __pyx_L1_error) __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; __pyx_v_crosse1e2 = __pyx_t_6; __pyx_t_6.memview = NULL; __pyx_t_6.data = NULL;
041:
+042: N1 = len(curve1)
__pyx_t_7 = __Pyx_MemoryView_Len(__pyx_v_curve1);
__pyx_v_N1 = __pyx_t_7;
+043: N2 = len(curve2)
__pyx_t_7 = __Pyx_MemoryView_Len(__pyx_v_curve2);
__pyx_v_N2 = __pyx_t_7;
+044: minoneover4pi = -0.07957747154594767
__pyx_v_minoneover4pi = -0.07957747154594767;
045: # N1 = curve1.shape[0]
046: # N2 = curve2.shape[0]
047:
+048: lk = 0
__pyx_v_lk = 0.0;
+049: for i in range(N1):
__pyx_t_2 = __pyx_v_N1;
__pyx_t_8 = __pyx_t_2;
for (__pyx_t_9 = 0; __pyx_t_9 < __pyx_t_8; __pyx_t_9+=1) {
__pyx_v_i = __pyx_t_9;
+050: for j in range(N2):
__pyx_t_10 = __pyx_v_N2;
__pyx_t_11 = __pyx_t_10;
for (__pyx_t_12 = 0; __pyx_t_12 < __pyx_t_11; __pyx_t_12+=1) {
__pyx_v_j = __pyx_t_12;
+051: p1 = curve1[i]
__pyx_t_6.data = __pyx_v_curve1.data;
__pyx_t_6.memview = __pyx_v_curve1.memview;
__PYX_INC_MEMVIEW(&__pyx_t_6, 1);
{
Py_ssize_t __pyx_tmp_idx = __pyx_v_i;
Py_ssize_t __pyx_tmp_stride = __pyx_v_curve1.strides[0];
__pyx_t_6.data += __pyx_tmp_idx * __pyx_tmp_stride;
}
__pyx_t_6.shape[0] = __pyx_v_curve1.shape[1];
__pyx_t_6.strides[0] = __pyx_v_curve1.strides[1];
__pyx_t_6.suboffsets[0] = -1;
__PYX_XCLEAR_MEMVIEW(&__pyx_v_p1, 1);
__pyx_v_p1 = __pyx_t_6;
__pyx_t_6.memview = NULL;
__pyx_t_6.data = NULL;
+052: p2 = curve1[(i+1)%N1]
__pyx_t_6.data = __pyx_v_curve1.data;
__pyx_t_6.memview = __pyx_v_curve1.memview;
__PYX_INC_MEMVIEW(&__pyx_t_6, 1);
{
Py_ssize_t __pyx_tmp_idx = ((__pyx_v_i + 1) % __pyx_v_N1);
Py_ssize_t __pyx_tmp_stride = __pyx_v_curve1.strides[0];
__pyx_t_6.data += __pyx_tmp_idx * __pyx_tmp_stride;
}
__pyx_t_6.shape[0] = __pyx_v_curve1.shape[1];
__pyx_t_6.strides[0] = __pyx_v_curve1.strides[1];
__pyx_t_6.suboffsets[0] = -1;
__PYX_XCLEAR_MEMVIEW(&__pyx_v_p2, 1);
__pyx_v_p2 = __pyx_t_6;
__pyx_t_6.memview = NULL;
__pyx_t_6.data = NULL;
+053: q1 = curve2[j]
__pyx_t_6.data = __pyx_v_curve2.data;
__pyx_t_6.memview = __pyx_v_curve2.memview;
__PYX_INC_MEMVIEW(&__pyx_t_6, 1);
{
Py_ssize_t __pyx_tmp_idx = __pyx_v_j;
Py_ssize_t __pyx_tmp_stride = __pyx_v_curve2.strides[0];
__pyx_t_6.data += __pyx_tmp_idx * __pyx_tmp_stride;
}
__pyx_t_6.shape[0] = __pyx_v_curve2.shape[1];
__pyx_t_6.strides[0] = __pyx_v_curve2.strides[1];
__pyx_t_6.suboffsets[0] = -1;
__PYX_XCLEAR_MEMVIEW(&__pyx_v_q1, 1);
__pyx_v_q1 = __pyx_t_6;
__pyx_t_6.memview = NULL;
__pyx_t_6.data = NULL;
+054: q2 = curve2[(j+1)%N2]
__pyx_t_6.data = __pyx_v_curve2.data;
__pyx_t_6.memview = __pyx_v_curve2.memview;
__PYX_INC_MEMVIEW(&__pyx_t_6, 1);
{
Py_ssize_t __pyx_tmp_idx = ((__pyx_v_j + 1) % __pyx_v_N2);
Py_ssize_t __pyx_tmp_stride = __pyx_v_curve2.strides[0];
__pyx_t_6.data += __pyx_tmp_idx * __pyx_tmp_stride;
}
__pyx_t_6.shape[0] = __pyx_v_curve2.shape[1];
__pyx_t_6.strides[0] = __pyx_v_curve2.strides[1];
__pyx_t_6.suboffsets[0] = -1;
__PYX_XCLEAR_MEMVIEW(&__pyx_v_q2, 1);
__pyx_v_q2 = __pyx_t_6;
__pyx_t_6.memview = NULL;
__pyx_t_6.data = NULL;
055:
+056: v[0] = p2[0]-p1[0]
__pyx_t_13 = 0;
__pyx_t_14 = 0;
__pyx_t_15 = 0;
*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_v.data) + __pyx_t_15)) )) = ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_p2.data) + __pyx_t_13)) ))) - (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_p1.data) + __pyx_t_14)) ))));
+057: v[1] = p2[1]-p1[1]
__pyx_t_14 = 1;
__pyx_t_13 = 1;
__pyx_t_15 = 1;
*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_v.data) + __pyx_t_15)) )) = ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_p2.data) + __pyx_t_14)) ))) - (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_p1.data) + __pyx_t_13)) ))));
+058: v[2] = p2[2]-p1[2]
__pyx_t_13 = 2;
__pyx_t_14 = 2;
__pyx_t_15 = 2;
*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_v.data) + __pyx_t_15)) )) = ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_p2.data) + __pyx_t_13)) ))) - (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_p1.data) + __pyx_t_14)) ))));
059:
+060: w[0] = q2[0]-q1[0]
__pyx_t_14 = 0;
__pyx_t_13 = 0;
__pyx_t_15 = 0;
*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_w.data) + __pyx_t_15)) )) = ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_q2.data) + __pyx_t_14)) ))) - (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_q1.data) + __pyx_t_13)) ))));
+061: w[1] = q2[1]-q1[1]
__pyx_t_13 = 1;
__pyx_t_14 = 1;
__pyx_t_15 = 1;
*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_w.data) + __pyx_t_15)) )) = ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_q2.data) + __pyx_t_13)) ))) - (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_q1.data) + __pyx_t_14)) ))));
+062: w[2] = q2[2]-q1[2]
__pyx_t_14 = 2;
__pyx_t_13 = 2;
__pyx_t_15 = 2;
*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_w.data) + __pyx_t_15)) )) = ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_q2.data) + __pyx_t_14)) ))) - (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_q1.data) + __pyx_t_13)) ))));
063:
+064: nv = sqrt(v[0]**2+v[1]**2+v[2]**2)
__pyx_t_13 = 0;
__pyx_t_14 = 1;
__pyx_t_15 = 2;
__pyx_v_nv = sqrt(((pow((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_v.data) + __pyx_t_13)) ))), 2.0) + pow((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_v.data) + __pyx_t_14)) ))), 2.0)) + pow((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_v.data) + __pyx_t_15)) ))), 2.0)));
+065: nw = sqrt(w[0]**2+w[1]**2+w[2]**2)
__pyx_t_15 = 0;
__pyx_t_14 = 1;
__pyx_t_13 = 2;
__pyx_v_nw = sqrt(((pow((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_w.data) + __pyx_t_15)) ))), 2.0) + pow((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_w.data) + __pyx_t_14)) ))), 2.0)) + pow((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_w.data) + __pyx_t_13)) ))), 2.0)));
066:
+067: e1[0] = v[0]/nv
__pyx_t_13 = 0;
__pyx_t_14 = 0;
*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e1.data) + __pyx_t_14)) )) = ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_v.data) + __pyx_t_13)) ))) / __pyx_v_nv);
+068: e1[1] = v[1]/nv
__pyx_t_13 = 1;
__pyx_t_14 = 1;
*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e1.data) + __pyx_t_14)) )) = ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_v.data) + __pyx_t_13)) ))) / __pyx_v_nv);
+069: e1[2] = v[2]/nv
__pyx_t_13 = 2;
__pyx_t_14 = 2;
*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e1.data) + __pyx_t_14)) )) = ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_v.data) + __pyx_t_13)) ))) / __pyx_v_nv);
+070: e2[0] = w[0]/nw
__pyx_t_13 = 0;
__pyx_t_14 = 0;
*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2.data) + __pyx_t_14)) )) = ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_w.data) + __pyx_t_13)) ))) / __pyx_v_nw);
+071: e2[1] = w[1]/nw
__pyx_t_13 = 1;
__pyx_t_14 = 1;
*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2.data) + __pyx_t_14)) )) = ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_w.data) + __pyx_t_13)) ))) / __pyx_v_nw);
+072: e2[2] = w[2]/nw
__pyx_t_13 = 2;
__pyx_t_14 = 2;
*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2.data) + __pyx_t_14)) )) = ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_w.data) + __pyx_t_13)) ))) / __pyx_v_nw);
073:
+074: r12[0] = q1[0]-p1[0]
__pyx_t_13 = 0;
__pyx_t_14 = 0;
__pyx_t_15 = 0;
*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_r12.data) + __pyx_t_15)) )) = ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_q1.data) + __pyx_t_13)) ))) - (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_p1.data) + __pyx_t_14)) ))));
+075: r12[1] = q1[1]-p1[1]
__pyx_t_14 = 1;
__pyx_t_13 = 1;
__pyx_t_15 = 1;
*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_r12.data) + __pyx_t_15)) )) = ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_q1.data) + __pyx_t_14)) ))) - (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_p1.data) + __pyx_t_13)) ))));
+076: r12[2] = q1[2]-p1[2]
__pyx_t_13 = 2;
__pyx_t_14 = 2;
__pyx_t_15 = 2;
*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_r12.data) + __pyx_t_15)) )) = ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_q1.data) + __pyx_t_13)) ))) - (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_p1.data) + __pyx_t_14)) ))));
077:
+078: cosbeta = e1[0]*e2[0] + e1[1]*e2[1] + e1[2]*e2[2]
__pyx_t_14 = 0;
__pyx_t_13 = 0;
__pyx_t_15 = 1;
__pyx_t_16 = 1;
__pyx_t_17 = 2;
__pyx_t_18 = 2;
__pyx_v_cosbeta = ((((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e1.data) + __pyx_t_14)) ))) * (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2.data) + __pyx_t_13)) )))) + ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e1.data) + __pyx_t_15)) ))) * (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2.data) + __pyx_t_16)) ))))) + ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e1.data) + __pyx_t_17)) ))) * (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2.data) + __pyx_t_18)) )))));
+079: sinbetasq = 1 - cosbeta**2
__pyx_v_sinbetasq = (1.0 - pow(__pyx_v_cosbeta, 2.0));
080:
+081: e2cosbetamine1[0] = e2[0]*cosbeta - e1[0]
__pyx_t_18 = 0;
__pyx_t_17 = 0;
__pyx_t_16 = 0;
*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2cosbetamine1.data) + __pyx_t_16)) )) = (((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2.data) + __pyx_t_18)) ))) * __pyx_v_cosbeta) - (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e1.data) + __pyx_t_17)) ))));
+082: e2cosbetamine1[1] = e2[1]*cosbeta - e1[1]
__pyx_t_17 = 1;
__pyx_t_18 = 1;
__pyx_t_16 = 1;
*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2cosbetamine1.data) + __pyx_t_16)) )) = (((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2.data) + __pyx_t_17)) ))) * __pyx_v_cosbeta) - (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e1.data) + __pyx_t_18)) ))));
+083: e2cosbetamine1[2] = e2[2]*cosbeta - e1[2]
__pyx_t_18 = 2;
__pyx_t_17 = 2;
__pyx_t_16 = 2;
*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2cosbetamine1.data) + __pyx_t_16)) )) = (((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2.data) + __pyx_t_18)) ))) * __pyx_v_cosbeta) - (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e1.data) + __pyx_t_17)) ))));
084:
+085: e2mine1cosbeta[0] = e2[0] - e1[0]*cosbeta
__pyx_t_17 = 0;
__pyx_t_18 = 0;
__pyx_t_16 = 0;
*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2mine1cosbeta.data) + __pyx_t_16)) )) = ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2.data) + __pyx_t_17)) ))) - ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e1.data) + __pyx_t_18)) ))) * __pyx_v_cosbeta));
+086: e2mine1cosbeta[1] = e2[1] - e1[1]*cosbeta
__pyx_t_18 = 1;
__pyx_t_17 = 1;
__pyx_t_16 = 1;
*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2mine1cosbeta.data) + __pyx_t_16)) )) = ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2.data) + __pyx_t_18)) ))) - ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e1.data) + __pyx_t_17)) ))) * __pyx_v_cosbeta));
+087: e2mine1cosbeta[2] = e2[2] - e1[2]*cosbeta
__pyx_t_17 = 2;
__pyx_t_18 = 2;
__pyx_t_16 = 2;
*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2mine1cosbeta.data) + __pyx_t_16)) )) = ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2.data) + __pyx_t_17)) ))) - ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e1.data) + __pyx_t_18)) ))) * __pyx_v_cosbeta));
088:
+089: crosse1e2[0] = e1[1]*e2[2] - e1[2]*e2[1]
__pyx_t_18 = 1;
__pyx_t_17 = 2;
__pyx_t_16 = 2;
__pyx_t_15 = 1;
__pyx_t_13 = 0;
*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_crosse1e2.data) + __pyx_t_13)) )) = (((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e1.data) + __pyx_t_18)) ))) * (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2.data) + __pyx_t_17)) )))) - ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e1.data) + __pyx_t_16)) ))) * (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2.data) + __pyx_t_15)) )))));
+090: crosse1e2[1] = e1[2]*e2[0] - e1[0]*e2[2]
__pyx_t_15 = 2;
__pyx_t_16 = 0;
__pyx_t_17 = 0;
__pyx_t_18 = 2;
__pyx_t_13 = 1;
*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_crosse1e2.data) + __pyx_t_13)) )) = (((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e1.data) + __pyx_t_15)) ))) * (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2.data) + __pyx_t_16)) )))) - ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e1.data) + __pyx_t_17)) ))) * (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2.data) + __pyx_t_18)) )))));
+091: crosse1e2[2] = e1[0]*e2[1] - e1[1]*e2[0]
__pyx_t_18 = 0;
__pyx_t_17 = 1;
__pyx_t_16 = 1;
__pyx_t_15 = 0;
__pyx_t_13 = 2;
*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_crosse1e2.data) + __pyx_t_13)) )) = (((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e1.data) + __pyx_t_18)) ))) * (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2.data) + __pyx_t_17)) )))) - ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e1.data) + __pyx_t_16)) ))) * (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2.data) + __pyx_t_15)) )))));
092:
+093: a1 = ( r12[0]*e2cosbetamine1[0] + r12[1]*e2cosbetamine1[1] + r12[2]*e2cosbetamine1[2] ) / sinbetasq
__pyx_t_15 = 0;
__pyx_t_16 = 0;
__pyx_t_17 = 1;
__pyx_t_18 = 1;
__pyx_t_13 = 2;
__pyx_t_14 = 2;
__pyx_v_a1 = (((((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_r12.data) + __pyx_t_15)) ))) * (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2cosbetamine1.data) + __pyx_t_16)) )))) + ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_r12.data) + __pyx_t_17)) ))) * (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2cosbetamine1.data) + __pyx_t_18)) ))))) + ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_r12.data) + __pyx_t_13)) ))) * (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2cosbetamine1.data) + __pyx_t_14)) ))))) / __pyx_v_sinbetasq);
+094: a2 = ( r12[0]*e2mine1cosbeta[0] + r12[1]*e2mine1cosbeta[1] + r12[2]*e2mine1cosbeta[2] ) / sinbetasq
__pyx_t_14 = 0;
__pyx_t_13 = 0;
__pyx_t_18 = 1;
__pyx_t_17 = 1;
__pyx_t_16 = 2;
__pyx_t_15 = 2;
__pyx_v_a2 = (((((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_r12.data) + __pyx_t_14)) ))) * (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2mine1cosbeta.data) + __pyx_t_13)) )))) + ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_r12.data) + __pyx_t_18)) ))) * (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2mine1cosbeta.data) + __pyx_t_17)) ))))) + ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_r12.data) + __pyx_t_16)) ))) * (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_e2mine1cosbeta.data) + __pyx_t_15)) ))))) / __pyx_v_sinbetasq);
+095: a0 = ( r12[0]*crosse1e2[0] + r12[1]*crosse1e2[1] + r12[2]*crosse1e2[2] ) / sinbetasq
__pyx_t_15 = 0;
__pyx_t_16 = 0;
__pyx_t_17 = 1;
__pyx_t_18 = 1;
__pyx_t_13 = 2;
__pyx_t_14 = 2;
__pyx_v_a0 = (((((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_r12.data) + __pyx_t_15)) ))) * (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_crosse1e2.data) + __pyx_t_16)) )))) + ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_r12.data) + __pyx_t_17)) ))) * (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_crosse1e2.data) + __pyx_t_18)) ))))) + ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_r12.data) + __pyx_t_13)) ))) * (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_crosse1e2.data) + __pyx_t_14)) ))))) / __pyx_v_sinbetasq);
096:
097: # dlk = _F(a1+nv,a2+nw,a0,cosbeta,sinbetasq)
098: # dlk = dlk - _F(a1+nv,a2,a0,cosbeta,sinbetasq)
099: # dlk = dlk - _F(a1,a2+nw,a0,cosbeta,sinbetasq)
100: # dlk = dlk + _F(a1,a2,a0,cosbeta,sinbetasq)
101:
+102: dlk = 0
__pyx_v_dlk = 0.0;
+103: a0sq = a0**2
__pyx_v_a0sq = pow(__pyx_v_a0, 2.0);
+104: a0sqcosbeta = a0sq*cosbeta
__pyx_v_a0sqcosbeta = (__pyx_v_a0sq * __pyx_v_cosbeta);
+105: a0sqsinbetasq = a0sq*sinbetasq
__pyx_v_a0sqsinbetasq = (__pyx_v_a0sq * __pyx_v_sinbetasq);
106:
+107: t1=a1+nv
__pyx_v_t1 = (__pyx_v_a1 + __pyx_v_nv);
+108: t2=a2+nw
__pyx_v_t2 = (__pyx_v_a2 + __pyx_v_nw);
+109: num = t1*t2+a0sq*cosbeta
__pyx_v_num = ((__pyx_v_t1 * __pyx_v_t2) + (__pyx_v_a0sq * __pyx_v_cosbeta));
+110: denom = a0 * sqrt( t1**2 + t2**2 - 2*t1*t2*cosbeta + a0sqsinbetasq )
__pyx_v_denom = (__pyx_v_a0 * sqrt((((pow(__pyx_v_t1, 2.0) + pow(__pyx_v_t2, 2.0)) - (((2.0 * __pyx_v_t1) * __pyx_v_t2) * __pyx_v_cosbeta)) + __pyx_v_a0sqsinbetasq)));
+111: val = minoneover4pi * atan(num/denom)
__pyx_v_val = (__pyx_v_minoneover4pi * atan((__pyx_v_num / __pyx_v_denom)));
+112: if (not isnan(val)):
__pyx_t_19 = (!isnan(__pyx_v_val));
if (__pyx_t_19) {
/* … */
}
+113: dlk = dlk + val
__pyx_v_dlk = (__pyx_v_dlk + __pyx_v_val);
114:
+115: t1=a1+nv
__pyx_v_t1 = (__pyx_v_a1 + __pyx_v_nv);
+116: t2=a2
__pyx_v_t2 = __pyx_v_a2;
+117: num = t1*t2+a0sq*cosbeta
__pyx_v_num = ((__pyx_v_t1 * __pyx_v_t2) + (__pyx_v_a0sq * __pyx_v_cosbeta));
+118: denom = a0 * sqrt( t1**2 + t2**2 - 2*t1*t2*cosbeta + a0sqsinbetasq )
__pyx_v_denom = (__pyx_v_a0 * sqrt((((pow(__pyx_v_t1, 2.0) + pow(__pyx_v_t2, 2.0)) - (((2.0 * __pyx_v_t1) * __pyx_v_t2) * __pyx_v_cosbeta)) + __pyx_v_a0sqsinbetasq)));
+119: val = minoneover4pi * atan(num/denom)
__pyx_v_val = (__pyx_v_minoneover4pi * atan((__pyx_v_num / __pyx_v_denom)));
+120: if (not isnan(val)):
__pyx_t_19 = (!isnan(__pyx_v_val));
if (__pyx_t_19) {
/* … */
}
+121: dlk = dlk - val
__pyx_v_dlk = (__pyx_v_dlk - __pyx_v_val);
122:
+123: t1=a1
__pyx_v_t1 = __pyx_v_a1;
+124: t2=a2+nw
__pyx_v_t2 = (__pyx_v_a2 + __pyx_v_nw);
+125: num = t1*t2+a0sq*cosbeta
__pyx_v_num = ((__pyx_v_t1 * __pyx_v_t2) + (__pyx_v_a0sq * __pyx_v_cosbeta));
+126: denom = a0 * sqrt( t1**2 + t2**2 - 2*t1*t2*cosbeta + a0sqsinbetasq )
__pyx_v_denom = (__pyx_v_a0 * sqrt((((pow(__pyx_v_t1, 2.0) + pow(__pyx_v_t2, 2.0)) - (((2.0 * __pyx_v_t1) * __pyx_v_t2) * __pyx_v_cosbeta)) + __pyx_v_a0sqsinbetasq)));
+127: val = minoneover4pi * atan(num/denom)
__pyx_v_val = (__pyx_v_minoneover4pi * atan((__pyx_v_num / __pyx_v_denom)));
+128: if (not isnan(val)):
__pyx_t_19 = (!isnan(__pyx_v_val));
if (__pyx_t_19) {
/* … */
}
+129: dlk = dlk - val
__pyx_v_dlk = (__pyx_v_dlk - __pyx_v_val);
130:
+131: t1=a1
__pyx_v_t1 = __pyx_v_a1;
+132: t2=a2
__pyx_v_t2 = __pyx_v_a2;
+133: num = t1*t2+a0sq*cosbeta
__pyx_v_num = ((__pyx_v_t1 * __pyx_v_t2) + (__pyx_v_a0sq * __pyx_v_cosbeta));
+134: denom = a0 * sqrt( t1**2 + t2**2 - 2*t1*t2*cosbeta + a0sqsinbetasq )
__pyx_v_denom = (__pyx_v_a0 * sqrt((((pow(__pyx_v_t1, 2.0) + pow(__pyx_v_t2, 2.0)) - (((2.0 * __pyx_v_t1) * __pyx_v_t2) * __pyx_v_cosbeta)) + __pyx_v_a0sqsinbetasq)));
+135: val = minoneover4pi * atan(num/denom)
__pyx_v_val = (__pyx_v_minoneover4pi * atan((__pyx_v_num / __pyx_v_denom)));
+136: if (not isnan(val)):
__pyx_t_19 = (!isnan(__pyx_v_val));
if (__pyx_t_19) {
/* … */
}
+137: dlk = dlk + val
__pyx_v_dlk = (__pyx_v_dlk + __pyx_v_val);
138:
139: # dlk = 0
140: # val = _F(a1+nv,a2+nw,a0,cosbeta,sinbetasq)
141: # if val == val:
142: # dlk = dlk + val
143: # val = _F(a1+nv,a2,a0,cosbeta,sinbetasq)
144: # if val == val:
145: # dlk = dlk - val
146: # val = _F(a1,a2+nw,a0,cosbeta,sinbetasq)
147: # if val == val:
148: # dlk = dlk - val
149: # val = _F(a1,a2,a0,cosbeta,sinbetasq)
150: # if val == val:
151: # dlk = dlk + val
152:
+153: if (not isnan(dlk)):
__pyx_t_19 = (!isnan(__pyx_v_dlk));
if (__pyx_t_19) {
/* … */
}
}
}
+154: lk = lk + dlk
__pyx_v_lk = (__pyx_v_lk + __pyx_v_dlk);
+155: return lk
__pyx_r = __pyx_v_lk; goto __pyx_L0;
156:
157:
158: # @cython.boundscheck(False)
159: # @cython.wraparound(False)
160: # @cython.nonecheck(False)
161: # @cython.profile(True)
162: # @cython.cdivision(True)
163: # cpdef double _F(double t1, double t2, double a0, double cosbeta, double sinbetasq):
164: # cdef double minoneover4pi, a0sq, num, denom
165: # minoneover4pi = -0.07957747154594767
166: # a0sq = a0**2
167: # num = t1*t2+a0sq*cosbeta
168: # denom = a0 * sqrt( t1**2 + t2**2 - 2*t1*t2*cosbeta + a0sq*sinbetasq )
169: # return minoneover4pi * atan(num/denom)