getting started with cython by using cdef
Question: I dont understand how to setup and run code with cython.

I added `cdef`, `double`, etc to pertinent pieces of my code.

setup.py of course the name hello isn't being used. [cython
doc](http://docs.cython.org/src/quickstart/build.html)

    
    
    from distutils.core import setup
    from distutils.extension import Extension
    from Cython.Distutils import build_ext
    
    ext_modules = [Extension("hello", ["hello.pyx"])]
    
    setup(
      name = 'Hello world app',
      cmdclass = {'build_ext': build_ext},
      ext_modules = ext_modules
    )
    

**Edit**

Redefining the question to make it more of a concrete example.

I want to run this system of ODEs through Cython. I am not sure if `double` is
the best choice but the numbers are on the order of magnitudes of 10 to the 10
and are non-integer.

So I want to call this in a bigger code where `mus` is a variable defined in
the script calling the module. I have my `setup.py` file to compile the pyx
file. I am not sure with what I need to do so I can call this ode now. Say I
name the module 3bodyproblem. I would then call it in the scipt as `import
3bodyproblem` and then do `3bodyproblem.3bodyproblem(what would this input
be)'

I have be reading [intro to cython for
odes](http://hplgit.github.io/teamods/cyode/cyode-sphinx/main_cyode.html) but
I am not sure how to use their example with mine. Also, if it needs to be in
rk format, see the code below the first code.

**Code 1**

    
    
    cdef deriv(double u, dt):
        cdef double u[0], u[1], u[2]                                                 
        return [u[3],
                u[4],
                u[5],
                -mus * u[0] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5,
                -mus * u[1] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5,
                -mus * u[2] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5]
    
    
    dt = np.linspace(0.0, t12sec, t12sec)
    u = odeint(deriv, u0, dt, atol = 1e-13, rtol = 1e-13)
    x, y, z, vx, vy, vz = u.T
    

**Code 2**

    
    
    cdef deriv(dt, double u):
            cdef double u[0], u[1], u[2], u[3], u[4], u[5]
            return [u[3],
                    u[4],
                    u[5],
                    -mus * u[0] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5,
                    -mus * u[1] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5,
                    -mus * u[2] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5]
    
    solver = ode(deriv).set_integrator('dop853')
    solver.set_initial_value(u0)
    z = np.zeros(300000)
    
    
    for ii in range(300000):
        z[ii] = solver.integrate(dt[ii])[0]
    
    
    x, y, z, x2, y2, z2 = z.T
    

If I just try to compile code 1, cython doesn't know about the variables
defined the larger `.py` file. I want to avoid setting variables in the cython
code so it can be used else where with out re compiling. Here are the errors:

    
    
    Error compiling Cython file:
    ------------------------------------------------------------
    ...
    
    import numpy as np
    from scipy.integrate import odeint
    
    cdef deriv(double u, dt):
        cdef double u[0], u[1], u[2]                                                 
                    ^
    ------------------------------------------------------------
    
    ODEcython.pyx:7:17: 'u' redeclared 
    
    Error compiling Cython file:
    ------------------------------------------------------------
    ...
    
    import numpy as np
    from scipy.integrate import odeint
    
    cdef deriv(double u, dt):
        cdef double u[0], u[1], u[2]                                                 
                          ^
    ------------------------------------------------------------
    
    ODEcython.pyx:7:23: 'u' redeclared 
    
    Error compiling Cython file:
    ------------------------------------------------------------
    ...
    
    import numpy as np
    from scipy.integrate import odeint
    
    cdef deriv(double u, dt):
        cdef double u[0], u[1], u[2]                                                 
                                ^
    ------------------------------------------------------------
    
    ODEcython.pyx:7:29: 'u' redeclared 
    
    Error compiling Cython file:
    ------------------------------------------------------------
    ...
                -mus * u[0] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5,
                -mus * u[1] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5,
                -mus * u[2] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5]
    
    
    dt = np.linspace(0.0, t12sec, t12sec)
                               ^
    ------------------------------------------------------------
    
    ODEcython.pyx:16:28: undeclared name not builtin: t12sec
    
    Error compiling Cython file:
    ------------------------------------------------------------
    ...
                -mus * u[1] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5,
                -mus * u[2] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5]
    
    
    dt = np.linspace(0.0, t12sec, t12sec)
    u = odeint(deriv, u0, dt, atol = 1e-13, rtol = 1e-13)
                   ^
    ------------------------------------------------------------
    
    ODEcython.pyx:17:16: Cannot convert 'object (double, object)' to Python object
    
    Error compiling Cython file:
    ------------------------------------------------------------
    ...
                -mus * u[1] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5,
                -mus * u[2] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5]
    
    
    dt = np.linspace(0.0, t12sec, t12sec)
    u = odeint(deriv, u0, dt, atol = 1e-13, rtol = 1e-13)
                       ^
    ------------------------------------------------------------
    
    ODEcython.pyx:17:20: undeclared name not builtin: u0
    
    Error compiling Cython file:
    ------------------------------------------------------------
    ...
    cdef deriv(double u, dt):
        cdef double u[0], u[1], u[2]                                                 
        return [u[3],
                u[4],
                u[5],
                -mus * u[0] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5,
                    ^
    ------------------------------------------------------------
    
    ODEcython.pyx:11:17: undeclared name not builtin: mus
    building 'ODEcython' extension
    creating build
    creating build/temp.linux-x86_64-2.7
    x86_64-linux-gnu-gcc -pthread -fno-strict-aliasing -DNDEBUG -g -fwrapv -O2 -Wall -Wstrict-prototypes -fPIC -I/usr/include/python2.7 -c ODEcython.c -o build/temp.linux-x86_64-2.7/ODEcython.o
    ODEcython.c:1:2: error: #error Do not use this file, it is the result of a failed Cython compilation.
    error: command 'x86_64-linux-gnu-gcc' failed with exit status 1
    

**Maybe this will help**

So I want to use the module in a file like this: (different deriv function but
just ignore that the idea is the same)

    
    
    import numpy as np
    from scipy.integrate import ode
    import pylab
    from mpl_toolkits.mplot3d import Axes3D
    from scipy.optimize import brentq
    from scipy.optimize import fsolve
    
    me = 5.974 * 10 ** 24  #  mass of the earth                                         
    mm = 7.348 * 10 ** 22  #  mass of the moon                                          
    G = 6.67259 * 10 ** -20  #  gravitational parameter                                 
    re = 6378.0  #  radius of the earth in km                                           
    rm = 1737.0  #  radius of the moon in km                                            
    r12 = 384400.0  #  distance between the CoM of the earth and moon                   
    rs = 66100.0  #  distance to the moon SOI                                           
    Lambda = np.pi / 6  #  angle at arrival to SOI                                      
    M = me + mm
    d = 300  #  distance the spacecraft is above the Earth                              
    pi1 = me / M
    pi2 = mm / M
    mue = 398600.0  #  gravitational parameter of earth km^3/sec^2                      
    mum = G * mm  #  grav param of the moon                                             
    mu = mue + mum
    omega = np.sqrt(mu / r12 ** 3)
    #  distance from the earth to Lambda on the SOI                                     
    r1 = np.sqrt(r12 ** 2 + rs ** 2 - 2 * r12 * rs * np.cos(Lambda))
    vbo = 10.85  #  velocity at burnout                                                 
    h = (re + d) * vbo  #  angular momentum                                             
    energy = vbo ** 2 / 2 - mue / (re + d)  #  energy                                   
    v1 = np.sqrt(2.0 * (energy + mue / r1))  #  refer to the close up of moon diagram   
    #  refer to diagram for angles                                                      
    theta1 = np.arccos(h / (r1 * v1))
    phi1 = np.arcsin(rs * np.sin(Lambda) / r1)
    #                                                                                   
    p = h ** 2 / mue  #  semi-latus rectum                                              
    a = -mue / (2 * energy)  #  semi-major axis                                         
    eccen = np.sqrt(1 - p / a)  #  eccentricity 
    
    nu0 = 0
    nu1 = np.arccos((p - r1) / (eccen * r1))
    
    
    #  Solving for the eccentric anomaly                                                
    def f(E0):
        return np.tan(E0 / 2) - np.sqrt((1 - eccen) / (1 + eccen)) * np.tan(0)
    
    
    E0 = brentq(f, 0, 5)
    
    
    def g(E1):
        return np.tan(E1 / 2) - np.sqrt((1 - eccen) / (1 + eccen)) * np.tan(nu1 / 2)
    
    
    E1 = fsolve(g, 0)
    
    
    #  Time of flight from r0 to SOI                                                    
    deltat = (np.sqrt(a ** 3 / mue) * (E1 - eccen * np.sin(E1)
                                       - (E0 - eccen * np.sin(E0))))
    
    
    
    #  Solve for the initial phase angle 
    def s(phi0):
        return phi0 + deltat * 2 * np.pi / (27.32 * 86400) + phi1 - nu1
    
    
    phi0 = fsolve(s, 0)
    
    
    nu = -phi0                                                                   
    
    gamma = 0 * np.pi / 180  #  angle in radians of the flight path                     
    
    vx = vbo * (np.sin(gamma) * np.cos(nu) - np.cos(gamma) * np.sin(nu))
    #  velocity of the bo in the x direction                                            
    vy = vbo * (np.sin(gamma) * np.sin(nu) + np.cos(gamma) * np.cos(nu))
    #  velocity of the bo in the y direction                                            
    
    xrel = (re + 300.0) * np.cos(nu) - pi2 * r12
    
    
    yrel = (re + 300.0) * np.sin(nu)
    
    
    u0 = [xrel, yrel, 0, vx, vy, 0]
    
    
    def deriv(u, dt):
        return [u[3],  #  dotu[0] = u[3]                                                
                u[4],  #  dotu[1] = u[4]                                                
                u[5],  #  dotu[2] = u[5]                                                
                (2 * omega * u[4] + omega ** 2 * u[0] - mue * (u[0] + pi2 * r12) /
                 np.sqrt(((u[0] + pi2 * r12) ** 2 + u[1] ** 2) ** 3) - mum *
                 (u[0] - pi1 * r12) /
                 np.sqrt(((u[0] - pi1 * r12) ** 2 + u[1] ** 2) ** 3)),
                #  dotu[3] = that                                                       
                (-2 * omega * u[3] + omega ** 2 * u[1] - mue * u[1] /
                 np.sqrt(((u[0] + pi2 * r12) ** 2 + u[1] ** 2) ** 3) - mum * u[1] /
                 np.sqrt(((u[0] - pi1 * r12) ** 2 + u[1] ** 2) ** 3)),
                #  dotu[4] = that                                                       
                0]  #  dotu[5] = 0                                                      
    
    
    dt = np.linspace(0.0, 259200.0, 259200.0)  #  secs to run the simulation            
    u = odeint(deriv, u0, dt)
    x, y, z, x2, y2, z2 = u.T
    
    
    fig = pylab.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.plot(x, y, z, color = 'r')
    #  adding the moon                                                                  
    phi = np.linspace(0, 2 * np.pi, 100)
    theta = np.linspace(0, np.pi, 100)
    xm = rm * np.outer(np.cos(phi), np.sin(theta)) + r12 - pi2 * r12
    ym = rm * np.outer(np.sin(phi), np.sin(theta))
    zm = rm * np.outer(np.ones(np.size(phi)), np.cos(theta))
    ax.plot_surface(xm, ym, zm, color = '#696969', linewidth = 0)
    ax.auto_scale_xyz([-8000, 385000], [-8000, 385000], [-8000, 385000])
    #  adding the earth                                                                 
    xe = re * np.outer(np.cos(phi), np.sin(theta)) - pi2 * r12
    ye = re * np.outer(np.sin(phi), np.sin(theta))
    ze = re * np.outer(np.ones(np.size(phi)), np.cos(theta))
    ax.plot_surface(xe, ye, ze, color = '#4169E1', linewidth = 0)
    ax.auto_scale_xyz([-8000, 385000], [-8000, 385000], [-8000, 385000])
    
    pylab.show()
    


Answer: try something like this:

    
    
    # hello.pyx
    cimport numpy as np
    import numpy as np
    
    def deriv(np.ndarray u, double mus):
        return [u[3],
                u[4],
                u[5],
                -mus * u[0] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5,
                -mus * u[1] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5,
                -mus * u[2] / (u[0] ** 2 + u[1] ** 2 + u[2] ** 2) ** 1.5]
    

build as before, than use it in python like this:

    
    
    import hello
    [...]
    u_new = hello.deriv(u, mus)  # or put it into odeint
    

`cdef` functions are only callable from the c-side, also you must put all
parameters in the function definition.

