+001: cimport cython
__pyx_t_1 = PyDict_New(); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 1, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
if (PyDict_SetItem(__pyx_d, __pyx_n_s_test, __pyx_t_1) < 0) __PYX_ERR(0, 1, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+002: import numpy as np
__pyx_t_1 = __Pyx_Import(__pyx_n_s_numpy, 0, -1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 2, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
if (PyDict_SetItem(__pyx_d, __pyx_n_s_np, __pyx_t_1) < 0) __PYX_ERR(0, 2, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
003: cimport numpy as cnp
+004: import math
__pyx_t_1 = __Pyx_Import(__pyx_n_s_math, 0, -1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 4, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
if (PyDict_SetItem(__pyx_d, __pyx_n_s_math, __pyx_t_1) < 0) __PYX_ERR(0, 4, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
005: from cpython cimport bool
006: from libc.math cimport sqrt as Csqrt, cos as Ccos, sin as Csin, atan2 as Catan2
007:
008: @cython.cdivision(True)
009: @cython.wraparound(False)
010: @cython.boundscheck(False)
+011: def Calc_LOS_PInOut_cython(double [:,::1] Ds, double [:,::1] us, double [:,::1] VPoly, double [:,::1] vIn,
/* Python wrapper */
static PyObject *__pyx_pw_46_cython_magic_163cfc3d2ce7aad39b88b29178aa2d8c_1Calc_LOS_PInOut_cython(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
static PyMethodDef __pyx_mdef_46_cython_magic_163cfc3d2ce7aad39b88b29178aa2d8c_1Calc_LOS_PInOut_cython = {"Calc_LOS_PInOut_cython", (PyCFunction)__pyx_pw_46_cython_magic_163cfc3d2ce7aad39b88b29178aa2d8c_1Calc_LOS_PInOut_cython, METH_VARARGS|METH_KEYWORDS, 0};
static PyObject *__pyx_pw_46_cython_magic_163cfc3d2ce7aad39b88b29178aa2d8c_1Calc_LOS_PInOut_cython(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
__Pyx_memviewslice __pyx_v_Ds = { 0, 0, { 0 }, { 0 }, { 0 } };
__Pyx_memviewslice __pyx_v_us = { 0, 0, { 0 }, { 0 }, { 0 } };
__Pyx_memviewslice __pyx_v_VPoly = { 0, 0, { 0 }, { 0 }, { 0 } };
__Pyx_memviewslice __pyx_v_vIn = { 0, 0, { 0 }, { 0 }, { 0 } };
PyBoolObject *__pyx_v_Forbid = 0;
PyObject *__pyx_v_RMin = 0;
double __pyx_v_EpsUz;
double __pyx_v_EpsVz;
double __pyx_v_EpsA;
double __pyx_v_EpsB;
PyObject *__pyx_r = 0;
__Pyx_RefNannyDeclarations
__Pyx_RefNannySetupContext("Calc_LOS_PInOut_cython (wrapper)", 0);
{
static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_Ds,&__pyx_n_s_us,&__pyx_n_s_VPoly,&__pyx_n_s_vIn,&__pyx_n_s_Forbid,&__pyx_n_s_RMin,&__pyx_n_s_EpsUz,&__pyx_n_s_EpsVz,&__pyx_n_s_EpsA,&__pyx_n_s_EpsB,0};
PyObject* values[10] = {0,0,0,0,0,0,0,0,0,0};
/* … */
/* function exit code */
goto __pyx_L0;
__pyx_L1_error:;
__pyx_r = NULL;
__pyx_L0:;
__Pyx_RefNannyFinishContext();
return __pyx_r;
}
static PyObject *__pyx_pf_46_cython_magic_163cfc3d2ce7aad39b88b29178aa2d8c_Calc_LOS_PInOut_cython(CYTHON_UNUSED PyObject *__pyx_self, __Pyx_memviewslice __pyx_v_Ds, __Pyx_memviewslice __pyx_v_us, __Pyx_memviewslice __pyx_v_VPoly, __Pyx_memviewslice __pyx_v_vIn, PyBoolObject *__pyx_v_Forbid, PyObject *__pyx_v_RMin, double __pyx_v_EpsUz, double __pyx_v_EpsVz, double __pyx_v_EpsA, double __pyx_v_EpsB) {
int __pyx_v_ii;
int __pyx_v_jj;
int __pyx_v_Nl;
int __pyx_v_Ns;
double __pyx_v_Rmin;
double __pyx_v_upscaDp;
double __pyx_v_upar2;
double __pyx_v_Dpar2;
double __pyx_v_Crit2;
double __pyx_v_kout;
double __pyx_v_kin;
int __pyx_v_indout;
int __pyx_v_Done;
double __pyx_v_L;
double __pyx_v_S1X;
double __pyx_v_S1Y;
double __pyx_v_S2X;
double __pyx_v_S2Y;
double __pyx_v_sca;
double __pyx_v_sca0;
double __pyx_v_sca1;
double __pyx_v_sca2;
double __pyx_v_q;
double __pyx_v_C;
double __pyx_v_delta;
double __pyx_v_sqd;
double __pyx_v_k;
double __pyx_v_sol0;
double __pyx_v_sol1;
double __pyx_v_v0;
double __pyx_v_v1;
double __pyx_v_A;
double __pyx_v_B;
int __pyx_v_Forbidbis;
int __pyx_v_Forbid0;
PyArrayObject *__pyx_v_SIn_ = 0;
PyArrayObject *__pyx_v_SOut_ = 0;
PyArrayObject *__pyx_v_VPerp_ = 0;
PyArrayObject *__pyx_v_indOut_ = 0;
__Pyx_memviewslice __pyx_v_SIn = { 0, 0, { 0 }, { 0 }, { 0 } };
__Pyx_memviewslice __pyx_v_SOut = { 0, 0, { 0 }, { 0 }, { 0 } };
__Pyx_memviewslice __pyx_v_VPerp = { 0, 0, { 0 }, { 0 }, { 0 } };
__Pyx_memviewslice __pyx_v_indOut = { 0, 0, { 0 }, { 0 }, { 0 } };
double __pyx_v_phi;
__Pyx_LocalBuf_ND __pyx_pybuffernd_SIn_;
__Pyx_Buffer __pyx_pybuffer_SIn_;
__Pyx_LocalBuf_ND __pyx_pybuffernd_SOut_;
__Pyx_Buffer __pyx_pybuffer_SOut_;
__Pyx_LocalBuf_ND __pyx_pybuffernd_VPerp_;
__Pyx_Buffer __pyx_pybuffer_VPerp_;
__Pyx_LocalBuf_ND __pyx_pybuffernd_indOut_;
__Pyx_Buffer __pyx_pybuffer_indOut_;
PyObject *__pyx_r = NULL;
__Pyx_RefNannyDeclarations
__Pyx_RefNannySetupContext("Calc_LOS_PInOut_cython", 0);
__pyx_pybuffer_SIn_.pybuffer.buf = NULL;
__pyx_pybuffer_SIn_.refcount = 0;
__pyx_pybuffernd_SIn_.data = NULL;
__pyx_pybuffernd_SIn_.rcbuffer = &__pyx_pybuffer_SIn_;
__pyx_pybuffer_SOut_.pybuffer.buf = NULL;
__pyx_pybuffer_SOut_.refcount = 0;
__pyx_pybuffernd_SOut_.data = NULL;
__pyx_pybuffernd_SOut_.rcbuffer = &__pyx_pybuffer_SOut_;
__pyx_pybuffer_VPerp_.pybuffer.buf = NULL;
__pyx_pybuffer_VPerp_.refcount = 0;
__pyx_pybuffernd_VPerp_.data = NULL;
__pyx_pybuffernd_VPerp_.rcbuffer = &__pyx_pybuffer_VPerp_;
__pyx_pybuffer_indOut_.pybuffer.buf = NULL;
__pyx_pybuffer_indOut_.refcount = 0;
__pyx_pybuffernd_indOut_.data = NULL;
__pyx_pybuffernd_indOut_.rcbuffer = &__pyx_pybuffer_indOut_;
/* … */
/* function exit code */
__pyx_L1_error:;
__Pyx_XDECREF(__pyx_t_1);
__Pyx_XDECREF(__pyx_t_2);
__Pyx_XDECREF(__pyx_t_3);
__Pyx_XDECREF(__pyx_t_4);
__Pyx_XDECREF(__pyx_t_5);
__Pyx_XDECREF(__pyx_t_6);
__PYX_XDEC_MEMVIEW(&__pyx_t_10, 1);
__PYX_XDEC_MEMVIEW(&__pyx_t_11, 1);
__Pyx_XDECREF(__pyx_t_14);
__Pyx_XDECREF(__pyx_t_16);
__Pyx_XDECREF(__pyx_t_254);
__Pyx_XDECREF(__pyx_t_255);
{ 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_SIn_.rcbuffer->pybuffer);
__Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_SOut_.rcbuffer->pybuffer);
__Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_VPerp_.rcbuffer->pybuffer);
__Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_indOut_.rcbuffer->pybuffer);
__Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
__Pyx_AddTraceback("_cython_magic_163cfc3d2ce7aad39b88b29178aa2d8c.Calc_LOS_PInOut_cython", __pyx_clineno, __pyx_lineno, __pyx_filename);
__pyx_r = NULL;
goto __pyx_L2;
__pyx_L0:;
__Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_SIn_.rcbuffer->pybuffer);
__Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_SOut_.rcbuffer->pybuffer);
__Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_VPerp_.rcbuffer->pybuffer);
__Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_indOut_.rcbuffer->pybuffer);
__pyx_L2:;
__Pyx_XDECREF((PyObject *)__pyx_v_SIn_);
__Pyx_XDECREF((PyObject *)__pyx_v_SOut_);
__Pyx_XDECREF((PyObject *)__pyx_v_VPerp_);
__Pyx_XDECREF((PyObject *)__pyx_v_indOut_);
__PYX_XDEC_MEMVIEW(&__pyx_v_SIn, 1);
__PYX_XDEC_MEMVIEW(&__pyx_v_SOut, 1);
__PYX_XDEC_MEMVIEW(&__pyx_v_VPerp, 1);
__PYX_XDEC_MEMVIEW(&__pyx_v_indOut, 1);
__PYX_XDEC_MEMVIEW(&__pyx_v_Ds, 1);
__PYX_XDEC_MEMVIEW(&__pyx_v_us, 1);
__PYX_XDEC_MEMVIEW(&__pyx_v_VPoly, 1);
__PYX_XDEC_MEMVIEW(&__pyx_v_vIn, 1);
__Pyx_XGIVEREF(__pyx_r);
__Pyx_RefNannyFinishContext();
return __pyx_r;
}
/* … */
__pyx_tuple__23 = PyTuple_Pack(54, __pyx_n_s_Ds, __pyx_n_s_us, __pyx_n_s_VPoly, __pyx_n_s_vIn, __pyx_n_s_Forbid, __pyx_n_s_RMin, __pyx_n_s_EpsUz, __pyx_n_s_EpsVz, __pyx_n_s_EpsA, __pyx_n_s_EpsB, __pyx_n_s_ii, __pyx_n_s_jj, __pyx_n_s_Nl, __pyx_n_s_Ns, __pyx_n_s_Rmin, __pyx_n_s_upscaDp, __pyx_n_s_upar2, __pyx_n_s_Dpar2, __pyx_n_s_Crit2, __pyx_n_s_kout, __pyx_n_s_kin, __pyx_n_s_indout, __pyx_n_s_Done, __pyx_n_s_L, __pyx_n_s_S1X, __pyx_n_s_S1Y, __pyx_n_s_S2X, __pyx_n_s_S2Y, __pyx_n_s_sca, __pyx_n_s_sca0, __pyx_n_s_sca1, __pyx_n_s_sca2, __pyx_n_s_q, __pyx_n_s_C, __pyx_n_s_delta, __pyx_n_s_sqd, __pyx_n_s_k, __pyx_n_s_sol0, __pyx_n_s_sol1, __pyx_n_s_v0, __pyx_n_s_v1, __pyx_n_s_A, __pyx_n_s_B, __pyx_n_s_Forbidbis, __pyx_n_s_Forbid0, __pyx_n_s_SIn, __pyx_n_s_SOut, __pyx_n_s_VPerp, __pyx_n_s_indOut, __pyx_n_s_SIn_2, __pyx_n_s_SOut_2, __pyx_n_s_VPerp_2, __pyx_n_s_indOut_2, __pyx_n_s_phi); if (unlikely(!__pyx_tuple__23)) __PYX_ERR(0, 11, __pyx_L1_error)
__Pyx_GOTREF(__pyx_tuple__23);
__Pyx_GIVEREF(__pyx_tuple__23);
/* … */
__pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_46_cython_magic_163cfc3d2ce7aad39b88b29178aa2d8c_1Calc_LOS_PInOut_cython, NULL, __pyx_n_s_cython_magic_163cfc3d2ce7aad39b); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 11, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
if (PyDict_SetItem(__pyx_d, __pyx_n_s_Calc_LOS_PInOut_cython, __pyx_t_1) < 0) __PYX_ERR(0, 11, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__pyx_codeobj__24 = (PyObject*)__Pyx_PyCode_New(10, 0, 54, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__23, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_Home_DV226270_cache_ipython_cyt, __pyx_n_s_Calc_LOS_PInOut_cython, 11, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__24)) __PYX_ERR(0, 11, __pyx_L1_error)
+012: bool Forbid=True, RMin=None, double EpsUz=1.e-6, double EpsVz=1.e-9, double EpsA=1.e-9, double EpsB=1.e-9):
values[4] = (PyObject *)((PyBoolObject *)Py_True);
values[5] = ((PyObject *)Py_None);
if (unlikely(__pyx_kwds)) {
Py_ssize_t kw_args;
const Py_ssize_t pos_args = PyTuple_GET_SIZE(__pyx_args);
switch (pos_args) {
case 10: values[9] = PyTuple_GET_ITEM(__pyx_args, 9);
case 9: values[8] = PyTuple_GET_ITEM(__pyx_args, 8);
case 8: values[7] = PyTuple_GET_ITEM(__pyx_args, 7);
case 7: values[6] = PyTuple_GET_ITEM(__pyx_args, 6);
case 6: values[5] = PyTuple_GET_ITEM(__pyx_args, 5);
case 5: values[4] = PyTuple_GET_ITEM(__pyx_args, 4);
case 4: values[3] = PyTuple_GET_ITEM(__pyx_args, 3);
case 3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
case 2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
case 1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
case 0: break;
default: goto __pyx_L5_argtuple_error;
}
kw_args = PyDict_Size(__pyx_kwds);
switch (pos_args) {
case 0:
if (likely((values[0] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_Ds)) != 0)) kw_args--;
else goto __pyx_L5_argtuple_error;
case 1:
if (likely((values[1] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_us)) != 0)) kw_args--;
else {
__Pyx_RaiseArgtupleInvalid("Calc_LOS_PInOut_cython", 0, 4, 10, 1); __PYX_ERR(0, 11, __pyx_L3_error)
}
case 2:
if (likely((values[2] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_VPoly)) != 0)) kw_args--;
else {
__Pyx_RaiseArgtupleInvalid("Calc_LOS_PInOut_cython", 0, 4, 10, 2); __PYX_ERR(0, 11, __pyx_L3_error)
}
case 3:
if (likely((values[3] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_vIn)) != 0)) kw_args--;
else {
__Pyx_RaiseArgtupleInvalid("Calc_LOS_PInOut_cython", 0, 4, 10, 3); __PYX_ERR(0, 11, __pyx_L3_error)
}
case 4:
if (kw_args > 0) {
PyObject* value = PyDict_GetItem(__pyx_kwds, __pyx_n_s_Forbid);
if (value) { values[4] = value; kw_args--; }
}
case 5:
if (kw_args > 0) {
PyObject* value = PyDict_GetItem(__pyx_kwds, __pyx_n_s_RMin);
if (value) { values[5] = value; kw_args--; }
}
case 6:
if (kw_args > 0) {
PyObject* value = PyDict_GetItem(__pyx_kwds, __pyx_n_s_EpsUz);
if (value) { values[6] = value; kw_args--; }
}
case 7:
if (kw_args > 0) {
PyObject* value = PyDict_GetItem(__pyx_kwds, __pyx_n_s_EpsVz);
if (value) { values[7] = value; kw_args--; }
}
case 8:
if (kw_args > 0) {
PyObject* value = PyDict_GetItem(__pyx_kwds, __pyx_n_s_EpsA);
if (value) { values[8] = value; kw_args--; }
}
case 9:
if (kw_args > 0) {
PyObject* value = PyDict_GetItem(__pyx_kwds, __pyx_n_s_EpsB);
if (value) { values[9] = value; kw_args--; }
}
}
if (unlikely(kw_args > 0)) {
if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "Calc_LOS_PInOut_cython") < 0)) __PYX_ERR(0, 11, __pyx_L3_error)
}
} else {
switch (PyTuple_GET_SIZE(__pyx_args)) {
case 10: values[9] = PyTuple_GET_ITEM(__pyx_args, 9);
case 9: values[8] = PyTuple_GET_ITEM(__pyx_args, 8);
case 8: values[7] = PyTuple_GET_ITEM(__pyx_args, 7);
case 7: values[6] = PyTuple_GET_ITEM(__pyx_args, 6);
case 6: values[5] = PyTuple_GET_ITEM(__pyx_args, 5);
case 5: values[4] = PyTuple_GET_ITEM(__pyx_args, 4);
case 4: values[3] = PyTuple_GET_ITEM(__pyx_args, 3);
values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
break;
default: goto __pyx_L5_argtuple_error;
}
}
__pyx_v_Ds = __Pyx_PyObject_to_MemoryviewSlice_d_dc_double(values[0]); if (unlikely(!__pyx_v_Ds.memview)) __PYX_ERR(0, 11, __pyx_L3_error)
__pyx_v_us = __Pyx_PyObject_to_MemoryviewSlice_d_dc_double(values[1]); if (unlikely(!__pyx_v_us.memview)) __PYX_ERR(0, 11, __pyx_L3_error)
__pyx_v_VPoly = __Pyx_PyObject_to_MemoryviewSlice_d_dc_double(values[2]); if (unlikely(!__pyx_v_VPoly.memview)) __PYX_ERR(0, 11, __pyx_L3_error)
__pyx_v_vIn = __Pyx_PyObject_to_MemoryviewSlice_d_dc_double(values[3]); if (unlikely(!__pyx_v_vIn.memview)) __PYX_ERR(0, 11, __pyx_L3_error)
__pyx_v_Forbid = ((PyBoolObject *)values[4]);
__pyx_v_RMin = values[5];
if (values[6]) {
__pyx_v_EpsUz = __pyx_PyFloat_AsDouble(values[6]); if (unlikely((__pyx_v_EpsUz == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 12, __pyx_L3_error)
} else {
__pyx_v_EpsUz = ((double)1.e-6);
}
if (values[7]) {
__pyx_v_EpsVz = __pyx_PyFloat_AsDouble(values[7]); if (unlikely((__pyx_v_EpsVz == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 12, __pyx_L3_error)
} else {
__pyx_v_EpsVz = ((double)1.e-9);
}
if (values[8]) {
__pyx_v_EpsA = __pyx_PyFloat_AsDouble(values[8]); if (unlikely((__pyx_v_EpsA == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 12, __pyx_L3_error)
} else {
__pyx_v_EpsA = ((double)1.e-9);
}
if (values[9]) {
__pyx_v_EpsB = __pyx_PyFloat_AsDouble(values[9]); if (unlikely((__pyx_v_EpsB == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 12, __pyx_L3_error)
} else {
__pyx_v_EpsB = ((double)1.e-9);
}
}
goto __pyx_L4_argument_unpacking_done;
__pyx_L5_argtuple_error:;
__Pyx_RaiseArgtupleInvalid("Calc_LOS_PInOut_cython", 0, 4, 10, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 11, __pyx_L3_error)
__pyx_L3_error:;
__Pyx_AddTraceback("_cython_magic_163cfc3d2ce7aad39b88b29178aa2d8c.Calc_LOS_PInOut_cython", __pyx_clineno, __pyx_lineno, __pyx_filename);
__Pyx_RefNannyFinishContext();
return NULL;
__pyx_L4_argument_unpacking_done:;
if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_Forbid), __pyx_ptype_7cpython_4bool_bool, 1, "Forbid", 0))) __PYX_ERR(0, 12, __pyx_L1_error)
__pyx_r = __pyx_pf_46_cython_magic_163cfc3d2ce7aad39b88b29178aa2d8c_Calc_LOS_PInOut_cython(__pyx_self, __pyx_v_Ds, __pyx_v_us, __pyx_v_VPoly, __pyx_v_vIn, __pyx_v_Forbid, __pyx_v_RMin, __pyx_v_EpsUz, __pyx_v_EpsVz, __pyx_v_EpsA, __pyx_v_EpsB);
013:
+014: cdef int ii, jj, Nl=Ds.shape[1], Ns=vIn.shape[1]
__pyx_v_Nl = (__pyx_v_Ds.shape[1]);
__pyx_v_Ns = (__pyx_v_vIn.shape[1]);
015: cdef double Rmin, upscaDp, upar2, Dpar2, Crit2, kout, kin
016: cdef int indout, Done
017: cdef double L, S1X, S1Y, S2X, S2Y, sca, sca0, sca1, sca2
018: cdef double q, C, delta, sqd, k, sol0, sol1
019: cdef double v0, v1, A, B
020: cdef int Forbidbis, Forbid0
+021: cdef cnp.ndarray[double,ndim=2] SIn_=np.nan*np.ones((3,Nl)), SOut_=np.nan*np.ones((3,Nl))
__pyx_t_1 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 21, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_nan); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 21, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__pyx_t_3 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 21, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_ones); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 21, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__pyx_t_3 = __Pyx_PyInt_From_int(__pyx_v_Nl); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 21, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__pyx_t_5 = PyTuple_New(2); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 21, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__Pyx_INCREF(__pyx_int_3);
__Pyx_GIVEREF(__pyx_int_3);
PyTuple_SET_ITEM(__pyx_t_5, 0, __pyx_int_3);
__Pyx_GIVEREF(__pyx_t_3);
PyTuple_SET_ITEM(__pyx_t_5, 1, __pyx_t_3);
__pyx_t_3 = 0;
__pyx_t_3 = NULL;
if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_4))) {
__pyx_t_3 = PyMethod_GET_SELF(__pyx_t_4);
if (likely(__pyx_t_3)) {
PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_4);
__Pyx_INCREF(__pyx_t_3);
__Pyx_INCREF(function);
__Pyx_DECREF_SET(__pyx_t_4, function);
}
}
if (!__pyx_t_3) {
__pyx_t_1 = __Pyx_PyObject_CallOneArg(__pyx_t_4, __pyx_t_5); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 21, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
__Pyx_GOTREF(__pyx_t_1);
} else {
#if CYTHON_FAST_PYCALL
if (PyFunction_Check(__pyx_t_4)) {
PyObject *__pyx_temp[2] = {__pyx_t_3, __pyx_t_5};
__pyx_t_1 = __Pyx_PyFunction_FastCall(__pyx_t_4, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 21, __pyx_L1_error)
__Pyx_XDECREF(__pyx_t_3); __pyx_t_3 = 0;
__Pyx_GOTREF(__pyx_t_1);
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
} else
#endif
#if CYTHON_FAST_PYCCALL
if (__Pyx_PyFastCFunction_Check(__pyx_t_4)) {
PyObject *__pyx_temp[2] = {__pyx_t_3, __pyx_t_5};
__pyx_t_1 = __Pyx_PyCFunction_FastCall(__pyx_t_4, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 21, __pyx_L1_error)
__Pyx_XDECREF(__pyx_t_3); __pyx_t_3 = 0;
__Pyx_GOTREF(__pyx_t_1);
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
} else
#endif
{
__pyx_t_6 = PyTuple_New(1+1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 21, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_6);
__Pyx_GIVEREF(__pyx_t_3); PyTuple_SET_ITEM(__pyx_t_6, 0, __pyx_t_3); __pyx_t_3 = NULL;
__Pyx_GIVEREF(__pyx_t_5);
PyTuple_SET_ITEM(__pyx_t_6, 0+1, __pyx_t_5);
__pyx_t_5 = 0;
__pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_4, __pyx_t_6, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 21, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
}
}
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
__pyx_t_4 = PyNumber_Multiply(__pyx_t_2, __pyx_t_1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 21, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
if (!(likely(((__pyx_t_4) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_4, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 21, __pyx_L1_error)
__pyx_t_7 = ((PyArrayObject *)__pyx_t_4);
{
__Pyx_BufFmt_StackElem __pyx_stack[1];
if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_SIn_.rcbuffer->pybuffer, (PyObject*)__pyx_t_7, &__Pyx_TypeInfo_double, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
__pyx_v_SIn_ = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_SIn_.rcbuffer->pybuffer.buf = NULL;
__PYX_ERR(0, 21, __pyx_L1_error)
} else {__pyx_pybuffernd_SIn_.diminfo[0].strides = __pyx_pybuffernd_SIn_.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_SIn_.diminfo[0].shape = __pyx_pybuffernd_SIn_.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_SIn_.diminfo[1].strides = __pyx_pybuffernd_SIn_.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_SIn_.diminfo[1].shape = __pyx_pybuffernd_SIn_.rcbuffer->pybuffer.shape[1];
}
}
__pyx_t_7 = 0;
__pyx_v_SIn_ = ((PyArrayObject *)__pyx_t_4);
__pyx_t_4 = 0;
__pyx_t_4 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 21, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_4, __pyx_n_s_nan); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 21, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
__pyx_t_2 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 21, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__pyx_t_6 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_ones); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 21, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_6);
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
__pyx_t_2 = __Pyx_PyInt_From_int(__pyx_v_Nl); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 21, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__pyx_t_5 = PyTuple_New(2); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 21, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__Pyx_INCREF(__pyx_int_3);
__Pyx_GIVEREF(__pyx_int_3);
PyTuple_SET_ITEM(__pyx_t_5, 0, __pyx_int_3);
__Pyx_GIVEREF(__pyx_t_2);
PyTuple_SET_ITEM(__pyx_t_5, 1, __pyx_t_2);
__pyx_t_2 = 0;
__pyx_t_2 = NULL;
if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_6))) {
__pyx_t_2 = PyMethod_GET_SELF(__pyx_t_6);
if (likely(__pyx_t_2)) {
PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_6);
__Pyx_INCREF(__pyx_t_2);
__Pyx_INCREF(function);
__Pyx_DECREF_SET(__pyx_t_6, function);
}
}
if (!__pyx_t_2) {
__pyx_t_4 = __Pyx_PyObject_CallOneArg(__pyx_t_6, __pyx_t_5); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 21, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
__Pyx_GOTREF(__pyx_t_4);
} else {
#if CYTHON_FAST_PYCALL
if (PyFunction_Check(__pyx_t_6)) {
PyObject *__pyx_temp[2] = {__pyx_t_2, __pyx_t_5};
__pyx_t_4 = __Pyx_PyFunction_FastCall(__pyx_t_6, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 21, __pyx_L1_error)
__Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
__Pyx_GOTREF(__pyx_t_4);
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
} else
#endif
#if CYTHON_FAST_PYCCALL
if (__Pyx_PyFastCFunction_Check(__pyx_t_6)) {
PyObject *__pyx_temp[2] = {__pyx_t_2, __pyx_t_5};
__pyx_t_4 = __Pyx_PyCFunction_FastCall(__pyx_t_6, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 21, __pyx_L1_error)
__Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
__Pyx_GOTREF(__pyx_t_4);
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
} else
#endif
{
__pyx_t_3 = PyTuple_New(1+1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 21, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__Pyx_GIVEREF(__pyx_t_2); PyTuple_SET_ITEM(__pyx_t_3, 0, __pyx_t_2); __pyx_t_2 = NULL;
__Pyx_GIVEREF(__pyx_t_5);
PyTuple_SET_ITEM(__pyx_t_3, 0+1, __pyx_t_5);
__pyx_t_5 = 0;
__pyx_t_4 = __Pyx_PyObject_Call(__pyx_t_6, __pyx_t_3, NULL); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 21, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
}
}
__Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
__pyx_t_6 = PyNumber_Multiply(__pyx_t_1, __pyx_t_4); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 21, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_6);
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
if (!(likely(((__pyx_t_6) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_6, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 21, __pyx_L1_error)
__pyx_t_7 = ((PyArrayObject *)__pyx_t_6);
{
__Pyx_BufFmt_StackElem __pyx_stack[1];
if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_SOut_.rcbuffer->pybuffer, (PyObject*)__pyx_t_7, &__Pyx_TypeInfo_double, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
__pyx_v_SOut_ = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_SOut_.rcbuffer->pybuffer.buf = NULL;
__PYX_ERR(0, 21, __pyx_L1_error)
} else {__pyx_pybuffernd_SOut_.diminfo[0].strides = __pyx_pybuffernd_SOut_.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_SOut_.diminfo[0].shape = __pyx_pybuffernd_SOut_.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_SOut_.diminfo[1].strides = __pyx_pybuffernd_SOut_.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_SOut_.diminfo[1].shape = __pyx_pybuffernd_SOut_.rcbuffer->pybuffer.shape[1];
}
}
__pyx_t_7 = 0;
__pyx_v_SOut_ = ((PyArrayObject *)__pyx_t_6);
__pyx_t_6 = 0;
+022: cdef cnp.ndarray[double,ndim=2] VPerp_=np.nan*np.ones((3,Nl))
__pyx_t_6 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 22, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_6);
__pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_6, __pyx_n_s_nan); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 22, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
__pyx_t_1 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 22, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_ones); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 22, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__pyx_t_1 = __Pyx_PyInt_From_int(__pyx_v_Nl); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 22, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__pyx_t_5 = PyTuple_New(2); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 22, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__Pyx_INCREF(__pyx_int_3);
__Pyx_GIVEREF(__pyx_int_3);
PyTuple_SET_ITEM(__pyx_t_5, 0, __pyx_int_3);
__Pyx_GIVEREF(__pyx_t_1);
PyTuple_SET_ITEM(__pyx_t_5, 1, __pyx_t_1);
__pyx_t_1 = 0;
__pyx_t_1 = NULL;
if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_3))) {
__pyx_t_1 = PyMethod_GET_SELF(__pyx_t_3);
if (likely(__pyx_t_1)) {
PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_3);
__Pyx_INCREF(__pyx_t_1);
__Pyx_INCREF(function);
__Pyx_DECREF_SET(__pyx_t_3, function);
}
}
if (!__pyx_t_1) {
__pyx_t_6 = __Pyx_PyObject_CallOneArg(__pyx_t_3, __pyx_t_5); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 22, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
__Pyx_GOTREF(__pyx_t_6);
} else {
#if CYTHON_FAST_PYCALL
if (PyFunction_Check(__pyx_t_3)) {
PyObject *__pyx_temp[2] = {__pyx_t_1, __pyx_t_5};
__pyx_t_6 = __Pyx_PyFunction_FastCall(__pyx_t_3, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 22, __pyx_L1_error)
__Pyx_XDECREF(__pyx_t_1); __pyx_t_1 = 0;
__Pyx_GOTREF(__pyx_t_6);
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
} else
#endif
#if CYTHON_FAST_PYCCALL
if (__Pyx_PyFastCFunction_Check(__pyx_t_3)) {
PyObject *__pyx_temp[2] = {__pyx_t_1, __pyx_t_5};
__pyx_t_6 = __Pyx_PyCFunction_FastCall(__pyx_t_3, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 22, __pyx_L1_error)
__Pyx_XDECREF(__pyx_t_1); __pyx_t_1 = 0;
__Pyx_GOTREF(__pyx_t_6);
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
} else
#endif
{
__pyx_t_2 = PyTuple_New(1+1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 22, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__Pyx_GIVEREF(__pyx_t_1); PyTuple_SET_ITEM(__pyx_t_2, 0, __pyx_t_1); __pyx_t_1 = NULL;
__Pyx_GIVEREF(__pyx_t_5);
PyTuple_SET_ITEM(__pyx_t_2, 0+1, __pyx_t_5);
__pyx_t_5 = 0;
__pyx_t_6 = __Pyx_PyObject_Call(__pyx_t_3, __pyx_t_2, NULL); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 22, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_6);
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
}
}
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__pyx_t_3 = PyNumber_Multiply(__pyx_t_4, __pyx_t_6); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 22, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
__Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
if (!(likely(((__pyx_t_3) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_3, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 22, __pyx_L1_error)
__pyx_t_8 = ((PyArrayObject *)__pyx_t_3);
{
__Pyx_BufFmt_StackElem __pyx_stack[1];
if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_VPerp_.rcbuffer->pybuffer, (PyObject*)__pyx_t_8, &__Pyx_TypeInfo_double, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
__pyx_v_VPerp_ = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_VPerp_.rcbuffer->pybuffer.buf = NULL;
__PYX_ERR(0, 22, __pyx_L1_error)
} else {__pyx_pybuffernd_VPerp_.diminfo[0].strides = __pyx_pybuffernd_VPerp_.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_VPerp_.diminfo[0].shape = __pyx_pybuffernd_VPerp_.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_VPerp_.diminfo[1].strides = __pyx_pybuffernd_VPerp_.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_VPerp_.diminfo[1].shape = __pyx_pybuffernd_VPerp_.rcbuffer->pybuffer.shape[1];
}
}
__pyx_t_8 = 0;
__pyx_v_VPerp_ = ((PyArrayObject *)__pyx_t_3);
__pyx_t_3 = 0;
+023: cdef cnp.ndarray[double,ndim=1] indOut_=np.nan*np.ones((Nl,))
__pyx_t_3 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 23, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__pyx_t_6 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_nan); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 23, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_6);
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__pyx_t_4 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 23, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_4, __pyx_n_s_ones); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 23, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
__pyx_t_4 = __Pyx_PyInt_From_int(__pyx_v_Nl); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 23, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__pyx_t_5 = PyTuple_New(1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 23, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__Pyx_GIVEREF(__pyx_t_4);
PyTuple_SET_ITEM(__pyx_t_5, 0, __pyx_t_4);
__pyx_t_4 = 0;
__pyx_t_4 = NULL;
if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_2))) {
__pyx_t_4 = PyMethod_GET_SELF(__pyx_t_2);
if (likely(__pyx_t_4)) {
PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_2);
__Pyx_INCREF(__pyx_t_4);
__Pyx_INCREF(function);
__Pyx_DECREF_SET(__pyx_t_2, function);
}
}
if (!__pyx_t_4) {
__pyx_t_3 = __Pyx_PyObject_CallOneArg(__pyx_t_2, __pyx_t_5); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 23, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
__Pyx_GOTREF(__pyx_t_3);
} else {
#if CYTHON_FAST_PYCALL
if (PyFunction_Check(__pyx_t_2)) {
PyObject *__pyx_temp[2] = {__pyx_t_4, __pyx_t_5};
__pyx_t_3 = __Pyx_PyFunction_FastCall(__pyx_t_2, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 23, __pyx_L1_error)
__Pyx_XDECREF(__pyx_t_4); __pyx_t_4 = 0;
__Pyx_GOTREF(__pyx_t_3);
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
} else
#endif
#if CYTHON_FAST_PYCCALL
if (__Pyx_PyFastCFunction_Check(__pyx_t_2)) {
PyObject *__pyx_temp[2] = {__pyx_t_4, __pyx_t_5};
__pyx_t_3 = __Pyx_PyCFunction_FastCall(__pyx_t_2, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 23, __pyx_L1_error)
__Pyx_XDECREF(__pyx_t_4); __pyx_t_4 = 0;
__Pyx_GOTREF(__pyx_t_3);
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
} else
#endif
{
__pyx_t_1 = PyTuple_New(1+1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 23, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__Pyx_GIVEREF(__pyx_t_4); PyTuple_SET_ITEM(__pyx_t_1, 0, __pyx_t_4); __pyx_t_4 = NULL;
__Pyx_GIVEREF(__pyx_t_5);
PyTuple_SET_ITEM(__pyx_t_1, 0+1, __pyx_t_5);
__pyx_t_5 = 0;
__pyx_t_3 = __Pyx_PyObject_Call(__pyx_t_2, __pyx_t_1, NULL); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 23, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
}
}
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
__pyx_t_2 = PyNumber_Multiply(__pyx_t_6, __pyx_t_3); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 23, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
if (!(likely(((__pyx_t_2) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_2, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 23, __pyx_L1_error)
__pyx_t_9 = ((PyArrayObject *)__pyx_t_2);
{
__Pyx_BufFmt_StackElem __pyx_stack[1];
if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_indOut_.rcbuffer->pybuffer, (PyObject*)__pyx_t_9, &__Pyx_TypeInfo_double, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {
__pyx_v_indOut_ = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_indOut_.rcbuffer->pybuffer.buf = NULL;
__PYX_ERR(0, 23, __pyx_L1_error)
} else {__pyx_pybuffernd_indOut_.diminfo[0].strides = __pyx_pybuffernd_indOut_.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_indOut_.diminfo[0].shape = __pyx_pybuffernd_indOut_.rcbuffer->pybuffer.shape[0];
}
}
__pyx_t_9 = 0;
__pyx_v_indOut_ = ((PyArrayObject *)__pyx_t_2);
__pyx_t_2 = 0;
024:
+025: cdef double[:,::1] SIn=SIn_, SOut=SOut_, VPerp=VPerp_
__pyx_t_10 = __Pyx_PyObject_to_MemoryviewSlice_d_dc_double(((PyObject *)__pyx_v_SIn_));
if (unlikely(!__pyx_t_10.memview)) __PYX_ERR(0, 25, __pyx_L1_error)
__pyx_v_SIn = __pyx_t_10;
__pyx_t_10.memview = NULL;
__pyx_t_10.data = NULL;
__pyx_t_10 = __Pyx_PyObject_to_MemoryviewSlice_d_dc_double(((PyObject *)__pyx_v_SOut_));
if (unlikely(!__pyx_t_10.memview)) __PYX_ERR(0, 25, __pyx_L1_error)
__pyx_v_SOut = __pyx_t_10;
__pyx_t_10.memview = NULL;
__pyx_t_10.data = NULL;
__pyx_t_10 = __Pyx_PyObject_to_MemoryviewSlice_d_dc_double(((PyObject *)__pyx_v_VPerp_));
if (unlikely(!__pyx_t_10.memview)) __PYX_ERR(0, 25, __pyx_L1_error)
__pyx_v_VPerp = __pyx_t_10;
__pyx_t_10.memview = NULL;
__pyx_t_10.data = NULL;
+026: cdef double[::1] indOut=indOut_
__pyx_t_11 = __Pyx_PyObject_to_MemoryviewSlice_dc_double(((PyObject *)__pyx_v_indOut_));
if (unlikely(!__pyx_t_11.memview)) __PYX_ERR(0, 26, __pyx_L1_error)
__pyx_v_indOut = __pyx_t_11;
__pyx_t_11.memview = NULL;
__pyx_t_11.data = NULL;
027:
028: ################
029: # Prepare input
+030: if RMin is None:
__pyx_t_12 = (__pyx_v_RMin == Py_None);
__pyx_t_13 = (__pyx_t_12 != 0);
if (__pyx_t_13) {
/* … */
goto __pyx_L3;
}
+031: Rmin = 0.95*min(np.min(VPoly[0,:]),np.min(np.hypot(Ds[0,:],Ds[1,:])))
__pyx_t_3 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 31, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__pyx_t_6 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_min); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 31, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_6);
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__pyx_t_1 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 31, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_hypot); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 31, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__pyx_t_11.data = __pyx_v_Ds.data;
__pyx_t_11.memview = __pyx_v_Ds.memview;
__PYX_INC_MEMVIEW(&__pyx_t_11, 0);
{
Py_ssize_t __pyx_tmp_idx = 0;
Py_ssize_t __pyx_tmp_shape = __pyx_v_Ds.shape[0];
Py_ssize_t __pyx_tmp_stride = __pyx_v_Ds.strides[0];
if (0 && (__pyx_tmp_idx < 0))
__pyx_tmp_idx += __pyx_tmp_shape;
if (0 && (__pyx_tmp_idx < 0 || __pyx_tmp_idx >= __pyx_tmp_shape)) {
PyErr_SetString(PyExc_IndexError, "Index out of bounds (axis 0)");
__PYX_ERR(0, 31, __pyx_L1_error)
}
__pyx_t_11.data += __pyx_tmp_idx * __pyx_tmp_stride;
}
__pyx_t_11.shape[0] = __pyx_v_Ds.shape[1];
__pyx_t_11.strides[0] = __pyx_v_Ds.strides[1];
__pyx_t_11.suboffsets[0] = -1;
__pyx_t_1 = __pyx_memoryview_fromslice(__pyx_t_11, 1, (PyObject *(*)(char *)) __pyx_memview_get_double, (int (*)(char *, PyObject *)) __pyx_memview_set_double, 0);; if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 31, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__PYX_XDEC_MEMVIEW(&__pyx_t_11, 1);
__pyx_t_11.memview = NULL;
__pyx_t_11.data = NULL;
__pyx_t_11.data = __pyx_v_Ds.data;
__pyx_t_11.memview = __pyx_v_Ds.memview;
__PYX_INC_MEMVIEW(&__pyx_t_11, 0);
{
Py_ssize_t __pyx_tmp_idx = 1;
Py_ssize_t __pyx_tmp_shape = __pyx_v_Ds.shape[0];
Py_ssize_t __pyx_tmp_stride = __pyx_v_Ds.strides[0];
if (0 && (__pyx_tmp_idx < 0))
__pyx_tmp_idx += __pyx_tmp_shape;
if (0 && (__pyx_tmp_idx < 0 || __pyx_tmp_idx >= __pyx_tmp_shape)) {
PyErr_SetString(PyExc_IndexError, "Index out of bounds (axis 0)");
__PYX_ERR(0, 31, __pyx_L1_error)
}
__pyx_t_11.data += __pyx_tmp_idx * __pyx_tmp_stride;
}
__pyx_t_11.shape[0] = __pyx_v_Ds.shape[1];
__pyx_t_11.strides[0] = __pyx_v_Ds.strides[1];
__pyx_t_11.suboffsets[0] = -1;
__pyx_t_4 = __pyx_memoryview_fromslice(__pyx_t_11, 1, (PyObject *(*)(char *)) __pyx_memview_get_double, (int (*)(char *, PyObject *)) __pyx_memview_set_double, 0);; if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 31, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__PYX_XDEC_MEMVIEW(&__pyx_t_11, 1);
__pyx_t_11.memview = NULL;
__pyx_t_11.data = NULL;
__pyx_t_14 = NULL;
__pyx_t_15 = 0;
if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_5))) {
__pyx_t_14 = PyMethod_GET_SELF(__pyx_t_5);
if (likely(__pyx_t_14)) {
PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_5);
__Pyx_INCREF(__pyx_t_14);
__Pyx_INCREF(function);
__Pyx_DECREF_SET(__pyx_t_5, function);
__pyx_t_15 = 1;
}
}
#if CYTHON_FAST_PYCALL
if (PyFunction_Check(__pyx_t_5)) {
PyObject *__pyx_temp[3] = {__pyx_t_14, __pyx_t_1, __pyx_t_4};
__pyx_t_3 = __Pyx_PyFunction_FastCall(__pyx_t_5, __pyx_temp+1-__pyx_t_15, 2+__pyx_t_15); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 31, __pyx_L1_error)
__Pyx_XDECREF(__pyx_t_14); __pyx_t_14 = 0;
__Pyx_GOTREF(__pyx_t_3);
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
} else
#endif
#if CYTHON_FAST_PYCCALL
if (__Pyx_PyFastCFunction_Check(__pyx_t_5)) {
PyObject *__pyx_temp[3] = {__pyx_t_14, __pyx_t_1, __pyx_t_4};
__pyx_t_3 = __Pyx_PyCFunction_FastCall(__pyx_t_5, __pyx_temp+1-__pyx_t_15, 2+__pyx_t_15); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 31, __pyx_L1_error)
__Pyx_XDECREF(__pyx_t_14); __pyx_t_14 = 0;
__Pyx_GOTREF(__pyx_t_3);
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
} else
#endif
{
__pyx_t_16 = PyTuple_New(2+__pyx_t_15); if (unlikely(!__pyx_t_16)) __PYX_ERR(0, 31, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_16);
if (__pyx_t_14) {
__Pyx_GIVEREF(__pyx_t_14); PyTuple_SET_ITEM(__pyx_t_16, 0, __pyx_t_14); __pyx_t_14 = NULL;
}
__Pyx_GIVEREF(__pyx_t_1);
PyTuple_SET_ITEM(__pyx_t_16, 0+__pyx_t_15, __pyx_t_1);
__Pyx_GIVEREF(__pyx_t_4);
PyTuple_SET_ITEM(__pyx_t_16, 1+__pyx_t_15, __pyx_t_4);
__pyx_t_1 = 0;
__pyx_t_4 = 0;
__pyx_t_3 = __Pyx_PyObject_Call(__pyx_t_5, __pyx_t_16, NULL); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 31, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__Pyx_DECREF(__pyx_t_16); __pyx_t_16 = 0;
}
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
__pyx_t_5 = NULL;
if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_6))) {
__pyx_t_5 = PyMethod_GET_SELF(__pyx_t_6);
if (likely(__pyx_t_5)) {
PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_6);
__Pyx_INCREF(__pyx_t_5);
__Pyx_INCREF(function);
__Pyx_DECREF_SET(__pyx_t_6, function);
}
}
if (!__pyx_t_5) {
__pyx_t_2 = __Pyx_PyObject_CallOneArg(__pyx_t_6, __pyx_t_3); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 31, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__Pyx_GOTREF(__pyx_t_2);
} else {
#if CYTHON_FAST_PYCALL
if (PyFunction_Check(__pyx_t_6)) {
PyObject *__pyx_temp[2] = {__pyx_t_5, __pyx_t_3};
__pyx_t_2 = __Pyx_PyFunction_FastCall(__pyx_t_6, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 31, __pyx_L1_error)
__Pyx_XDECREF(__pyx_t_5); __pyx_t_5 = 0;
__Pyx_GOTREF(__pyx_t_2);
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
} else
#endif
#if CYTHON_FAST_PYCCALL
if (__Pyx_PyFastCFunction_Check(__pyx_t_6)) {
PyObject *__pyx_temp[2] = {__pyx_t_5, __pyx_t_3};
__pyx_t_2 = __Pyx_PyCFunction_FastCall(__pyx_t_6, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 31, __pyx_L1_error)
__Pyx_XDECREF(__pyx_t_5); __pyx_t_5 = 0;
__Pyx_GOTREF(__pyx_t_2);
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
} else
#endif
{
__pyx_t_16 = PyTuple_New(1+1); if (unlikely(!__pyx_t_16)) __PYX_ERR(0, 31, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_16);
__Pyx_GIVEREF(__pyx_t_5); PyTuple_SET_ITEM(__pyx_t_16, 0, __pyx_t_5); __pyx_t_5 = NULL;
__Pyx_GIVEREF(__pyx_t_3);
PyTuple_SET_ITEM(__pyx_t_16, 0+1, __pyx_t_3);
__pyx_t_3 = 0;
__pyx_t_2 = __Pyx_PyObject_Call(__pyx_t_6, __pyx_t_16, NULL); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 31, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__Pyx_DECREF(__pyx_t_16); __pyx_t_16 = 0;
}
}
__Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
__pyx_t_16 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_16)) __PYX_ERR(0, 31, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_16);
__pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_16, __pyx_n_s_min); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 31, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__Pyx_DECREF(__pyx_t_16); __pyx_t_16 = 0;
__pyx_t_11.data = __pyx_v_VPoly.data;
__pyx_t_11.memview = __pyx_v_VPoly.memview;
__PYX_INC_MEMVIEW(&__pyx_t_11, 0);
{
Py_ssize_t __pyx_tmp_idx = 0;
Py_ssize_t __pyx_tmp_shape = __pyx_v_VPoly.shape[0];
Py_ssize_t __pyx_tmp_stride = __pyx_v_VPoly.strides[0];
if (0 && (__pyx_tmp_idx < 0))
__pyx_tmp_idx += __pyx_tmp_shape;
if (0 && (__pyx_tmp_idx < 0 || __pyx_tmp_idx >= __pyx_tmp_shape)) {
PyErr_SetString(PyExc_IndexError, "Index out of bounds (axis 0)");
__PYX_ERR(0, 31, __pyx_L1_error)
}
__pyx_t_11.data += __pyx_tmp_idx * __pyx_tmp_stride;
}
__pyx_t_11.shape[0] = __pyx_v_VPoly.shape[1];
__pyx_t_11.strides[0] = __pyx_v_VPoly.strides[1];
__pyx_t_11.suboffsets[0] = -1;
__pyx_t_16 = __pyx_memoryview_fromslice(__pyx_t_11, 1, (PyObject *(*)(char *)) __pyx_memview_get_double, (int (*)(char *, PyObject *)) __pyx_memview_set_double, 0);; if (unlikely(!__pyx_t_16)) __PYX_ERR(0, 31, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_16);
__PYX_XDEC_MEMVIEW(&__pyx_t_11, 1);
__pyx_t_11.memview = NULL;
__pyx_t_11.data = NULL;
__pyx_t_5 = NULL;
if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_3))) {
__pyx_t_5 = PyMethod_GET_SELF(__pyx_t_3);
if (likely(__pyx_t_5)) {
PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_3);
__Pyx_INCREF(__pyx_t_5);
__Pyx_INCREF(function);
__Pyx_DECREF_SET(__pyx_t_3, function);
}
}
if (!__pyx_t_5) {
__pyx_t_6 = __Pyx_PyObject_CallOneArg(__pyx_t_3, __pyx_t_16); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 31, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_16); __pyx_t_16 = 0;
__Pyx_GOTREF(__pyx_t_6);
} else {
#if CYTHON_FAST_PYCALL
if (PyFunction_Check(__pyx_t_3)) {
PyObject *__pyx_temp[2] = {__pyx_t_5, __pyx_t_16};
__pyx_t_6 = __Pyx_PyFunction_FastCall(__pyx_t_3, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 31, __pyx_L1_error)
__Pyx_XDECREF(__pyx_t_5); __pyx_t_5 = 0;
__Pyx_GOTREF(__pyx_t_6);
__Pyx_DECREF(__pyx_t_16); __pyx_t_16 = 0;
} else
#endif
#if CYTHON_FAST_PYCCALL
if (__Pyx_PyFastCFunction_Check(__pyx_t_3)) {
PyObject *__pyx_temp[2] = {__pyx_t_5, __pyx_t_16};
__pyx_t_6 = __Pyx_PyCFunction_FastCall(__pyx_t_3, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 31, __pyx_L1_error)
__Pyx_XDECREF(__pyx_t_5); __pyx_t_5 = 0;
__Pyx_GOTREF(__pyx_t_6);
__Pyx_DECREF(__pyx_t_16); __pyx_t_16 = 0;
} else
#endif
{
__pyx_t_4 = PyTuple_New(1+1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 31, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__Pyx_GIVEREF(__pyx_t_5); PyTuple_SET_ITEM(__pyx_t_4, 0, __pyx_t_5); __pyx_t_5 = NULL;
__Pyx_GIVEREF(__pyx_t_16);
PyTuple_SET_ITEM(__pyx_t_4, 0+1, __pyx_t_16);
__pyx_t_16 = 0;
__pyx_t_6 = __Pyx_PyObject_Call(__pyx_t_3, __pyx_t_4, NULL); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 31, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_6);
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
}
}
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__pyx_t_4 = PyObject_RichCompare(__pyx_t_2, __pyx_t_6, Py_LT); __Pyx_XGOTREF(__pyx_t_4); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 31, __pyx_L1_error)
__pyx_t_13 = __Pyx_PyObject_IsTrue(__pyx_t_4); if (unlikely(__pyx_t_13 < 0)) __PYX_ERR(0, 31, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
if (__pyx_t_13) {
__Pyx_INCREF(__pyx_t_2);
__pyx_t_3 = __pyx_t_2;
} else {
__Pyx_INCREF(__pyx_t_6);
__pyx_t_3 = __pyx_t_6;
}
__Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
__pyx_t_2 = PyNumber_Multiply(__pyx_float_0_95, __pyx_t_3); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 31, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__pyx_t_17 = __pyx_PyFloat_AsDouble(__pyx_t_2); if (unlikely((__pyx_t_17 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 31, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
__pyx_v_Rmin = __pyx_t_17;
032: else:
+033: Rmin = RMin
/*else*/ {
__pyx_t_17 = __pyx_PyFloat_AsDouble(__pyx_v_RMin); if (unlikely((__pyx_t_17 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 33, __pyx_L1_error)
__pyx_v_Rmin = __pyx_t_17;
}
__pyx_L3:;
034:
035: ################
036: # Compute
+037: if Forbid:
__pyx_t_13 = __Pyx_PyObject_IsTrue(((PyObject *)__pyx_v_Forbid)); if (unlikely(__pyx_t_13 < 0)) __PYX_ERR(0, 37, __pyx_L1_error)
if (__pyx_t_13) {
/* … */
goto __pyx_L4;
}
+038: Forbid0, Forbidbis = 1, 1
__pyx_t_15 = 1;
__pyx_t_18 = 1;
__pyx_v_Forbid0 = __pyx_t_15;
__pyx_v_Forbidbis = __pyx_t_18;
039: else:
+040: Forbid0, Forbidbis = 0, 0
/*else*/ {
__pyx_t_18 = 0;
__pyx_t_15 = 0;
__pyx_v_Forbid0 = __pyx_t_18;
__pyx_v_Forbidbis = __pyx_t_15;
}
__pyx_L4:;
+041: for ii in range(0,Nl):
__pyx_t_15 = __pyx_v_Nl;
for (__pyx_t_18 = 0; __pyx_t_18 < __pyx_t_15; __pyx_t_18+=1) {
__pyx_v_ii = __pyx_t_18;
+042: upscaDp = us[0,ii]*Ds[0,ii] + us[1,ii]*Ds[1,ii]
__pyx_t_19 = 0;
__pyx_t_20 = __pyx_v_ii;
__pyx_t_21 = 0;
__pyx_t_22 = __pyx_v_ii;
__pyx_t_23 = 1;
__pyx_t_24 = __pyx_v_ii;
__pyx_t_25 = 1;
__pyx_t_26 = __pyx_v_ii;
__pyx_v_upscaDp = (((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_19 * __pyx_v_us.strides[0]) )) + __pyx_t_20)) ))) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_21 * __pyx_v_Ds.strides[0]) )) + __pyx_t_22)) )))) + ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_23 * __pyx_v_us.strides[0]) )) + __pyx_t_24)) ))) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_25 * __pyx_v_Ds.strides[0]) )) + __pyx_t_26)) )))));
+043: upar2 = us[0,ii]**2 + us[1,ii]**2
__pyx_t_27 = 0;
__pyx_t_28 = __pyx_v_ii;
__pyx_t_29 = 1;
__pyx_t_30 = __pyx_v_ii;
__pyx_v_upar2 = (pow((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_27 * __pyx_v_us.strides[0]) )) + __pyx_t_28)) ))), 2.0) + pow((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_29 * __pyx_v_us.strides[0]) )) + __pyx_t_30)) ))), 2.0));
+044: Dpar2 = Ds[0,ii]**2 + Ds[1,ii]**2
__pyx_t_31 = 0;
__pyx_t_32 = __pyx_v_ii;
__pyx_t_33 = 1;
__pyx_t_34 = __pyx_v_ii;
__pyx_v_Dpar2 = (pow((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_31 * __pyx_v_Ds.strides[0]) )) + __pyx_t_32)) ))), 2.0) + pow((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_33 * __pyx_v_Ds.strides[0]) )) + __pyx_t_34)) ))), 2.0));
045: # Prepare in case Forbid is True
+046: if Forbid0 and not Dpar2>0:
__pyx_t_12 = (__pyx_v_Forbid0 != 0);
if (__pyx_t_12) {
} else {
__pyx_t_13 = __pyx_t_12;
goto __pyx_L8_bool_binop_done;
}
__pyx_t_12 = ((!((__pyx_v_Dpar2 > 0.0) != 0)) != 0);
__pyx_t_13 = __pyx_t_12;
__pyx_L8_bool_binop_done:;
if (__pyx_t_13) {
/* … */
}
+047: Forbidbis = 0
__pyx_v_Forbidbis = 0;
+048: if Forbidbis:
__pyx_t_13 = (__pyx_v_Forbidbis != 0);
if (__pyx_t_13) {
/* … */
}
049: # Compute coordinates of the 2 points where the tangents touch the inner circle
+050: L = Csqrt(Dpar2-Rmin**2)
__pyx_v_L = sqrt((__pyx_v_Dpar2 - pow(__pyx_v_Rmin, 2.0)));
+051: S1X = (Rmin**2*Ds[0,ii]+Rmin*Ds[1,ii]*L)/Dpar2
__pyx_t_35 = 0;
__pyx_t_36 = __pyx_v_ii;
__pyx_t_37 = 1;
__pyx_t_38 = __pyx_v_ii;
__pyx_v_S1X = (((pow(__pyx_v_Rmin, 2.0) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_35 * __pyx_v_Ds.strides[0]) )) + __pyx_t_36)) )))) + ((__pyx_v_Rmin * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_37 * __pyx_v_Ds.strides[0]) )) + __pyx_t_38)) )))) * __pyx_v_L)) / __pyx_v_Dpar2);
+052: S1Y = (Rmin**2*Ds[1,ii]-Rmin*Ds[0,ii]*L)/Dpar2
__pyx_t_39 = 1;
__pyx_t_40 = __pyx_v_ii;
__pyx_t_41 = 0;
__pyx_t_42 = __pyx_v_ii;
__pyx_v_S1Y = (((pow(__pyx_v_Rmin, 2.0) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_39 * __pyx_v_Ds.strides[0]) )) + __pyx_t_40)) )))) - ((__pyx_v_Rmin * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_41 * __pyx_v_Ds.strides[0]) )) + __pyx_t_42)) )))) * __pyx_v_L)) / __pyx_v_Dpar2);
+053: S2X = (Rmin**2*Ds[0,ii]-Rmin*Ds[1,ii]*L)/Dpar2
__pyx_t_43 = 0;
__pyx_t_44 = __pyx_v_ii;
__pyx_t_45 = 1;
__pyx_t_46 = __pyx_v_ii;
__pyx_v_S2X = (((pow(__pyx_v_Rmin, 2.0) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_43 * __pyx_v_Ds.strides[0]) )) + __pyx_t_44)) )))) - ((__pyx_v_Rmin * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_45 * __pyx_v_Ds.strides[0]) )) + __pyx_t_46)) )))) * __pyx_v_L)) / __pyx_v_Dpar2);
+054: S2Y = (Rmin**2*Ds[1,ii]+Rmin*Ds[0,ii]*L)/Dpar2
__pyx_t_47 = 1;
__pyx_t_48 = __pyx_v_ii;
__pyx_t_49 = 0;
__pyx_t_50 = __pyx_v_ii;
__pyx_v_S2Y = (((pow(__pyx_v_Rmin, 2.0) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_47 * __pyx_v_Ds.strides[0]) )) + __pyx_t_48)) )))) + ((__pyx_v_Rmin * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_49 * __pyx_v_Ds.strides[0]) )) + __pyx_t_50)) )))) * __pyx_v_L)) / __pyx_v_Dpar2);
055:
056: # Compute all solutions
057: # Set tolerance value for us[2,ii]
058: # EpsUz is the tolerated DZ across 20m (max Tokamak size)
+059: Crit2 = EpsUz**2*upar2/400.
__pyx_v_Crit2 = ((pow(__pyx_v_EpsUz, 2.0) * __pyx_v_upar2) / 400.);
+060: kout, kin, Done = 1.e12, 1e12, 0
__pyx_t_17 = 1.e12;
__pyx_t_51 = 1e12;
__pyx_t_52 = 0;
__pyx_v_kout = __pyx_t_17;
__pyx_v_kin = __pyx_t_51;
__pyx_v_Done = __pyx_t_52;
061: # Case with horizontal semi-line
+062: if us[2,ii]**2<Crit2:
__pyx_t_53 = 2;
__pyx_t_54 = __pyx_v_ii;
__pyx_t_13 = ((pow((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_53 * __pyx_v_us.strides[0]) )) + __pyx_t_54)) ))), 2.0) < __pyx_v_Crit2) != 0);
if (__pyx_t_13) {
/* … */
goto __pyx_L11;
}
+063: for jj in range(0,Ns):
__pyx_t_52 = __pyx_v_Ns;
for (__pyx_t_55 = 0; __pyx_t_55 < __pyx_t_52; __pyx_t_55+=1) {
__pyx_v_jj = __pyx_t_55;
064: # Solutions exist only in the case with non-horizontal segment (i.e.: cone, not plane)
+065: if (VPoly[1,jj+1]-VPoly[1,jj])**2>EpsVz**2:
__pyx_t_56 = 1;
__pyx_t_57 = (__pyx_v_jj + 1);
__pyx_t_58 = 1;
__pyx_t_59 = __pyx_v_jj;
__pyx_t_13 = ((pow(((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_VPoly.data + __pyx_t_56 * __pyx_v_VPoly.strides[0]) )) + __pyx_t_57)) ))) - (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_VPoly.data + __pyx_t_58 * __pyx_v_VPoly.strides[0]) )) + __pyx_t_59)) )))), 2.0) > pow(__pyx_v_EpsVz, 2.0)) != 0);
if (__pyx_t_13) {
/* … */
}
}
+066: q = (Ds[2,ii]-VPoly[1,jj])/(VPoly[1,jj+1]-VPoly[1,jj])
__pyx_t_60 = 2;
__pyx_t_61 = __pyx_v_ii;
__pyx_t_62 = 1;
__pyx_t_63 = __pyx_v_jj;
__pyx_t_64 = 1;
__pyx_t_65 = (__pyx_v_jj + 1);
__pyx_t_66 = 1;
__pyx_t_67 = __pyx_v_jj;
__pyx_v_q = (((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_60 * __pyx_v_Ds.strides[0]) )) + __pyx_t_61)) ))) - (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_VPoly.data + __pyx_t_62 * __pyx_v_VPoly.strides[0]) )) + __pyx_t_63)) )))) / ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_VPoly.data + __pyx_t_64 * __pyx_v_VPoly.strides[0]) )) + __pyx_t_65)) ))) - (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_VPoly.data + __pyx_t_66 * __pyx_v_VPoly.strides[0]) )) + __pyx_t_67)) )))));
067: # The intersection must stand on the segment
+068: if q>=0 and q<1:
__pyx_t_12 = ((__pyx_v_q >= 0.0) != 0);
if (__pyx_t_12) {
} else {
__pyx_t_13 = __pyx_t_12;
goto __pyx_L16_bool_binop_done;
}
__pyx_t_12 = ((__pyx_v_q < 1.0) != 0);
__pyx_t_13 = __pyx_t_12;
__pyx_L16_bool_binop_done:;
if (__pyx_t_13) {
/* … */
}
+069: C = q**2*(VPoly[0,jj+1]-VPoly[0,jj])**2 + 2.*q*VPoly[0,jj]*(VPoly[0,jj+1]-VPoly[0,jj]) + VPoly[0,jj]**2
__pyx_t_68 = 0;
__pyx_t_69 = (__pyx_v_jj + 1);
__pyx_t_70 = 0;
__pyx_t_71 = __pyx_v_jj;
__pyx_t_72 = 0;
__pyx_t_73 = __pyx_v_jj;
__pyx_t_74 = 0;
__pyx_t_75 = (__pyx_v_jj + 1);
__pyx_t_76 = 0;
__pyx_t_77 = __pyx_v_jj;
__pyx_t_78 = 0;
__pyx_t_79 = __pyx_v_jj;
__pyx_v_C = (((pow(__pyx_v_q, 2.0) * pow(((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_VPoly.data + __pyx_t_68 * __pyx_v_VPoly.strides[0]) )) + __pyx_t_69)) ))) - (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_VPoly.data + __pyx_t_70 * __pyx_v_VPoly.strides[0]) )) + __pyx_t_71)) )))), 2.0)) + (((2. * __pyx_v_q) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_VPoly.data + __pyx_t_72 * __pyx_v_VPoly.strides[0]) )) + __pyx_t_73)) )))) * ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_VPoly.data + __pyx_t_74 * __pyx_v_VPoly.strides[0]) )) + __pyx_t_75)) ))) - (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_VPoly.data + __pyx_t_76 * __pyx_v_VPoly.strides[0]) )) + __pyx_t_77)) )))))) + pow((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_VPoly.data + __pyx_t_78 * __pyx_v_VPoly.strides[0]) )) + __pyx_t_79)) ))), 2.0));
+070: delta = upscaDp**2 - upar2*(Dpar2-C)
__pyx_v_delta = (pow(__pyx_v_upscaDp, 2.0) - (__pyx_v_upar2 * (__pyx_v_Dpar2 - __pyx_v_C)));
+071: if delta>0.:
__pyx_t_13 = ((__pyx_v_delta > 0.) != 0);
if (__pyx_t_13) {
/* … */
}
+072: sqd = Csqrt(delta)
__pyx_v_sqd = sqrt(__pyx_v_delta);
073: # The intersection must be on the semi-line (i.e.: k>=0)
074: # First solution
+075: if -upscaDp - sqd >=0:
__pyx_t_13 = ((((-__pyx_v_upscaDp) - __pyx_v_sqd) >= 0.0) != 0);
if (__pyx_t_13) {
/* … */
}
+076: k = (-upscaDp - sqd)/upar2
__pyx_v_k = (((-__pyx_v_upscaDp) - __pyx_v_sqd) / __pyx_v_upar2);
+077: sol0, sol1 = Ds[0,ii] + k*us[0,ii], Ds[1,ii] + k*us[1,ii]
__pyx_t_80 = 0;
__pyx_t_81 = __pyx_v_ii;
__pyx_t_82 = 0;
__pyx_t_83 = __pyx_v_ii;
__pyx_t_51 = ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_80 * __pyx_v_Ds.strides[0]) )) + __pyx_t_81)) ))) + (__pyx_v_k * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_82 * __pyx_v_us.strides[0]) )) + __pyx_t_83)) )))));
__pyx_t_84 = 1;
__pyx_t_85 = __pyx_v_ii;
__pyx_t_86 = 1;
__pyx_t_87 = __pyx_v_ii;
__pyx_t_17 = ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_84 * __pyx_v_Ds.strides[0]) )) + __pyx_t_85)) ))) + (__pyx_v_k * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_86 * __pyx_v_us.strides[0]) )) + __pyx_t_87)) )))));
__pyx_v_sol0 = __pyx_t_51;
__pyx_v_sol1 = __pyx_t_17;
+078: if Forbidbis:
__pyx_t_13 = (__pyx_v_Forbidbis != 0);
if (__pyx_t_13) {
/* … */
}
+079: sca0 = (sol0-S1X)*Ds[0,ii] + (sol1-S1Y)*Ds[1,ii]
__pyx_t_88 = 0;
__pyx_t_89 = __pyx_v_ii;
__pyx_t_90 = 1;
__pyx_t_91 = __pyx_v_ii;
__pyx_v_sca0 = (((__pyx_v_sol0 - __pyx_v_S1X) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_88 * __pyx_v_Ds.strides[0]) )) + __pyx_t_89)) )))) + ((__pyx_v_sol1 - __pyx_v_S1Y) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_90 * __pyx_v_Ds.strides[0]) )) + __pyx_t_91)) )))));
+080: sca1 = (sol0-S1X)*S1X + (sol1-S1Y)*S1Y
__pyx_v_sca1 = (((__pyx_v_sol0 - __pyx_v_S1X) * __pyx_v_S1X) + ((__pyx_v_sol1 - __pyx_v_S1Y) * __pyx_v_S1Y));
+081: sca2 = (sol0-S2X)*S2X + (sol1-S2Y)*S2Y
__pyx_v_sca2 = (((__pyx_v_sol0 - __pyx_v_S2X) * __pyx_v_S2X) + ((__pyx_v_sol1 - __pyx_v_S2Y) * __pyx_v_S2Y));
+082: if not Forbidbis or (Forbidbis and not (sca0<0 and sca1<0 and sca2<0)):
__pyx_t_12 = ((!(__pyx_v_Forbidbis != 0)) != 0);
if (!__pyx_t_12) {
} else {
__pyx_t_13 = __pyx_t_12;
goto __pyx_L22_bool_binop_done;
}
__pyx_t_12 = (__pyx_v_Forbidbis != 0);
if (__pyx_t_12) {
} else {
__pyx_t_13 = __pyx_t_12;
goto __pyx_L22_bool_binop_done;
}
__pyx_t_92 = ((__pyx_v_sca0 < 0.0) != 0);
if (__pyx_t_92) {
} else {
__pyx_t_12 = __pyx_t_92;
goto __pyx_L25_bool_binop_done;
}
__pyx_t_92 = ((__pyx_v_sca1 < 0.0) != 0);
if (__pyx_t_92) {
} else {
__pyx_t_12 = __pyx_t_92;
goto __pyx_L25_bool_binop_done;
}
__pyx_t_92 = ((__pyx_v_sca2 < 0.0) != 0);
__pyx_t_12 = __pyx_t_92;
__pyx_L25_bool_binop_done:;
__pyx_t_92 = ((!__pyx_t_12) != 0);
__pyx_t_13 = __pyx_t_92;
__pyx_L22_bool_binop_done:;
if (__pyx_t_13) {
/* … */
}
083: # Get the normalized perpendicular vector at intersection
+084: phi = Catan2(sol1,sol0)
__pyx_v_phi = atan2(__pyx_v_sol1, __pyx_v_sol0);
085: # Get the scalar product to determine entry or exit point
+086: sca = Ccos(phi)*vIn[0,jj]*us[0,ii] + Csin(phi)*vIn[0,jj]*us[1,ii] + vIn[1,jj]*us[2,ii]
__pyx_t_93 = 0;
__pyx_t_94 = __pyx_v_jj;
__pyx_t_95 = 0;
__pyx_t_96 = __pyx_v_ii;
__pyx_t_97 = 0;
__pyx_t_98 = __pyx_v_jj;
__pyx_t_99 = 1;
__pyx_t_100 = __pyx_v_ii;
__pyx_t_101 = 1;
__pyx_t_102 = __pyx_v_jj;
__pyx_t_103 = 2;
__pyx_t_104 = __pyx_v_ii;
__pyx_v_sca = ((((cos(__pyx_v_phi) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_vIn.data + __pyx_t_93 * __pyx_v_vIn.strides[0]) )) + __pyx_t_94)) )))) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_95 * __pyx_v_us.strides[0]) )) + __pyx_t_96)) )))) + ((sin(__pyx_v_phi) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_vIn.data + __pyx_t_97 * __pyx_v_vIn.strides[0]) )) + __pyx_t_98)) )))) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_99 * __pyx_v_us.strides[0]) )) + __pyx_t_100)) ))))) + ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_vIn.data + __pyx_t_101 * __pyx_v_vIn.strides[0]) )) + __pyx_t_102)) ))) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_103 * __pyx_v_us.strides[0]) )) + __pyx_t_104)) )))));
+087: if sca<=0 and k<kout:
__pyx_t_92 = ((__pyx_v_sca <= 0.0) != 0);
if (__pyx_t_92) {
} else {
__pyx_t_13 = __pyx_t_92;
goto __pyx_L29_bool_binop_done;
}
__pyx_t_92 = ((__pyx_v_k < __pyx_v_kout) != 0);
__pyx_t_13 = __pyx_t_92;
__pyx_L29_bool_binop_done:;
if (__pyx_t_13) {
/* … */
goto __pyx_L28;
}
+088: kout = k
__pyx_v_kout = __pyx_v_k;
+089: indout = jj
__pyx_v_indout = __pyx_v_jj;
+090: Done = 1
__pyx_v_Done = 1;
+091: print 1, k
__pyx_t_2 = PyFloat_FromDouble(__pyx_v_k); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 91, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__pyx_t_3 = PyTuple_New(2); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 91, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__Pyx_INCREF(__pyx_int_1);
__Pyx_GIVEREF(__pyx_int_1);
PyTuple_SET_ITEM(__pyx_t_3, 0, __pyx_int_1);
__Pyx_GIVEREF(__pyx_t_2);
PyTuple_SET_ITEM(__pyx_t_3, 1, __pyx_t_2);
__pyx_t_2 = 0;
if (__Pyx_Print(0, __pyx_t_3, 1) < 0) __PYX_ERR(0, 91, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
+092: elif sca>=0 and k<min(kin,kout):
__pyx_t_92 = ((__pyx_v_sca >= 0.0) != 0);
if (__pyx_t_92) {
} else {
__pyx_t_13 = __pyx_t_92;
goto __pyx_L31_bool_binop_done;
}
__pyx_t_17 = __pyx_v_kout;
__pyx_t_51 = __pyx_v_kin;
if (((__pyx_t_17 < __pyx_t_51) != 0)) {
__pyx_t_105 = __pyx_t_17;
} else {
__pyx_t_105 = __pyx_t_51;
}
__pyx_t_92 = ((__pyx_v_k < __pyx_t_105) != 0);
__pyx_t_13 = __pyx_t_92;
__pyx_L31_bool_binop_done:;
if (__pyx_t_13) {
/* … */
}
__pyx_L28:;
+093: kin = k
__pyx_v_kin = __pyx_v_k;
+094: print 2, k
__pyx_t_3 = PyFloat_FromDouble(__pyx_v_k); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 94, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__pyx_t_2 = PyTuple_New(2); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 94, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__Pyx_INCREF(__pyx_int_2);
__Pyx_GIVEREF(__pyx_int_2);
PyTuple_SET_ITEM(__pyx_t_2, 0, __pyx_int_2);
__Pyx_GIVEREF(__pyx_t_3);
PyTuple_SET_ITEM(__pyx_t_2, 1, __pyx_t_3);
__pyx_t_3 = 0;
if (__Pyx_Print(0, __pyx_t_2, 1) < 0) __PYX_ERR(0, 94, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
095:
096: # Second solution
+097: if -upscaDp + sqd >=0:
__pyx_t_13 = ((((-__pyx_v_upscaDp) + __pyx_v_sqd) >= 0.0) != 0);
if (__pyx_t_13) {
/* … */
}
+098: k = (-upscaDp + sqd)/upar2
__pyx_v_k = (((-__pyx_v_upscaDp) + __pyx_v_sqd) / __pyx_v_upar2);
+099: sol0, sol1 = Ds[0,ii] + k*us[0,ii], Ds[1,ii] + k*us[1,ii]
__pyx_t_106 = 0;
__pyx_t_107 = __pyx_v_ii;
__pyx_t_108 = 0;
__pyx_t_109 = __pyx_v_ii;
__pyx_t_105 = ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_106 * __pyx_v_Ds.strides[0]) )) + __pyx_t_107)) ))) + (__pyx_v_k * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_108 * __pyx_v_us.strides[0]) )) + __pyx_t_109)) )))));
__pyx_t_110 = 1;
__pyx_t_111 = __pyx_v_ii;
__pyx_t_112 = 1;
__pyx_t_113 = __pyx_v_ii;
__pyx_t_17 = ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_110 * __pyx_v_Ds.strides[0]) )) + __pyx_t_111)) ))) + (__pyx_v_k * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_112 * __pyx_v_us.strides[0]) )) + __pyx_t_113)) )))));
__pyx_v_sol0 = __pyx_t_105;
__pyx_v_sol1 = __pyx_t_17;
+100: if Forbidbis:
__pyx_t_13 = (__pyx_v_Forbidbis != 0);
if (__pyx_t_13) {
/* … */
}
+101: sca0 = (sol0-S1X)*Ds[0,ii] + (sol1-S1Y)*Ds[1,ii]
__pyx_t_114 = 0;
__pyx_t_115 = __pyx_v_ii;
__pyx_t_116 = 1;
__pyx_t_117 = __pyx_v_ii;
__pyx_v_sca0 = (((__pyx_v_sol0 - __pyx_v_S1X) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_114 * __pyx_v_Ds.strides[0]) )) + __pyx_t_115)) )))) + ((__pyx_v_sol1 - __pyx_v_S1Y) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_116 * __pyx_v_Ds.strides[0]) )) + __pyx_t_117)) )))));
+102: sca1 = (sol0-S1X)*S1X + (sol1-S1Y)*S1Y
__pyx_v_sca1 = (((__pyx_v_sol0 - __pyx_v_S1X) * __pyx_v_S1X) + ((__pyx_v_sol1 - __pyx_v_S1Y) * __pyx_v_S1Y));
+103: sca2 = (sol0-S2X)*S2X + (sol1-S2Y)*S2Y
__pyx_v_sca2 = (((__pyx_v_sol0 - __pyx_v_S2X) * __pyx_v_S2X) + ((__pyx_v_sol1 - __pyx_v_S2Y) * __pyx_v_S2Y));
+104: if not Forbidbis or (Forbidbis and not (sca0<0 and sca1<0 and sca2<0)):
__pyx_t_92 = ((!(__pyx_v_Forbidbis != 0)) != 0);
if (!__pyx_t_92) {
} else {
__pyx_t_13 = __pyx_t_92;
goto __pyx_L36_bool_binop_done;
}
__pyx_t_92 = (__pyx_v_Forbidbis != 0);
if (__pyx_t_92) {
} else {
__pyx_t_13 = __pyx_t_92;
goto __pyx_L36_bool_binop_done;
}
__pyx_t_12 = ((__pyx_v_sca0 < 0.0) != 0);
if (__pyx_t_12) {
} else {
__pyx_t_92 = __pyx_t_12;
goto __pyx_L39_bool_binop_done;
}
__pyx_t_12 = ((__pyx_v_sca1 < 0.0) != 0);
if (__pyx_t_12) {
} else {
__pyx_t_92 = __pyx_t_12;
goto __pyx_L39_bool_binop_done;
}
__pyx_t_12 = ((__pyx_v_sca2 < 0.0) != 0);
__pyx_t_92 = __pyx_t_12;
__pyx_L39_bool_binop_done:;
__pyx_t_12 = ((!__pyx_t_92) != 0);
__pyx_t_13 = __pyx_t_12;
__pyx_L36_bool_binop_done:;
if (__pyx_t_13) {
/* … */
}
105: # Get the normalized perpendicular vector at intersection
+106: phi = Catan2(sol1,sol0)
__pyx_v_phi = atan2(__pyx_v_sol1, __pyx_v_sol0);
107: # Get the scalar product to determine entry or exit point
+108: sca = Ccos(phi)*vIn[0,jj]*us[0,ii] + Csin(phi)*vIn[0,jj]*us[1,ii] + vIn[1,jj]*us[2,ii]
__pyx_t_118 = 0;
__pyx_t_119 = __pyx_v_jj;
__pyx_t_120 = 0;
__pyx_t_121 = __pyx_v_ii;
__pyx_t_122 = 0;
__pyx_t_123 = __pyx_v_jj;
__pyx_t_124 = 1;
__pyx_t_125 = __pyx_v_ii;
__pyx_t_126 = 1;
__pyx_t_127 = __pyx_v_jj;
__pyx_t_128 = 2;
__pyx_t_129 = __pyx_v_ii;
__pyx_v_sca = ((((cos(__pyx_v_phi) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_vIn.data + __pyx_t_118 * __pyx_v_vIn.strides[0]) )) + __pyx_t_119)) )))) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_120 * __pyx_v_us.strides[0]) )) + __pyx_t_121)) )))) + ((sin(__pyx_v_phi) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_vIn.data + __pyx_t_122 * __pyx_v_vIn.strides[0]) )) + __pyx_t_123)) )))) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_124 * __pyx_v_us.strides[0]) )) + __pyx_t_125)) ))))) + ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_vIn.data + __pyx_t_126 * __pyx_v_vIn.strides[0]) )) + __pyx_t_127)) ))) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_128 * __pyx_v_us.strides[0]) )) + __pyx_t_129)) )))));
+109: if sca<=0 and k<kout:
__pyx_t_12 = ((__pyx_v_sca <= 0.0) != 0);
if (__pyx_t_12) {
} else {
__pyx_t_13 = __pyx_t_12;
goto __pyx_L43_bool_binop_done;
}
__pyx_t_12 = ((__pyx_v_k < __pyx_v_kout) != 0);
__pyx_t_13 = __pyx_t_12;
__pyx_L43_bool_binop_done:;
if (__pyx_t_13) {
/* … */
goto __pyx_L42;
}
+110: kout = k
__pyx_v_kout = __pyx_v_k;
+111: indout = jj
__pyx_v_indout = __pyx_v_jj;
+112: Done = 1
__pyx_v_Done = 1;
+113: print 3, k
__pyx_t_2 = PyFloat_FromDouble(__pyx_v_k); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 113, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__pyx_t_3 = PyTuple_New(2); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 113, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__Pyx_INCREF(__pyx_int_3);
__Pyx_GIVEREF(__pyx_int_3);
PyTuple_SET_ITEM(__pyx_t_3, 0, __pyx_int_3);
__Pyx_GIVEREF(__pyx_t_2);
PyTuple_SET_ITEM(__pyx_t_3, 1, __pyx_t_2);
__pyx_t_2 = 0;
if (__Pyx_Print(0, __pyx_t_3, 1) < 0) __PYX_ERR(0, 113, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
+114: elif sca>=0 and k<min(kin,kout):
__pyx_t_12 = ((__pyx_v_sca >= 0.0) != 0);
if (__pyx_t_12) {
} else {
__pyx_t_13 = __pyx_t_12;
goto __pyx_L45_bool_binop_done;
}
__pyx_t_17 = __pyx_v_kout;
__pyx_t_105 = __pyx_v_kin;
if (((__pyx_t_17 < __pyx_t_105) != 0)) {
__pyx_t_51 = __pyx_t_17;
} else {
__pyx_t_51 = __pyx_t_105;
}
__pyx_t_12 = ((__pyx_v_k < __pyx_t_51) != 0);
__pyx_t_13 = __pyx_t_12;
__pyx_L45_bool_binop_done:;
if (__pyx_t_13) {
/* … */
}
__pyx_L42:;
+115: kin = k
__pyx_v_kin = __pyx_v_k;
+116: print 4, k
__pyx_t_3 = PyFloat_FromDouble(__pyx_v_k); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 116, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__pyx_t_2 = PyTuple_New(2); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 116, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__Pyx_INCREF(__pyx_int_4);
__Pyx_GIVEREF(__pyx_int_4);
PyTuple_SET_ITEM(__pyx_t_2, 0, __pyx_int_4);
__Pyx_GIVEREF(__pyx_t_3);
PyTuple_SET_ITEM(__pyx_t_2, 1, __pyx_t_3);
__pyx_t_3 = 0;
if (__Pyx_Print(0, __pyx_t_2, 1) < 0) __PYX_ERR(0, 116, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
117:
118: # More general non-horizontal semi-line case
119: else:
+120: for jj in range(Ns):
/*else*/ {
__pyx_t_52 = __pyx_v_Ns;
for (__pyx_t_55 = 0; __pyx_t_55 < __pyx_t_52; __pyx_t_55+=1) {
__pyx_v_jj = __pyx_t_55;
+121: v0, v1 = VPoly[0,jj+1]-VPoly[0,jj], VPoly[1,jj+1]-VPoly[1,jj]
__pyx_t_130 = 0;
__pyx_t_131 = (__pyx_v_jj + 1);
__pyx_t_132 = 0;
__pyx_t_133 = __pyx_v_jj;
__pyx_t_51 = ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_VPoly.data + __pyx_t_130 * __pyx_v_VPoly.strides[0]) )) + __pyx_t_131)) ))) - (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_VPoly.data + __pyx_t_132 * __pyx_v_VPoly.strides[0]) )) + __pyx_t_133)) ))));
__pyx_t_134 = 1;
__pyx_t_135 = (__pyx_v_jj + 1);
__pyx_t_136 = 1;
__pyx_t_137 = __pyx_v_jj;
__pyx_t_17 = ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_VPoly.data + __pyx_t_134 * __pyx_v_VPoly.strides[0]) )) + __pyx_t_135)) ))) - (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_VPoly.data + __pyx_t_136 * __pyx_v_VPoly.strides[0]) )) + __pyx_t_137)) ))));
__pyx_v_v0 = __pyx_t_51;
__pyx_v_v1 = __pyx_t_17;
+122: A = v0**2 - upar2*(v1/us[2,ii])**2
__pyx_t_138 = 2;
__pyx_t_139 = __pyx_v_ii;
__pyx_v_A = (pow(__pyx_v_v0, 2.0) - (__pyx_v_upar2 * pow((__pyx_v_v1 / (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_138 * __pyx_v_us.strides[0]) )) + __pyx_t_139)) )))), 2.0)));
+123: B = VPoly[0,jj]*v0 + v1*(Ds[2,ii]-VPoly[1,jj])*upar2/us[2,ii]**2 - upscaDp*v1/us[2,ii]
__pyx_t_140 = 0;
__pyx_t_141 = __pyx_v_jj;
__pyx_t_142 = 2;
__pyx_t_143 = __pyx_v_ii;
__pyx_t_144 = 1;
__pyx_t_145 = __pyx_v_jj;
__pyx_t_146 = 2;
__pyx_t_147 = __pyx_v_ii;
__pyx_t_148 = 2;
__pyx_t_149 = __pyx_v_ii;
__pyx_v_B = ((((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_VPoly.data + __pyx_t_140 * __pyx_v_VPoly.strides[0]) )) + __pyx_t_141)) ))) * __pyx_v_v0) + (((__pyx_v_v1 * ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_142 * __pyx_v_Ds.strides[0]) )) + __pyx_t_143)) ))) - (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_VPoly.data + __pyx_t_144 * __pyx_v_VPoly.strides[0]) )) + __pyx_t_145)) ))))) * __pyx_v_upar2) / pow((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_146 * __pyx_v_us.strides[0]) )) + __pyx_t_147)) ))), 2.0))) - ((__pyx_v_upscaDp * __pyx_v_v1) / (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_148 * __pyx_v_us.strides[0]) )) + __pyx_t_149)) )))));
+124: C = -upar2*(Ds[2,ii]-VPoly[1,jj])**2/us[2,ii]**2 + 2.*upscaDp*(Ds[2,ii]-VPoly[1,jj])/us[2,ii] - Dpar2 + VPoly[0,jj]**2
__pyx_t_150 = 2;
__pyx_t_151 = __pyx_v_ii;
__pyx_t_152 = 1;
__pyx_t_153 = __pyx_v_jj;
__pyx_t_154 = 2;
__pyx_t_155 = __pyx_v_ii;
__pyx_t_156 = 2;
__pyx_t_157 = __pyx_v_ii;
__pyx_t_158 = 1;
__pyx_t_159 = __pyx_v_jj;
__pyx_t_160 = 2;
__pyx_t_161 = __pyx_v_ii;
__pyx_t_162 = 0;
__pyx_t_163 = __pyx_v_jj;
__pyx_v_C = ((((((-__pyx_v_upar2) * pow(((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_150 * __pyx_v_Ds.strides[0]) )) + __pyx_t_151)) ))) - (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_VPoly.data + __pyx_t_152 * __pyx_v_VPoly.strides[0]) )) + __pyx_t_153)) )))), 2.0)) / pow((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_154 * __pyx_v_us.strides[0]) )) + __pyx_t_155)) ))), 2.0)) + (((2. * __pyx_v_upscaDp) * ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_156 * __pyx_v_Ds.strides[0]) )) + __pyx_t_157)) ))) - (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_VPoly.data + __pyx_t_158 * __pyx_v_VPoly.strides[0]) )) + __pyx_t_159)) ))))) / (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_160 * __pyx_v_us.strides[0]) )) + __pyx_t_161)) ))))) - __pyx_v_Dpar2) + pow((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_VPoly.data + __pyx_t_162 * __pyx_v_VPoly.strides[0]) )) + __pyx_t_163)) ))), 2.0));
125:
+126: if A**2<EpsA**2 and B**2>EpsB**2:
__pyx_t_12 = ((pow(__pyx_v_A, 2.0) < pow(__pyx_v_EpsA, 2.0)) != 0);
if (__pyx_t_12) {
} else {
__pyx_t_13 = __pyx_t_12;
goto __pyx_L50_bool_binop_done;
}
__pyx_t_12 = ((pow(__pyx_v_B, 2.0) > pow(__pyx_v_EpsB, 2.0)) != 0);
__pyx_t_13 = __pyx_t_12;
__pyx_L50_bool_binop_done:;
if (__pyx_t_13) {
/* … */
goto __pyx_L49;
}
+127: q = -C/(2.*B)
__pyx_v_q = ((-__pyx_v_C) / (2. * __pyx_v_B));
+128: if q>=0. and q<1.:
__pyx_t_12 = ((__pyx_v_q >= 0.) != 0);
if (__pyx_t_12) {
} else {
__pyx_t_13 = __pyx_t_12;
goto __pyx_L53_bool_binop_done;
}
__pyx_t_12 = ((__pyx_v_q < 1.) != 0);
__pyx_t_13 = __pyx_t_12;
__pyx_L53_bool_binop_done:;
if (__pyx_t_13) {
/* … */
}
+129: k = (q*v1 - (Ds[2,ii]-VPoly[2,jj]))/us[2,ii]
__pyx_t_164 = 2;
__pyx_t_165 = __pyx_v_ii;
__pyx_t_166 = 2;
__pyx_t_167 = __pyx_v_jj;
__pyx_t_168 = 2;
__pyx_t_169 = __pyx_v_ii;
__pyx_v_k = (((__pyx_v_q * __pyx_v_v1) - ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_164 * __pyx_v_Ds.strides[0]) )) + __pyx_t_165)) ))) - (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_VPoly.data + __pyx_t_166 * __pyx_v_VPoly.strides[0]) )) + __pyx_t_167)) ))))) / (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_168 * __pyx_v_us.strides[0]) )) + __pyx_t_169)) ))));
+130: if k>=0:
__pyx_t_13 = ((__pyx_v_k >= 0.0) != 0);
if (__pyx_t_13) {
/* … */
}
+131: sol0, sol1 = Ds[0,ii] + k*us[0,ii], Ds[1,ii] + k*us[1,ii]
__pyx_t_170 = 0;
__pyx_t_171 = __pyx_v_ii;
__pyx_t_172 = 0;
__pyx_t_173 = __pyx_v_ii;
__pyx_t_17 = ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_170 * __pyx_v_Ds.strides[0]) )) + __pyx_t_171)) ))) + (__pyx_v_k * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_172 * __pyx_v_us.strides[0]) )) + __pyx_t_173)) )))));
__pyx_t_174 = 1;
__pyx_t_175 = __pyx_v_ii;
__pyx_t_176 = 1;
__pyx_t_177 = __pyx_v_ii;
__pyx_t_51 = ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_174 * __pyx_v_Ds.strides[0]) )) + __pyx_t_175)) ))) + (__pyx_v_k * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_176 * __pyx_v_us.strides[0]) )) + __pyx_t_177)) )))));
__pyx_v_sol0 = __pyx_t_17;
__pyx_v_sol1 = __pyx_t_51;
+132: if Forbidbis:
__pyx_t_13 = (__pyx_v_Forbidbis != 0);
if (__pyx_t_13) {
/* … */
}
+133: sca0 = (sol0-S1X)*Ds[0,ii] + (sol1-S1Y)*Ds[1,ii]
__pyx_t_178 = 0;
__pyx_t_179 = __pyx_v_ii;
__pyx_t_180 = 1;
__pyx_t_181 = __pyx_v_ii;
__pyx_v_sca0 = (((__pyx_v_sol0 - __pyx_v_S1X) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_178 * __pyx_v_Ds.strides[0]) )) + __pyx_t_179)) )))) + ((__pyx_v_sol1 - __pyx_v_S1Y) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_180 * __pyx_v_Ds.strides[0]) )) + __pyx_t_181)) )))));
+134: sca1 = (sol0-S1X)*S1X + (sol1-S1Y)*S1Y
__pyx_v_sca1 = (((__pyx_v_sol0 - __pyx_v_S1X) * __pyx_v_S1X) + ((__pyx_v_sol1 - __pyx_v_S1Y) * __pyx_v_S1Y));
+135: sca2 = (sol0-S2X)*S2X + (sol1-S2Y)*S2Y
__pyx_v_sca2 = (((__pyx_v_sol0 - __pyx_v_S2X) * __pyx_v_S2X) + ((__pyx_v_sol1 - __pyx_v_S2Y) * __pyx_v_S2Y));
136: #print 1, k, kout, sca0, sca1, sca2
+137: if sca0<0 and sca1<0 and sca2<0:
__pyx_t_12 = ((__pyx_v_sca0 < 0.0) != 0);
if (__pyx_t_12) {
} else {
__pyx_t_13 = __pyx_t_12;
goto __pyx_L58_bool_binop_done;
}
__pyx_t_12 = ((__pyx_v_sca1 < 0.0) != 0);
if (__pyx_t_12) {
} else {
__pyx_t_13 = __pyx_t_12;
goto __pyx_L58_bool_binop_done;
}
__pyx_t_12 = ((__pyx_v_sca2 < 0.0) != 0);
__pyx_t_13 = __pyx_t_12;
__pyx_L58_bool_binop_done:;
if (__pyx_t_13) {
/* … */
}
+138: continue
goto __pyx_L47_continue;
139: # Get the normalized perpendicular vector at intersection
+140: phi = Catan2(sol1,sol0)
__pyx_v_phi = atan2(__pyx_v_sol1, __pyx_v_sol0);
141: # Get the scalar product to determine entry or exit point
+142: sca = Ccos(phi)*vIn[0,jj]*us[0,ii] + Csin(phi)*vIn[0,jj]*us[1,ii] + vIn[1,jj]*us[2,ii]
__pyx_t_182 = 0;
__pyx_t_183 = __pyx_v_jj;
__pyx_t_184 = 0;
__pyx_t_185 = __pyx_v_ii;
__pyx_t_186 = 0;
__pyx_t_187 = __pyx_v_jj;
__pyx_t_188 = 1;
__pyx_t_189 = __pyx_v_ii;
__pyx_t_190 = 1;
__pyx_t_191 = __pyx_v_jj;
__pyx_t_192 = 2;
__pyx_t_193 = __pyx_v_ii;
__pyx_v_sca = ((((cos(__pyx_v_phi) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_vIn.data + __pyx_t_182 * __pyx_v_vIn.strides[0]) )) + __pyx_t_183)) )))) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_184 * __pyx_v_us.strides[0]) )) + __pyx_t_185)) )))) + ((sin(__pyx_v_phi) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_vIn.data + __pyx_t_186 * __pyx_v_vIn.strides[0]) )) + __pyx_t_187)) )))) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_188 * __pyx_v_us.strides[0]) )) + __pyx_t_189)) ))))) + ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_vIn.data + __pyx_t_190 * __pyx_v_vIn.strides[0]) )) + __pyx_t_191)) ))) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_192 * __pyx_v_us.strides[0]) )) + __pyx_t_193)) )))));
+143: if sca<=0 and k<kout:
__pyx_t_12 = ((__pyx_v_sca <= 0.0) != 0);
if (__pyx_t_12) {
} else {
__pyx_t_13 = __pyx_t_12;
goto __pyx_L62_bool_binop_done;
}
__pyx_t_12 = ((__pyx_v_k < __pyx_v_kout) != 0);
__pyx_t_13 = __pyx_t_12;
__pyx_L62_bool_binop_done:;
if (__pyx_t_13) {
/* … */
goto __pyx_L61;
}
+144: kout = k
__pyx_v_kout = __pyx_v_k;
+145: indout = jj
__pyx_v_indout = __pyx_v_jj;
+146: Done = 1
__pyx_v_Done = 1;
+147: print 5, k
__pyx_t_2 = PyFloat_FromDouble(__pyx_v_k); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 147, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__pyx_t_3 = PyTuple_New(2); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 147, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__Pyx_INCREF(__pyx_int_5);
__Pyx_GIVEREF(__pyx_int_5);
PyTuple_SET_ITEM(__pyx_t_3, 0, __pyx_int_5);
__Pyx_GIVEREF(__pyx_t_2);
PyTuple_SET_ITEM(__pyx_t_3, 1, __pyx_t_2);
__pyx_t_2 = 0;
if (__Pyx_Print(0, __pyx_t_3, 1) < 0) __PYX_ERR(0, 147, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
+148: elif sca>=0 and k<min(kin,kout):
__pyx_t_12 = ((__pyx_v_sca >= 0.0) != 0);
if (__pyx_t_12) {
} else {
__pyx_t_13 = __pyx_t_12;
goto __pyx_L64_bool_binop_done;
}
__pyx_t_51 = __pyx_v_kout;
__pyx_t_17 = __pyx_v_kin;
if (((__pyx_t_51 < __pyx_t_17) != 0)) {
__pyx_t_105 = __pyx_t_51;
} else {
__pyx_t_105 = __pyx_t_17;
}
__pyx_t_12 = ((__pyx_v_k < __pyx_t_105) != 0);
__pyx_t_13 = __pyx_t_12;
__pyx_L64_bool_binop_done:;
if (__pyx_t_13) {
/* … */
}
__pyx_L61:;
+149: kin = k
__pyx_v_kin = __pyx_v_k;
+150: print 6, k
__pyx_t_3 = PyFloat_FromDouble(__pyx_v_k); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 150, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__pyx_t_2 = PyTuple_New(2); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 150, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__Pyx_INCREF(__pyx_int_6);
__Pyx_GIVEREF(__pyx_int_6);
PyTuple_SET_ITEM(__pyx_t_2, 0, __pyx_int_6);
__Pyx_GIVEREF(__pyx_t_3);
PyTuple_SET_ITEM(__pyx_t_2, 1, __pyx_t_3);
__pyx_t_3 = 0;
if (__Pyx_Print(0, __pyx_t_2, 1) < 0) __PYX_ERR(0, 150, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
151:
+152: elif A**2>=EpsA**2 and B**2>A*C:
__pyx_t_12 = ((pow(__pyx_v_A, 2.0) >= pow(__pyx_v_EpsA, 2.0)) != 0);
if (__pyx_t_12) {
} else {
__pyx_t_13 = __pyx_t_12;
goto __pyx_L66_bool_binop_done;
}
__pyx_t_12 = ((pow(__pyx_v_B, 2.0) > (__pyx_v_A * __pyx_v_C)) != 0);
__pyx_t_13 = __pyx_t_12;
__pyx_L66_bool_binop_done:;
if (__pyx_t_13) {
/* … */
}
__pyx_L49:;
__pyx_L47_continue:;
}
}
__pyx_L11:;
+153: sqd = Csqrt(B**2-A*C)
__pyx_v_sqd = sqrt((pow(__pyx_v_B, 2.0) - (__pyx_v_A * __pyx_v_C)));
154: # First solution
+155: q = (-B + sqd)/A
__pyx_v_q = (((-__pyx_v_B) + __pyx_v_sqd) / __pyx_v_A);
+156: if q>=0. and q<1.:
__pyx_t_12 = ((__pyx_v_q >= 0.) != 0);
if (__pyx_t_12) {
} else {
__pyx_t_13 = __pyx_t_12;
goto __pyx_L69_bool_binop_done;
}
__pyx_t_12 = ((__pyx_v_q < 1.) != 0);
__pyx_t_13 = __pyx_t_12;
__pyx_L69_bool_binop_done:;
if (__pyx_t_13) {
/* … */
}
+157: k = (q*v1 - (Ds[2,ii]-VPoly[1,jj]))/us[2,ii]
__pyx_t_194 = 2;
__pyx_t_195 = __pyx_v_ii;
__pyx_t_196 = 1;
__pyx_t_197 = __pyx_v_jj;
__pyx_t_198 = 2;
__pyx_t_199 = __pyx_v_ii;
__pyx_v_k = (((__pyx_v_q * __pyx_v_v1) - ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_194 * __pyx_v_Ds.strides[0]) )) + __pyx_t_195)) ))) - (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_VPoly.data + __pyx_t_196 * __pyx_v_VPoly.strides[0]) )) + __pyx_t_197)) ))))) / (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_198 * __pyx_v_us.strides[0]) )) + __pyx_t_199)) ))));
+158: if k>=0.:
__pyx_t_13 = ((__pyx_v_k >= 0.) != 0);
if (__pyx_t_13) {
/* … */
}
+159: sol0, sol1 = Ds[0,ii] + k*us[0,ii], Ds[1,ii] + k*us[1,ii]
__pyx_t_200 = 0;
__pyx_t_201 = __pyx_v_ii;
__pyx_t_202 = 0;
__pyx_t_203 = __pyx_v_ii;
__pyx_t_105 = ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_200 * __pyx_v_Ds.strides[0]) )) + __pyx_t_201)) ))) + (__pyx_v_k * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_202 * __pyx_v_us.strides[0]) )) + __pyx_t_203)) )))));
__pyx_t_204 = 1;
__pyx_t_205 = __pyx_v_ii;
__pyx_t_206 = 1;
__pyx_t_207 = __pyx_v_ii;
__pyx_t_51 = ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_204 * __pyx_v_Ds.strides[0]) )) + __pyx_t_205)) ))) + (__pyx_v_k * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_206 * __pyx_v_us.strides[0]) )) + __pyx_t_207)) )))));
__pyx_v_sol0 = __pyx_t_105;
__pyx_v_sol1 = __pyx_t_51;
+160: if Forbidbis:
__pyx_t_13 = (__pyx_v_Forbidbis != 0);
if (__pyx_t_13) {
/* … */
}
+161: sca0 = (sol0-S1X)*Ds[0,ii] + (sol1-S1Y)*Ds[1,ii]
__pyx_t_208 = 0;
__pyx_t_209 = __pyx_v_ii;
__pyx_t_210 = 1;
__pyx_t_211 = __pyx_v_ii;
__pyx_v_sca0 = (((__pyx_v_sol0 - __pyx_v_S1X) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_208 * __pyx_v_Ds.strides[0]) )) + __pyx_t_209)) )))) + ((__pyx_v_sol1 - __pyx_v_S1Y) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_210 * __pyx_v_Ds.strides[0]) )) + __pyx_t_211)) )))));
+162: sca1 = (sol0-S1X)*S1X + (sol1-S1Y)*S1Y
__pyx_v_sca1 = (((__pyx_v_sol0 - __pyx_v_S1X) * __pyx_v_S1X) + ((__pyx_v_sol1 - __pyx_v_S1Y) * __pyx_v_S1Y));
+163: sca2 = (sol0-S2X)*S2X + (sol1-S2Y)*S2Y
__pyx_v_sca2 = (((__pyx_v_sol0 - __pyx_v_S2X) * __pyx_v_S2X) + ((__pyx_v_sol1 - __pyx_v_S2Y) * __pyx_v_S2Y));
164: #print 2, k, kout, sca0, sca1, sca2
+165: if not Forbidbis or (Forbidbis and not (sca0<0 and sca1<0 and sca2<0)):
__pyx_t_12 = ((!(__pyx_v_Forbidbis != 0)) != 0);
if (!__pyx_t_12) {
} else {
__pyx_t_13 = __pyx_t_12;
goto __pyx_L74_bool_binop_done;
}
__pyx_t_12 = (__pyx_v_Forbidbis != 0);
if (__pyx_t_12) {
} else {
__pyx_t_13 = __pyx_t_12;
goto __pyx_L74_bool_binop_done;
}
__pyx_t_92 = ((__pyx_v_sca0 < 0.0) != 0);
if (__pyx_t_92) {
} else {
__pyx_t_12 = __pyx_t_92;
goto __pyx_L77_bool_binop_done;
}
__pyx_t_92 = ((__pyx_v_sca1 < 0.0) != 0);
if (__pyx_t_92) {
} else {
__pyx_t_12 = __pyx_t_92;
goto __pyx_L77_bool_binop_done;
}
__pyx_t_92 = ((__pyx_v_sca2 < 0.0) != 0);
__pyx_t_12 = __pyx_t_92;
__pyx_L77_bool_binop_done:;
__pyx_t_92 = ((!__pyx_t_12) != 0);
__pyx_t_13 = __pyx_t_92;
__pyx_L74_bool_binop_done:;
if (__pyx_t_13) {
/* … */
}
166: # Get the normalized perpendicular vector at intersection
+167: phi = Catan2(sol1,sol0)
__pyx_v_phi = atan2(__pyx_v_sol1, __pyx_v_sol0);
168: # Get the scalar product to determine entry or exit point
+169: sca = Ccos(phi)*vIn[0,jj]*us[0,ii] + Csin(phi)*vIn[0,jj]*us[1,ii] + vIn[1,jj]*us[2,ii]
__pyx_t_212 = 0;
__pyx_t_213 = __pyx_v_jj;
__pyx_t_214 = 0;
__pyx_t_215 = __pyx_v_ii;
__pyx_t_216 = 0;
__pyx_t_217 = __pyx_v_jj;
__pyx_t_218 = 1;
__pyx_t_219 = __pyx_v_ii;
__pyx_t_220 = 1;
__pyx_t_221 = __pyx_v_jj;
__pyx_t_222 = 2;
__pyx_t_223 = __pyx_v_ii;
__pyx_v_sca = ((((cos(__pyx_v_phi) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_vIn.data + __pyx_t_212 * __pyx_v_vIn.strides[0]) )) + __pyx_t_213)) )))) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_214 * __pyx_v_us.strides[0]) )) + __pyx_t_215)) )))) + ((sin(__pyx_v_phi) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_vIn.data + __pyx_t_216 * __pyx_v_vIn.strides[0]) )) + __pyx_t_217)) )))) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_218 * __pyx_v_us.strides[0]) )) + __pyx_t_219)) ))))) + ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_vIn.data + __pyx_t_220 * __pyx_v_vIn.strides[0]) )) + __pyx_t_221)) ))) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_222 * __pyx_v_us.strides[0]) )) + __pyx_t_223)) )))));
+170: if sca<=0 and k<kout:
__pyx_t_92 = ((__pyx_v_sca <= 0.0) != 0);
if (__pyx_t_92) {
} else {
__pyx_t_13 = __pyx_t_92;
goto __pyx_L81_bool_binop_done;
}
__pyx_t_92 = ((__pyx_v_k < __pyx_v_kout) != 0);
__pyx_t_13 = __pyx_t_92;
__pyx_L81_bool_binop_done:;
if (__pyx_t_13) {
/* … */
goto __pyx_L80;
}
+171: kout = k
__pyx_v_kout = __pyx_v_k;
+172: indout = jj
__pyx_v_indout = __pyx_v_jj;
+173: Done = 1
__pyx_v_Done = 1;
+174: print 7, k, q, A, B, C, sqd
__pyx_t_2 = PyFloat_FromDouble(__pyx_v_k); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 174, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__pyx_t_3 = PyFloat_FromDouble(__pyx_v_q); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 174, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__pyx_t_6 = PyFloat_FromDouble(__pyx_v_A); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 174, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_6);
__pyx_t_4 = PyFloat_FromDouble(__pyx_v_B); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 174, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__pyx_t_16 = PyFloat_FromDouble(__pyx_v_C); if (unlikely(!__pyx_t_16)) __PYX_ERR(0, 174, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_16);
__pyx_t_5 = PyFloat_FromDouble(__pyx_v_sqd); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 174, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__pyx_t_1 = PyTuple_New(7); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 174, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__Pyx_INCREF(__pyx_int_7);
__Pyx_GIVEREF(__pyx_int_7);
PyTuple_SET_ITEM(__pyx_t_1, 0, __pyx_int_7);
__Pyx_GIVEREF(__pyx_t_2);
PyTuple_SET_ITEM(__pyx_t_1, 1, __pyx_t_2);
__Pyx_GIVEREF(__pyx_t_3);
PyTuple_SET_ITEM(__pyx_t_1, 2, __pyx_t_3);
__Pyx_GIVEREF(__pyx_t_6);
PyTuple_SET_ITEM(__pyx_t_1, 3, __pyx_t_6);
__Pyx_GIVEREF(__pyx_t_4);
PyTuple_SET_ITEM(__pyx_t_1, 4, __pyx_t_4);
__Pyx_GIVEREF(__pyx_t_16);
PyTuple_SET_ITEM(__pyx_t_1, 5, __pyx_t_16);
__Pyx_GIVEREF(__pyx_t_5);
PyTuple_SET_ITEM(__pyx_t_1, 6, __pyx_t_5);
__pyx_t_2 = 0;
__pyx_t_3 = 0;
__pyx_t_6 = 0;
__pyx_t_4 = 0;
__pyx_t_16 = 0;
__pyx_t_5 = 0;
if (__Pyx_Print(0, __pyx_t_1, 1) < 0) __PYX_ERR(0, 174, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+175: elif sca>=0 and k<min(kin,kout):
__pyx_t_92 = ((__pyx_v_sca >= 0.0) != 0);
if (__pyx_t_92) {
} else {
__pyx_t_13 = __pyx_t_92;
goto __pyx_L83_bool_binop_done;
}
__pyx_t_51 = __pyx_v_kout;
__pyx_t_105 = __pyx_v_kin;
if (((__pyx_t_51 < __pyx_t_105) != 0)) {
__pyx_t_17 = __pyx_t_51;
} else {
__pyx_t_17 = __pyx_t_105;
}
__pyx_t_92 = ((__pyx_v_k < __pyx_t_17) != 0);
__pyx_t_13 = __pyx_t_92;
__pyx_L83_bool_binop_done:;
if (__pyx_t_13) {
/* … */
}
__pyx_L80:;
+176: kin = k
__pyx_v_kin = __pyx_v_k;
+177: print 8, k, jj
__pyx_t_1 = PyFloat_FromDouble(__pyx_v_k); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 177, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__pyx_t_5 = __Pyx_PyInt_From_int(__pyx_v_jj); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 177, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__pyx_t_16 = PyTuple_New(3); if (unlikely(!__pyx_t_16)) __PYX_ERR(0, 177, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_16);
__Pyx_INCREF(__pyx_int_8);
__Pyx_GIVEREF(__pyx_int_8);
PyTuple_SET_ITEM(__pyx_t_16, 0, __pyx_int_8);
__Pyx_GIVEREF(__pyx_t_1);
PyTuple_SET_ITEM(__pyx_t_16, 1, __pyx_t_1);
__Pyx_GIVEREF(__pyx_t_5);
PyTuple_SET_ITEM(__pyx_t_16, 2, __pyx_t_5);
__pyx_t_1 = 0;
__pyx_t_5 = 0;
if (__Pyx_Print(0, __pyx_t_16, 1) < 0) __PYX_ERR(0, 177, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_16); __pyx_t_16 = 0;
178:
179: # Second solution
+180: q = (-B - sqd)/A
__pyx_v_q = (((-__pyx_v_B) - __pyx_v_sqd) / __pyx_v_A);
+181: if q>=0. and q<1.:
__pyx_t_92 = ((__pyx_v_q >= 0.) != 0);
if (__pyx_t_92) {
} else {
__pyx_t_13 = __pyx_t_92;
goto __pyx_L86_bool_binop_done;
}
__pyx_t_92 = ((__pyx_v_q < 1.) != 0);
__pyx_t_13 = __pyx_t_92;
__pyx_L86_bool_binop_done:;
if (__pyx_t_13) {
/* … */
}
+182: k = (q*v1 - (Ds[2,ii]-VPoly[1,jj]))/us[2,ii]
__pyx_t_224 = 2;
__pyx_t_225 = __pyx_v_ii;
__pyx_t_226 = 1;
__pyx_t_227 = __pyx_v_jj;
__pyx_t_228 = 2;
__pyx_t_229 = __pyx_v_ii;
__pyx_v_k = (((__pyx_v_q * __pyx_v_v1) - ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_224 * __pyx_v_Ds.strides[0]) )) + __pyx_t_225)) ))) - (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_VPoly.data + __pyx_t_226 * __pyx_v_VPoly.strides[0]) )) + __pyx_t_227)) ))))) / (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_228 * __pyx_v_us.strides[0]) )) + __pyx_t_229)) ))));
183:
+184: if k>=0.:
__pyx_t_13 = ((__pyx_v_k >= 0.) != 0);
if (__pyx_t_13) {
/* … */
}
+185: sol0, sol1 = Ds[0,ii] + k*us[0,ii], Ds[1,ii] + k*us[1,ii]
__pyx_t_230 = 0;
__pyx_t_231 = __pyx_v_ii;
__pyx_t_232 = 0;
__pyx_t_233 = __pyx_v_ii;
__pyx_t_17 = ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_230 * __pyx_v_Ds.strides[0]) )) + __pyx_t_231)) ))) + (__pyx_v_k * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_232 * __pyx_v_us.strides[0]) )) + __pyx_t_233)) )))));
__pyx_t_234 = 1;
__pyx_t_235 = __pyx_v_ii;
__pyx_t_236 = 1;
__pyx_t_237 = __pyx_v_ii;
__pyx_t_51 = ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_234 * __pyx_v_Ds.strides[0]) )) + __pyx_t_235)) ))) + (__pyx_v_k * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_236 * __pyx_v_us.strides[0]) )) + __pyx_t_237)) )))));
__pyx_v_sol0 = __pyx_t_17;
__pyx_v_sol1 = __pyx_t_51;
+186: if Forbidbis:
__pyx_t_13 = (__pyx_v_Forbidbis != 0);
if (__pyx_t_13) {
/* … */
}
+187: sca0 = (sol0-S1X)*Ds[0,ii] + (sol1-S1Y)*Ds[1,ii]
__pyx_t_238 = 0;
__pyx_t_239 = __pyx_v_ii;
__pyx_t_240 = 1;
__pyx_t_241 = __pyx_v_ii;
__pyx_v_sca0 = (((__pyx_v_sol0 - __pyx_v_S1X) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_238 * __pyx_v_Ds.strides[0]) )) + __pyx_t_239)) )))) + ((__pyx_v_sol1 - __pyx_v_S1Y) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_240 * __pyx_v_Ds.strides[0]) )) + __pyx_t_241)) )))));
+188: sca1 = (sol0-S1X)*S1X + (sol1-S1Y)*S1Y
__pyx_v_sca1 = (((__pyx_v_sol0 - __pyx_v_S1X) * __pyx_v_S1X) + ((__pyx_v_sol1 - __pyx_v_S1Y) * __pyx_v_S1Y));
+189: sca2 = (sol0-S2X)*S2X + (sol1-S2Y)*S2Y
__pyx_v_sca2 = (((__pyx_v_sol0 - __pyx_v_S2X) * __pyx_v_S2X) + ((__pyx_v_sol1 - __pyx_v_S2Y) * __pyx_v_S2Y));
190: #print 3, k, kout, sca0, sca1, sca2
+191: if not Forbidbis or (Forbidbis and not (sca0<0 and sca1<0 and sca2<0)):
__pyx_t_92 = ((!(__pyx_v_Forbidbis != 0)) != 0);
if (!__pyx_t_92) {
} else {
__pyx_t_13 = __pyx_t_92;
goto __pyx_L91_bool_binop_done;
}
__pyx_t_92 = (__pyx_v_Forbidbis != 0);
if (__pyx_t_92) {
} else {
__pyx_t_13 = __pyx_t_92;
goto __pyx_L91_bool_binop_done;
}
__pyx_t_12 = ((__pyx_v_sca0 < 0.0) != 0);
if (__pyx_t_12) {
} else {
__pyx_t_92 = __pyx_t_12;
goto __pyx_L94_bool_binop_done;
}
__pyx_t_12 = ((__pyx_v_sca1 < 0.0) != 0);
if (__pyx_t_12) {
} else {
__pyx_t_92 = __pyx_t_12;
goto __pyx_L94_bool_binop_done;
}
__pyx_t_12 = ((__pyx_v_sca2 < 0.0) != 0);
__pyx_t_92 = __pyx_t_12;
__pyx_L94_bool_binop_done:;
__pyx_t_12 = ((!__pyx_t_92) != 0);
__pyx_t_13 = __pyx_t_12;
__pyx_L91_bool_binop_done:;
if (__pyx_t_13) {
/* … */
}
192: # Get the normalized perpendicular vector at intersection
+193: phi = Catan2(sol1,sol0)
__pyx_v_phi = atan2(__pyx_v_sol1, __pyx_v_sol0);
194: # Get the scalar product to determine entry or exit point
+195: sca = Ccos(phi)*vIn[0,jj]*us[0,ii] + Csin(phi)*vIn[0,jj]*us[1,ii] + vIn[1,jj]*us[2,ii]
__pyx_t_242 = 0;
__pyx_t_243 = __pyx_v_jj;
__pyx_t_244 = 0;
__pyx_t_245 = __pyx_v_ii;
__pyx_t_246 = 0;
__pyx_t_247 = __pyx_v_jj;
__pyx_t_248 = 1;
__pyx_t_249 = __pyx_v_ii;
__pyx_t_250 = 1;
__pyx_t_251 = __pyx_v_jj;
__pyx_t_252 = 2;
__pyx_t_253 = __pyx_v_ii;
__pyx_v_sca = ((((cos(__pyx_v_phi) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_vIn.data + __pyx_t_242 * __pyx_v_vIn.strides[0]) )) + __pyx_t_243)) )))) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_244 * __pyx_v_us.strides[0]) )) + __pyx_t_245)) )))) + ((sin(__pyx_v_phi) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_vIn.data + __pyx_t_246 * __pyx_v_vIn.strides[0]) )) + __pyx_t_247)) )))) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_248 * __pyx_v_us.strides[0]) )) + __pyx_t_249)) ))))) + ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_vIn.data + __pyx_t_250 * __pyx_v_vIn.strides[0]) )) + __pyx_t_251)) ))) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_252 * __pyx_v_us.strides[0]) )) + __pyx_t_253)) )))));
+196: if sca<=0 and k<kout:
__pyx_t_12 = ((__pyx_v_sca <= 0.0) != 0);
if (__pyx_t_12) {
} else {
__pyx_t_13 = __pyx_t_12;
goto __pyx_L98_bool_binop_done;
}
__pyx_t_12 = ((__pyx_v_k < __pyx_v_kout) != 0);
__pyx_t_13 = __pyx_t_12;
__pyx_L98_bool_binop_done:;
if (__pyx_t_13) {
/* … */
goto __pyx_L97;
}
+197: kout = k
__pyx_v_kout = __pyx_v_k;
+198: indout = jj
__pyx_v_indout = __pyx_v_jj;
+199: Done = 1
__pyx_v_Done = 1;
+200: print 9, k, jj
__pyx_t_16 = PyFloat_FromDouble(__pyx_v_k); if (unlikely(!__pyx_t_16)) __PYX_ERR(0, 200, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_16);
__pyx_t_5 = __Pyx_PyInt_From_int(__pyx_v_jj); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 200, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__pyx_t_1 = PyTuple_New(3); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 200, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__Pyx_INCREF(__pyx_int_9);
__Pyx_GIVEREF(__pyx_int_9);
PyTuple_SET_ITEM(__pyx_t_1, 0, __pyx_int_9);
__Pyx_GIVEREF(__pyx_t_16);
PyTuple_SET_ITEM(__pyx_t_1, 1, __pyx_t_16);
__Pyx_GIVEREF(__pyx_t_5);
PyTuple_SET_ITEM(__pyx_t_1, 2, __pyx_t_5);
__pyx_t_16 = 0;
__pyx_t_5 = 0;
if (__Pyx_Print(0, __pyx_t_1, 1) < 0) __PYX_ERR(0, 200, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+201: elif sca>=0 and k<min(kin,kout):
__pyx_t_12 = ((__pyx_v_sca >= 0.0) != 0);
if (__pyx_t_12) {
} else {
__pyx_t_13 = __pyx_t_12;
goto __pyx_L100_bool_binop_done;
}
__pyx_t_51 = __pyx_v_kout;
__pyx_t_17 = __pyx_v_kin;
if (((__pyx_t_51 < __pyx_t_17) != 0)) {
__pyx_t_105 = __pyx_t_51;
} else {
__pyx_t_105 = __pyx_t_17;
}
__pyx_t_12 = ((__pyx_v_k < __pyx_t_105) != 0);
__pyx_t_13 = __pyx_t_12;
__pyx_L100_bool_binop_done:;
if (__pyx_t_13) {
/* … */
}
__pyx_L97:;
+202: kin = k
__pyx_v_kin = __pyx_v_k;
+203: print 10, k, q, A, B, C, sqd, v0, v1, jj
__pyx_t_1 = PyFloat_FromDouble(__pyx_v_k); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 203, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__pyx_t_5 = PyFloat_FromDouble(__pyx_v_q); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 203, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__pyx_t_16 = PyFloat_FromDouble(__pyx_v_A); if (unlikely(!__pyx_t_16)) __PYX_ERR(0, 203, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_16);
__pyx_t_4 = PyFloat_FromDouble(__pyx_v_B); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 203, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__pyx_t_6 = PyFloat_FromDouble(__pyx_v_C); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 203, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_6);
__pyx_t_3 = PyFloat_FromDouble(__pyx_v_sqd); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 203, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__pyx_t_2 = PyFloat_FromDouble(__pyx_v_v0); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 203, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__pyx_t_14 = PyFloat_FromDouble(__pyx_v_v1); if (unlikely(!__pyx_t_14)) __PYX_ERR(0, 203, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_14);
__pyx_t_254 = __Pyx_PyInt_From_int(__pyx_v_jj); if (unlikely(!__pyx_t_254)) __PYX_ERR(0, 203, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_254);
__pyx_t_255 = PyTuple_New(10); if (unlikely(!__pyx_t_255)) __PYX_ERR(0, 203, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_255);
__Pyx_INCREF(__pyx_int_10);
__Pyx_GIVEREF(__pyx_int_10);
PyTuple_SET_ITEM(__pyx_t_255, 0, __pyx_int_10);
__Pyx_GIVEREF(__pyx_t_1);
PyTuple_SET_ITEM(__pyx_t_255, 1, __pyx_t_1);
__Pyx_GIVEREF(__pyx_t_5);
PyTuple_SET_ITEM(__pyx_t_255, 2, __pyx_t_5);
__Pyx_GIVEREF(__pyx_t_16);
PyTuple_SET_ITEM(__pyx_t_255, 3, __pyx_t_16);
__Pyx_GIVEREF(__pyx_t_4);
PyTuple_SET_ITEM(__pyx_t_255, 4, __pyx_t_4);
__Pyx_GIVEREF(__pyx_t_6);
PyTuple_SET_ITEM(__pyx_t_255, 5, __pyx_t_6);
__Pyx_GIVEREF(__pyx_t_3);
PyTuple_SET_ITEM(__pyx_t_255, 6, __pyx_t_3);
__Pyx_GIVEREF(__pyx_t_2);
PyTuple_SET_ITEM(__pyx_t_255, 7, __pyx_t_2);
__Pyx_GIVEREF(__pyx_t_14);
PyTuple_SET_ITEM(__pyx_t_255, 8, __pyx_t_14);
__Pyx_GIVEREF(__pyx_t_254);
PyTuple_SET_ITEM(__pyx_t_255, 9, __pyx_t_254);
__pyx_t_1 = 0;
__pyx_t_5 = 0;
__pyx_t_16 = 0;
__pyx_t_4 = 0;
__pyx_t_6 = 0;
__pyx_t_3 = 0;
__pyx_t_2 = 0;
__pyx_t_14 = 0;
__pyx_t_254 = 0;
if (__Pyx_Print(0, __pyx_t_255, 1) < 0) __PYX_ERR(0, 203, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_255); __pyx_t_255 = 0;
204:
+205: if Done==1:
__pyx_t_13 = ((__pyx_v_Done == 1) != 0);
if (__pyx_t_13) {
/* … */
}
+206: indOut[ii] = indout
__pyx_t_256 = __pyx_v_ii;
*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_indOut.data) + __pyx_t_256)) )) = __pyx_v_indout;
+207: SOut[0,ii] = Ds[0,ii] + kout*us[0,ii]
__pyx_t_257 = 0;
__pyx_t_258 = __pyx_v_ii;
__pyx_t_259 = 0;
__pyx_t_260 = __pyx_v_ii;
__pyx_t_261 = 0;
__pyx_t_262 = __pyx_v_ii;
*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_SOut.data + __pyx_t_261 * __pyx_v_SOut.strides[0]) )) + __pyx_t_262)) )) = ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_257 * __pyx_v_Ds.strides[0]) )) + __pyx_t_258)) ))) + (__pyx_v_kout * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_259 * __pyx_v_us.strides[0]) )) + __pyx_t_260)) )))));
+208: SOut[1,ii] = Ds[1,ii] + kout*us[1,ii]
__pyx_t_263 = 1;
__pyx_t_264 = __pyx_v_ii;
__pyx_t_265 = 1;
__pyx_t_266 = __pyx_v_ii;
__pyx_t_267 = 1;
__pyx_t_268 = __pyx_v_ii;
*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_SOut.data + __pyx_t_267 * __pyx_v_SOut.strides[0]) )) + __pyx_t_268)) )) = ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_263 * __pyx_v_Ds.strides[0]) )) + __pyx_t_264)) ))) + (__pyx_v_kout * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_265 * __pyx_v_us.strides[0]) )) + __pyx_t_266)) )))));
+209: SOut[2,ii] = Ds[2,ii] + kout*us[2,ii]
__pyx_t_269 = 2;
__pyx_t_270 = __pyx_v_ii;
__pyx_t_271 = 2;
__pyx_t_272 = __pyx_v_ii;
__pyx_t_273 = 2;
__pyx_t_274 = __pyx_v_ii;
*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_SOut.data + __pyx_t_273 * __pyx_v_SOut.strides[0]) )) + __pyx_t_274)) )) = ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_269 * __pyx_v_Ds.strides[0]) )) + __pyx_t_270)) ))) + (__pyx_v_kout * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_271 * __pyx_v_us.strides[0]) )) + __pyx_t_272)) )))));
+210: phi = Catan2(SOut[1,ii],SOut[0,ii])
__pyx_t_275 = 1;
__pyx_t_276 = __pyx_v_ii;
__pyx_t_277 = 0;
__pyx_t_278 = __pyx_v_ii;
__pyx_v_phi = atan2((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_SOut.data + __pyx_t_275 * __pyx_v_SOut.strides[0]) )) + __pyx_t_276)) ))), (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_SOut.data + __pyx_t_277 * __pyx_v_SOut.strides[0]) )) + __pyx_t_278)) ))));
+211: VPerp[0,ii] = Ccos(phi)*vIn[0,indout]
__pyx_t_279 = 0;
__pyx_t_280 = __pyx_v_indout;
__pyx_t_281 = 0;
__pyx_t_282 = __pyx_v_ii;
*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_VPerp.data + __pyx_t_281 * __pyx_v_VPerp.strides[0]) )) + __pyx_t_282)) )) = (cos(__pyx_v_phi) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_vIn.data + __pyx_t_279 * __pyx_v_vIn.strides[0]) )) + __pyx_t_280)) ))));
+212: VPerp[1,ii] = Csin(phi)*vIn[0,indout]
__pyx_t_283 = 0;
__pyx_t_284 = __pyx_v_indout;
__pyx_t_285 = 1;
__pyx_t_286 = __pyx_v_ii;
*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_VPerp.data + __pyx_t_285 * __pyx_v_VPerp.strides[0]) )) + __pyx_t_286)) )) = (sin(__pyx_v_phi) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_vIn.data + __pyx_t_283 * __pyx_v_vIn.strides[0]) )) + __pyx_t_284)) ))));
+213: VPerp[2,ii] = vIn[1,indout]
__pyx_t_287 = 1;
__pyx_t_288 = __pyx_v_indout;
__pyx_t_289 = 2;
__pyx_t_290 = __pyx_v_ii;
*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_VPerp.data + __pyx_t_289 * __pyx_v_VPerp.strides[0]) )) + __pyx_t_290)) )) = (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_vIn.data + __pyx_t_287 * __pyx_v_vIn.strides[0]) )) + __pyx_t_288)) )));
+214: if kin<kout:
__pyx_t_13 = ((__pyx_v_kin < __pyx_v_kout) != 0);
if (__pyx_t_13) {
/* … */
}
}
+215: SIn[0,ii] = Ds[0,ii] + kin*us[0,ii]
__pyx_t_291 = 0;
__pyx_t_292 = __pyx_v_ii;
__pyx_t_293 = 0;
__pyx_t_294 = __pyx_v_ii;
__pyx_t_295 = 0;
__pyx_t_296 = __pyx_v_ii;
*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_SIn.data + __pyx_t_295 * __pyx_v_SIn.strides[0]) )) + __pyx_t_296)) )) = ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_291 * __pyx_v_Ds.strides[0]) )) + __pyx_t_292)) ))) + (__pyx_v_kin * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_293 * __pyx_v_us.strides[0]) )) + __pyx_t_294)) )))));
+216: SIn[1,ii] = Ds[1,ii] + kin*us[1,ii]
__pyx_t_297 = 1;
__pyx_t_298 = __pyx_v_ii;
__pyx_t_299 = 1;
__pyx_t_300 = __pyx_v_ii;
__pyx_t_301 = 1;
__pyx_t_302 = __pyx_v_ii;
*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_SIn.data + __pyx_t_301 * __pyx_v_SIn.strides[0]) )) + __pyx_t_302)) )) = ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_297 * __pyx_v_Ds.strides[0]) )) + __pyx_t_298)) ))) + (__pyx_v_kin * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_299 * __pyx_v_us.strides[0]) )) + __pyx_t_300)) )))));
+217: SIn[2,ii] = Ds[2,ii] + kin*us[2,ii]
__pyx_t_303 = 2;
__pyx_t_304 = __pyx_v_ii;
__pyx_t_305 = 2;
__pyx_t_306 = __pyx_v_ii;
__pyx_t_307 = 2;
__pyx_t_308 = __pyx_v_ii;
*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_SIn.data + __pyx_t_307 * __pyx_v_SIn.strides[0]) )) + __pyx_t_308)) )) = ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_Ds.data + __pyx_t_303 * __pyx_v_Ds.strides[0]) )) + __pyx_t_304)) ))) + (__pyx_v_kin * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_us.data + __pyx_t_305 * __pyx_v_us.strides[0]) )) + __pyx_t_306)) )))));
+218: return np.asarray(SIn), np.asarray(SOut), np.asarray(VPerp), np.asarray(indOut)
__Pyx_XDECREF(__pyx_r);
__pyx_t_254 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_254)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_254);
__pyx_t_14 = __Pyx_PyObject_GetAttrStr(__pyx_t_254, __pyx_n_s_asarray); if (unlikely(!__pyx_t_14)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_14);
__Pyx_DECREF(__pyx_t_254); __pyx_t_254 = 0;
__pyx_t_254 = __pyx_memoryview_fromslice(__pyx_v_SIn, 2, (PyObject *(*)(char *)) __pyx_memview_get_double, (int (*)(char *, PyObject *)) __pyx_memview_set_double, 0);; if (unlikely(!__pyx_t_254)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_254);
__pyx_t_2 = NULL;
if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_14))) {
__pyx_t_2 = PyMethod_GET_SELF(__pyx_t_14);
if (likely(__pyx_t_2)) {
PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_14);
__Pyx_INCREF(__pyx_t_2);
__Pyx_INCREF(function);
__Pyx_DECREF_SET(__pyx_t_14, function);
}
}
if (!__pyx_t_2) {
__pyx_t_255 = __Pyx_PyObject_CallOneArg(__pyx_t_14, __pyx_t_254); if (unlikely(!__pyx_t_255)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_254); __pyx_t_254 = 0;
__Pyx_GOTREF(__pyx_t_255);
} else {
#if CYTHON_FAST_PYCALL
if (PyFunction_Check(__pyx_t_14)) {
PyObject *__pyx_temp[2] = {__pyx_t_2, __pyx_t_254};
__pyx_t_255 = __Pyx_PyFunction_FastCall(__pyx_t_14, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_255)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
__Pyx_GOTREF(__pyx_t_255);
__Pyx_DECREF(__pyx_t_254); __pyx_t_254 = 0;
} else
#endif
#if CYTHON_FAST_PYCCALL
if (__Pyx_PyFastCFunction_Check(__pyx_t_14)) {
PyObject *__pyx_temp[2] = {__pyx_t_2, __pyx_t_254};
__pyx_t_255 = __Pyx_PyCFunction_FastCall(__pyx_t_14, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_255)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
__Pyx_GOTREF(__pyx_t_255);
__Pyx_DECREF(__pyx_t_254); __pyx_t_254 = 0;
} else
#endif
{
__pyx_t_3 = PyTuple_New(1+1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__Pyx_GIVEREF(__pyx_t_2); PyTuple_SET_ITEM(__pyx_t_3, 0, __pyx_t_2); __pyx_t_2 = NULL;
__Pyx_GIVEREF(__pyx_t_254);
PyTuple_SET_ITEM(__pyx_t_3, 0+1, __pyx_t_254);
__pyx_t_254 = 0;
__pyx_t_255 = __Pyx_PyObject_Call(__pyx_t_14, __pyx_t_3, NULL); if (unlikely(!__pyx_t_255)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_255);
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
}
}
__Pyx_DECREF(__pyx_t_14); __pyx_t_14 = 0;
__pyx_t_3 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__pyx_t_254 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_asarray); if (unlikely(!__pyx_t_254)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_254);
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__pyx_t_3 = __pyx_memoryview_fromslice(__pyx_v_SOut, 2, (PyObject *(*)(char *)) __pyx_memview_get_double, (int (*)(char *, PyObject *)) __pyx_memview_set_double, 0);; if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__pyx_t_2 = NULL;
if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_254))) {
__pyx_t_2 = PyMethod_GET_SELF(__pyx_t_254);
if (likely(__pyx_t_2)) {
PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_254);
__Pyx_INCREF(__pyx_t_2);
__Pyx_INCREF(function);
__Pyx_DECREF_SET(__pyx_t_254, function);
}
}
if (!__pyx_t_2) {
__pyx_t_14 = __Pyx_PyObject_CallOneArg(__pyx_t_254, __pyx_t_3); if (unlikely(!__pyx_t_14)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__Pyx_GOTREF(__pyx_t_14);
} else {
#if CYTHON_FAST_PYCALL
if (PyFunction_Check(__pyx_t_254)) {
PyObject *__pyx_temp[2] = {__pyx_t_2, __pyx_t_3};
__pyx_t_14 = __Pyx_PyFunction_FastCall(__pyx_t_254, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_14)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
__Pyx_GOTREF(__pyx_t_14);
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
} else
#endif
#if CYTHON_FAST_PYCCALL
if (__Pyx_PyFastCFunction_Check(__pyx_t_254)) {
PyObject *__pyx_temp[2] = {__pyx_t_2, __pyx_t_3};
__pyx_t_14 = __Pyx_PyCFunction_FastCall(__pyx_t_254, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_14)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
__Pyx_GOTREF(__pyx_t_14);
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
} else
#endif
{
__pyx_t_6 = PyTuple_New(1+1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_6);
__Pyx_GIVEREF(__pyx_t_2); PyTuple_SET_ITEM(__pyx_t_6, 0, __pyx_t_2); __pyx_t_2 = NULL;
__Pyx_GIVEREF(__pyx_t_3);
PyTuple_SET_ITEM(__pyx_t_6, 0+1, __pyx_t_3);
__pyx_t_3 = 0;
__pyx_t_14 = __Pyx_PyObject_Call(__pyx_t_254, __pyx_t_6, NULL); if (unlikely(!__pyx_t_14)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_14);
__Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
}
}
__Pyx_DECREF(__pyx_t_254); __pyx_t_254 = 0;
__pyx_t_6 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_6);
__pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_6, __pyx_n_s_asarray); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
__pyx_t_6 = __pyx_memoryview_fromslice(__pyx_v_VPerp, 2, (PyObject *(*)(char *)) __pyx_memview_get_double, (int (*)(char *, PyObject *)) __pyx_memview_set_double, 0);; if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_6);
__pyx_t_2 = NULL;
if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_3))) {
__pyx_t_2 = PyMethod_GET_SELF(__pyx_t_3);
if (likely(__pyx_t_2)) {
PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_3);
__Pyx_INCREF(__pyx_t_2);
__Pyx_INCREF(function);
__Pyx_DECREF_SET(__pyx_t_3, function);
}
}
if (!__pyx_t_2) {
__pyx_t_254 = __Pyx_PyObject_CallOneArg(__pyx_t_3, __pyx_t_6); if (unlikely(!__pyx_t_254)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
__Pyx_GOTREF(__pyx_t_254);
} else {
#if CYTHON_FAST_PYCALL
if (PyFunction_Check(__pyx_t_3)) {
PyObject *__pyx_temp[2] = {__pyx_t_2, __pyx_t_6};
__pyx_t_254 = __Pyx_PyFunction_FastCall(__pyx_t_3, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_254)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
__Pyx_GOTREF(__pyx_t_254);
__Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
} else
#endif
#if CYTHON_FAST_PYCCALL
if (__Pyx_PyFastCFunction_Check(__pyx_t_3)) {
PyObject *__pyx_temp[2] = {__pyx_t_2, __pyx_t_6};
__pyx_t_254 = __Pyx_PyCFunction_FastCall(__pyx_t_3, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_254)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
__Pyx_GOTREF(__pyx_t_254);
__Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
} else
#endif
{
__pyx_t_4 = PyTuple_New(1+1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__Pyx_GIVEREF(__pyx_t_2); PyTuple_SET_ITEM(__pyx_t_4, 0, __pyx_t_2); __pyx_t_2 = NULL;
__Pyx_GIVEREF(__pyx_t_6);
PyTuple_SET_ITEM(__pyx_t_4, 0+1, __pyx_t_6);
__pyx_t_6 = 0;
__pyx_t_254 = __Pyx_PyObject_Call(__pyx_t_3, __pyx_t_4, NULL); if (unlikely(!__pyx_t_254)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_254);
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
}
}
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__pyx_t_4 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__pyx_t_6 = __Pyx_PyObject_GetAttrStr(__pyx_t_4, __pyx_n_s_asarray); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_6);
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
__pyx_t_4 = __pyx_memoryview_fromslice(__pyx_v_indOut, 1, (PyObject *(*)(char *)) __pyx_memview_get_double, (int (*)(char *, PyObject *)) __pyx_memview_set_double, 0);; if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__pyx_t_2 = NULL;
if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_6))) {
__pyx_t_2 = PyMethod_GET_SELF(__pyx_t_6);
if (likely(__pyx_t_2)) {
PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_6);
__Pyx_INCREF(__pyx_t_2);
__Pyx_INCREF(function);
__Pyx_DECREF_SET(__pyx_t_6, function);
}
}
if (!__pyx_t_2) {
__pyx_t_3 = __Pyx_PyObject_CallOneArg(__pyx_t_6, __pyx_t_4); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
__Pyx_GOTREF(__pyx_t_3);
} else {
#if CYTHON_FAST_PYCALL
if (PyFunction_Check(__pyx_t_6)) {
PyObject *__pyx_temp[2] = {__pyx_t_2, __pyx_t_4};
__pyx_t_3 = __Pyx_PyFunction_FastCall(__pyx_t_6, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
__Pyx_GOTREF(__pyx_t_3);
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
} else
#endif
#if CYTHON_FAST_PYCCALL
if (__Pyx_PyFastCFunction_Check(__pyx_t_6)) {
PyObject *__pyx_temp[2] = {__pyx_t_2, __pyx_t_4};
__pyx_t_3 = __Pyx_PyCFunction_FastCall(__pyx_t_6, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
__Pyx_GOTREF(__pyx_t_3);
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
} else
#endif
{
__pyx_t_16 = PyTuple_New(1+1); if (unlikely(!__pyx_t_16)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_16);
__Pyx_GIVEREF(__pyx_t_2); PyTuple_SET_ITEM(__pyx_t_16, 0, __pyx_t_2); __pyx_t_2 = NULL;
__Pyx_GIVEREF(__pyx_t_4);
PyTuple_SET_ITEM(__pyx_t_16, 0+1, __pyx_t_4);
__pyx_t_4 = 0;
__pyx_t_3 = __Pyx_PyObject_Call(__pyx_t_6, __pyx_t_16, NULL); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__Pyx_DECREF(__pyx_t_16); __pyx_t_16 = 0;
}
}
__Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
__pyx_t_6 = PyTuple_New(4); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 218, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_6);
__Pyx_GIVEREF(__pyx_t_255);
PyTuple_SET_ITEM(__pyx_t_6, 0, __pyx_t_255);
__Pyx_GIVEREF(__pyx_t_14);
PyTuple_SET_ITEM(__pyx_t_6, 1, __pyx_t_14);
__Pyx_GIVEREF(__pyx_t_254);
PyTuple_SET_ITEM(__pyx_t_6, 2, __pyx_t_254);
__Pyx_GIVEREF(__pyx_t_3);
PyTuple_SET_ITEM(__pyx_t_6, 3, __pyx_t_3);
__pyx_t_255 = 0;
__pyx_t_14 = 0;
__pyx_t_254 = 0;
__pyx_t_3 = 0;
__pyx_r = __pyx_t_6;
__pyx_t_6 = 0;
goto __pyx_L0;